Ambiguity-Fixing in Frequency-Varying Carrier Phase Measurements: Global Navigation Satellite System and Terrestrial Examples

  • NAVIGATION: Journal of the Institute of Navigation
  • June 2023,
  • 70
  • (2)
  • navi.580;
  • DOI: https://doi.org/10.33012/navi.580

Abstract

Carrier phase signals are considered among the key observations in global navigation satellite systems (GNSSs) and several other high-precision interferometric measurement systems. However, these ultra-precise measurements are not fully exploited when the integerness of their inherent ambiguities is discarded during the estimation process. Provided that the integer-estimable functions of their phase ambiguities are properly identified, integer ambiguity resolution (IAR) can be utilized to benefit their parameter solutions. For the GNSS code division multiple access systems with transmitters that broadcast carrier phase signals on identical frequencies, these integer-estimable functions have been characterized and are well-known as double differenced ambiguities. However, this is not the case with “frequency-varying” carrier phase signals that are broadcast by GLONASS satellites, Low-Earth-Orbiting communication satellites, or cellular long-term evolution (LTE) transmitters. This study aims to present full-rank models that can be used to identify integer-estimable ambiguity functions, thereby bringing the observation equations of frequency-varying carrier phase measurements into an IAR-applicable form. Our analytical results are supported by several numerical examples, including GNSS and terrestrial-based IAR as well as a new set of “inter-frequency” integer ambiguities that this study discovers in Galileo multi-frequency carrier phase signals.

Keywords

1 INTRODUCTION

Carrier phase signals have served as key observations in various high-precision measurement systems and remote-sensing applications, including global navigation satellite systems (GNSSs) (Teunissen & Montenbruck, 2017), very long baseline interferometry (VLBI) (Hobiger et al., 2009), interferometric synthetic aperture radar (Kampes & Hanssen, 2004), underwater acoustic positioning (Viegas & Cunha, 2007), radio interferometric positioning systems (RIPS) (Maróti et al., 2005), and opportunistic navigation with non-conventional sensors (Shamaei & Kassas, 2019). The provision of GNSS carrier phase signals can lead to millimeter-level positioning parameter solutions. GNSS carrier phase measurements have also been shown to improve the precision of several other important model parameters including those associated with atmospheric sounding, instrumental calibration, and time transfer (Khodabandeh & Teunissen, 2018).

The central role played by carrier phase measurements is particularly pronounced in cases in which the integerness of phase ambiguous cycles (i.e., the so-called ambiguities) can be fully exploited in the estimation process. Integer ambiguity resolution (IAR) is therefore frequently applied to resolve these phase ambiguities. Ideally, one aims to resolve as many ambiguities as possible in the measurement model. However, because the underlying system of observation equations lacks information content, only certain functions of the ambiguities can be considered integer-estimable. If the integer-estimability condition is not met, IAR may fail to deliver a correct solution (Teunissen, 2019). For GNSS code division multiple access (CDMA) systems whose transmitters broadcast carrier phase signals on identical frequencies, such integer-estimable functions are the well-known double differenced (DD) ambiguities (Khodabandeh & Teunissen, 2019). This is not the case with multi-carrier (or frequency-varying) phase signals, such as those broadcast by GLONASS satellites, Low-Earth-Orbiting (LEO) communication satellites (Khalife et al., 2020), or cellular long-term evolution (LTE) transmitters (Shamaei, 2020). Since the corresponding wavelengths of these signals are transmitter-dependent, the classical double-differencing technique generates non-integer combinations of the phase ambiguities which clearly cannot serve as input for IAR. Therefore, the goal of the present study is to present full-rank models that bring the observation equations of frequency-varying phase signals into a form that is amenable for IAR.

In our recent contribution (Teunissen & Khodabandeh, 2022), we employed the integer-estimability theory developed by Teunissen (2019) to investigate the (in) ability to enable single-receiver IAR for frequency-varying measurement systems in the context of precise point positioning (PPP)-real-time kinematics (RTK). This entails the construction of PPP-RTK corrections that avoid fixing the integer ambiguities to non-integer values. In this study, we plan to showcase yet another feature of integer-estimability theory, namely, the construction of IAR-applicable full-rank models that can be used in frequency-varying measurement systems. The interpretation of the estimable parameters involved in the these specific full-rank models is constructed in a step-by-step manner, and the ambiguity-resolution performance of a few leading examples is evaluated numerically.

Accordingly, we re-parameterize mixed-integer “rank-deficient” models into their “full-rank” integer-estimable parameterized versions. These full-rank models will, in turn, deliver integer-estimable functions of the original ambiguities (Teunissen, 1996; Teunissen & Odijk, 2003). Float solutions of the integer-estimable functions can then serve as input for IAR methods. GNSS applications of integer-estimability theory can be found in the literature (Hou et al., 2020; Khodabandeh & Teunissen, 2019; Teunissen & Khodabandeh, 2019; Zaminpardaz et al., 2021) or (Brack et al., 2021). The applicability of this theory is also illustrated when applied to other frequency-varying carrier phase measurement systems such as those involved in terrestrial LTE networks. This theory may also be used to discover integer-estimable inter-frequency ambiguities in Galileo carrier phase measurements.

The remaining sections of this paper are organized as follows. In Section 2, we will discuss the undifferenced observation equations associated with frequency-varying signals. Highlighting the underlying rank-deficiency of this system, we will first present full-rank models that have been parameterized in terms of real-valued ambiguities. The disadvantage of using such models is that the “integerness” of the carrier phase ambiguities will be discarded. Therefore, in Section 3, we show that the estimable forms of network phase biases can be separated from certain functions of the ambiguities via the use of admissible integer transformations. The float solutions of these integer-estimable ambiguity functions can serve as input for IAR methods, thereby documenting the ambiguity-fixing of frequency-varying carrier phase signals. To evaluate the corresponding integer-estimable parameterized model at work, we first present the results of a simulated network with transmitting/receiving cellular LTE signals in Section 4. Our findings reveal how IAR applied to an LTE network depends on several contributing factors, including measurement precision, auxiliary height-measurements, and the number of epochs involved. In Section 5 we take the applicability of integer-estimability theory one step further by describing our discovery of a new form of “inter-frequency” integer ambiguity detected in Galileo multi-frequency carrier phase signals. Finally, concluding remarks will be provided in Section 6.

2 FULL-RANK MODELS OF FREQUENCY-VARYING SIGNALS

As a point of departure, let sensor r (r = 1,…, n) be a GNSS receiver or a software-defined radio that tracks radio-frequency signals from transmitter s (s = 1,…, m) on frequency-band j (j = 1,…, f). The structure of the corresponding linearized undifferenced carrier-phase and pseudorange (code) observation equations (Leick et al., 2015; Teunissen & Montenbruck, 2017) can be expressed as

1 1

with Graphic and Graphic representing the “observed-minus-computed” phase and pseudorange (code) observables, respectively, and Graphic, the known transmitter-dependent wavelength of frequency Graphic given c as the speed of light in a vacuum. The unknown receiver and transmitter clock offsets (in units of range) are denoted by dtr and dts, respectively. Assuming that the transmitter position is known, the v × 1 vector Δxr contains the unknown receiver position increment vector (of size 3), and, optionally, the wet component of the Earth zenith tropospheric delay (ZTD). Accordingly, vector Δxr is linked to the observables by the 1 × v vector Graphic containing the receiver-to-transmitter line-of-sight unit vector and/or the corresponding ZTD mapping function. Thus, the scalar v is set to v = 4 (when both the position and ZTD are assumed unknown) or to v = 3 (when the ZTD is assumed to be known). The dry component of the tropospheric delay is assumed a priori to have been corrected. The Earth first-order slant ionospheric delays, which are experienced on the L-band frequency Graphic, are denoted by Graphic. While the code observables are biased by the receiver and transmitter code delays dr,j and Graphic, respectively, the phase observables are biased by the real-valued ambiguity term

2 2

in which the integer-valued ambiguities Graphic are accompanied by the non-integer receiver and transmitter phase delays δr,j and Graphic, respectively. Here and in the following sections, potential code-specific mis-modeled effects such as inter-channel biases are assumed absent or calibrated a priori (Aggrey & Bisnath, 2016; Banville et al., 2018; Henkel et al., 2016; Sleewaegen et al., 2012; Wanninger & Wallstab-Freitag, 2007). Working with equations that include undifferenced observations shown in Equation (1) implies that one has to account for rank deficiencies because all unknown parameters cannot be estimated in an unbiased manner. This is because—as described below—one set of parameters can be fully absorbed by those that remain. This set of parameters forms the S-ingularity-basis (or S-basis), and therefore would not be estimable under the system of equations shown in Equation (1) (Baarda, 1973; Koch, 1999; Teunissen, 1985). On the other hand, the remaining estimable parameters do not represent the original parameters as they are biased by the S-basis parameters (Odijk et al., 2015). This concept will be explained with greater clarity in the following section.

2.1 Base Frequencies and their Rational Ratios

To demonstrate explicitly how the stated S-basis parameters are absorbed by the remaining unknown parameters in Equation (1), we will first characterize the links between the underlying carrier frequencies. We assume that the transmitter frequencies Graphic are rational multiples of the base frequency f0; put another way

3 3

Thus, frequency-specific scalars ωj are rational numbers, whereas the transmitter-specific ratios rs are integers. Note, if the ratios rs are not integers, appropriate scaling of f0 and the relative frequency ratios will recover as shown in Equation (3). A prime example of this phenomenon as shown in Equation (3) is the set of frequencies under which GLONASS frequency division multiple access (FDMA) operates. The corresponding ratios are currently given by

4 4

with a base frequency of f0 = 9/16 MHz. Likewise, the GPS CDMA frequencies follow from Equation (3) as

5 5

with a base frequency f0 = 1575.42 MHz.

Substitution of Equation (3) into Graphic gives Graphic with λ0 = c/f0. Taking into account the ambiguity term from Equation (2), this implies that

6 6

Similarly, because the first-order slant ionospheric delays Graphic are inversely proportional to the term frequency-squared Graphic, they can be linked to their GPS counterparts, that are experienced on L1 of fL1 = 1575.42 MHz as follows

7 7

In this way, instead of f frequency-dependent parameters, Graphic, only one ionospheric parameter per receiver-transmitter pair, i.e. Graphic, will need to be estimated.

2.2 Full-Rank Models for Medium-to-Large Scale Networks

Given the parameterization shown in Equation (7), we are in a position to show how the code biases dr,j and Graphic on the first two frequency-bands (j = 1,2) can be lumped with the clock terms dtr and dts and the ionospheric parameter Graphic. Our claim follows from the decomposition of the stated code biases into their ionosphere-free (IF) and geometry-free (GF) combinations (Khodabandeh & Teunissen, 2015)

8 8

The definitions of the IF and GF combinations, together with the symbols Graphic and Graphic, are shown in Table 1. The decomposition defined in Equation (8) indicates that the IF combinations of the code biases, i.e., dr,IF and Graphic, have the same coefficient (i.e., 1) as those of the clock terms dtr and dts. Likewise, the GF combinations dr,GF and Graphic have the same coefficient μj(j = 1,2) as that of the scaled ionospheric parameter Graphic. Thus, the IF combinations of the code biases are absorbed by the clock terms, whereas their GF counterparts are absorbed by the ionospheric parameter. Since these combinations are defined on the basis of the first two frequency-bands j = 1,2, the code biases on the remaining frequency-bands (j > 2), symbolized by Graphic and Graphic, remain estimable. Here and in the following, the Graphic symbol distinguishes the parameters that are estimable from their original counterparts. Likewise, the original clock terms dtrdts and ionospheric parameter Graphic are replaced by their estimable counterparts Graphic and Graphic, respectively. Therefore, the code observation equations shown as Equation (1) would remain unchanged if the original clock and ionospheric parameters were replaced by their estimable versions.

View this table:
TABLE 1

A Medium-to-Large scale network model Shown in Equation (10): Estimable Parameters for a Choice of S-basis

We now consider the phase observation equations shown in Equation (1). These equations will also remain unchanged upon the replacement of Graphic, dtr and dts by their estimable versions if the code-bias term (given inside the box) in

9 9

is absorbed by the non-integer phase biases of the ambiguity vector Graphic, thereby forming the estimable ambiguity vector Graphic. In this way, we arrive at the full-rank version of Equation (1) as follows:

10 10

The interpretations of the estimable parameters involved in Equation (10) and the corresponding S-basis are shown in Table 1. Provided that a minimum number of transmitters (m ≥ 1 + v) are tracked by the network receivers, the receiver position vectors can be estimated in an unbiased fashion, that is Graphic. However, this is not the case with, for example, the receiver clock offsets, i.e., Graphic. In addition to the IF code-bias combinations, the clock offset of the first receiver dt1 is chosen as the S-basis to remove the rank-deficiency between the clock terms dtrdts. Thus, here the estimable clocks are to be determined relative to the time sensed by the first receiver r = 1. The estimable clock offset of the first receiver is then absent or (equivalently) set to zero (i.e., Graphic). Similarly, the estimable code biases on the first two frequency-bands are also absent, i.e., Graphic. It should be noted that, in the formulation shown in Equation (10), the phase biases δr,j and Graphic on all frequencies are chosen as S-basis and therefore included toegether with the ambiguities. This is why the estimable ambiguities Graphic are real-valued. In Section 3 we will show how the integerness of the phase ambiguities can be retrieved at the expense of replacing the S-basis parameters δr,j and Graphic by functions associated with the original ambiguities Graphic.

2.3 Full-Rank Models for Small Scale Networks

The full-rank model shown in Equation (10) applies to medium-to-large scale networks of receivers, e.g., those with inter-station distances larger than 50 km. For small scale networks with inter-station distances below 50 km, one may impose additional constraints on the ionospheric parameters Graphic; see for example Teunissen & Khodabandeh (2021). For instance, the ionospheric delays corresponding to satellite s experienced by all network receivers can be assumed to be identical albeit with some degree of uncertainty, that is, Graphic. Based on this assumption, the estimable ionospheric parameters can be expressed as (cf. Table 1)

11 11

Thus, upon imposing the constraints Graphic, the m times n unknowns Graphic are replaced by m + (n – 1) unknowns, i.e. Graphic (of size m) and Graphic (of size n – 1). Substitution of Equation (11) into Equation (10) provides a full-rank model for small scale networks as follows

12 12

with the estimable code-bias term Graphic. When using the formulation shown in Equation (12), the code-bias term Graphic will also be present in the phase observation equations. To avoid having any estimable code bias present in the phase observation equations, one can re-parameterize Graphic via the following one-to-one transformation:

13 13

Accordingly, the receiver clocks and ambiguities are reformulated as Graphic and Graphic. On the other hand, the (f – 2) estimable code biases Graphic, together with Graphic, are replaced by (f – 1) newly-defined code biases Graphic. By applying the parameterization shown in Equation (13) to Equation (12), one obtains another equivalent, full-rank model which is also applicable to small scale networks:

14 14

One can compare the full-rank model shown in Equation (14) with that shown in Equation (12). Both have the same number of unknown estimable parameters. However, in contrast to Equation (12), the estimable receiver code biases of Equation (14) are only present in the code observation equations. The interpretations of the estimable parameters of Equation (14) are shown in Table 2. In comparison to their counterparts in Table 1, the estimable parameters shown in Table 2 contain extra estimable receiver code biases on the second frequency-band Graphic. This is due to the additional information Graphic that reduces the number of S-basis parameters in the small scale network model shown in Equation (14). In other words, the receiver code biases on the second frequency-band dr,2 (r ≠ 1) are no longer needed to be maintained as S-basis.

View this table:
TABLE 2

Small Scale Network Model Shown in Equation (14): Estimable Parameters for a Choice of S-basis

It is worth noting that, given the short inter-station distances, the receivers experience near parallel line-of-sight vectors Graphic to the same GNSS satellite, as their elevation angles to the satellite are nearly identical. This is because of the high-altitude orbits of GNSS satellites relative to the receivers’ inter-station distances. Therefore, the positioning parameters Graphic (Table 2) may be poorly estimable. The poor estimability of Graphic may even hold for inter-station distances of approximately 500 km as described by Rocken et al. (2005). To avoid potential rank-deficiency, the position of a reference receiver, for example, Graphic, can be held fixed (Odijk et al., 2015).

Both the full-rank models shown in Equation (10) and Equation (14) have real-valued estimable ambiguity parameters (i.e. Graphic and Graphic), meaning that they do not allow a direct application of IAR for full exploitation of the carrier phase measurements Graphic. In the next section, we will show how certain functions of these ambiguities Graphic and Graphic can be identified that have float solutions that can serve as input for use in IAR methods.

3 INTEGER-ESTIMABILITY OF AMBIGUITIES

Thus far, we have shown how the choice of S-basis can lead to the removal of rank-deficiency, thus forming full-rank models. This includes the non-integer phase biases δr,j and Graphic which have been maintained as S-basis so that real-valued estimable ambiguities Graphic (Table 1) and Graphic (Table 2) can be realized. The interpretations of these estimable ambiguities are identical in structure. To facilitate the presentation, in the following section we use the notation Graphic (with a single tilde symbol) to refer to both Graphic and Graphic. Accordingly, the interpretation of the estimable real-valued ambiguities follows from Tables 1 and 2 as compared with Equation 6

15 15

and Graphic. The task is to retrieve the integer component of Graphic, i.e. Graphic, so that one will be in a position to perform IAR when working with full-rank models shown in Equation (10) and Equation (14). Because the non-integer component Graphic is not directly separable from the integer ambiguity Graphic, one must identify integer-estimable functions of Graphic instead. As shown below, this is equivalent to replacing the S-basis parameters δr,j and Graphic with functions associated with the integer ambiguities Graphic.

To show how integer-estimable functions of the ambiguities Graphic in Equation (15) can be identified, we provide an illustrative GNSS example in Figure 1. The figure presents a simple network in which three receivers (n = 3) track three satellites (m = 3). Each solid line in the figure indicates whether or not satellite s is tracked by receiver r and thus represents a receiver-to-satellite connection. While the number of these connections may reach a maximum value of m × n = 9 when each receiver tracks all the satellites, the total number of the connections (lines) in the figure is actually q = 7. This is because the receivers r = 2 and r = 3 do not track satellites s = 3 or s = 1, respectively. If we define the ambiguity vector Graphic as containing all the undifferenced estimable network ambiguities on frequency-band j, the vectorial form of Equation (15) for Figure 1 will be

16 16

with l = m + n – 1. The q × l matrix R contains either 0 or integer ratios ± rs(s = 1,…, m) as entries. Each row of R corresponds to a receiver-to-satellite connection. To specify the R-matrix, we can consider two cases: (1) CDMA satellite signals versus (2) GLONASS FDMA signals. For the latter, the channel numbers of the three satellites are assumed to be κ1 = 1, κ2 = −4, and κ3 = −7 (cf. Equation 4). Based on this assumption, the R-matrices corresponding to the two cases would read

17 17

FIGURE 1

A simple GNSS network in which three receivers (blue triangles) track three GNSS satellites (black squares). Each solid line indicates tracking of satellite s by receiver r.

Given these R-matrices, one then faces the task of eliminating the rank-deficiency between the integer ambiguity vector zj and the non-integer phase-bias vector b,j as previously described in Equation (16). To do so, one may be inclined to group b,j together with z,j. Rank-deficiency removal will lead to the same full-rank models as in Equation (10) and Equation (14), thereby discarding the integerness of z,j. To identify integer-estimable functions of the ambiguities z,j and thus bring Equation (16) into a form to which IAR can be applied, one needs to group the functions of z,j into the non-integer term b,j instead. To this end, we will parameterize the integer ambiguity vector z,j using the same square invertible matrix Graphic and Graphic) that has the inverse-transpose matrix Graphic. Thus Graphic (identity matrix).

Substitution into Equation (16) results in

18 18

with Graphic and Graphic.

Accordingly, the integer ambiguity vector Graphic is decomposed into two components Graphic and Graphic. As the dimension of Graphic is the same as that of Graphic is the candidate that may be added to the phase-bias vector b,j. However, this requires the corresponding design matrices R and Z2 to be linearly-dependent. Thus, the invertible matrix Z = [Z1, Z2] must satisfy the condition Graphic for some matrix Graphic, or, equivalently, the inverse-transpose matrix Graphic must satisfy

Condition 1: 19 19

Under conditions corresponding to Equation (19), the final expression in Equation (18) can be simplified as the full-rank version of Equation (16), that is

20 20

The system of equations described above is full-rank because, under Equation (19), the columns Z1 and Graphic are linearly-independent. With Equation (20), one can separate the estimable form of the phase ambiguities, i.e., Graphic, from the estimable phase-bias vector Graphic as per frequency-band j. However, it has not yet been determined whether or not the float solutions of these estimable functions can serve as input for IAR methods. Clearly, IAR-inputs are required to be “integer-mean”, i.e., the ambiguity vector Graphic must be an integer. This imposes another constraint on the inverse-transpose matrix Graphic, namely, that its entries also must be integers, i.e. Graphic. The integerness of Graphic is a necessary but not sufficient condition. To illustrate this point, one can suppose that the float solution of the estimable ambiguity vector Graphic, say Graphic, can be obtained via Equation (20). Upon employing an IAR method, the real-valued “float” solution Graphic is mapped to the “fixed” solution Graphic. This implies that an original integer ambiguity vector Graphic must exist that satisfies the relation Graphic. Not all integer matrices Graphic can satisfy this relationship, i.e., Graphic for every Graphic. For instance, consider the two-dimensional example Graphic. If the fixed solution happens to be an odd number (e.g., Graphic), no integer vector Graphic would exist that satisfies the equality Graphic. Only one specific class of integer matrices Graphic can satisfy Graphic for every Graphic, which is the class of admissible integer transformations (Teunissen, 1995a). An admissible integer transformation is an integer matrix whose inverse is also integer. Thus, the inverse-transpose matrix Graphic, next to the condition shown in Equation (19), must also satisfy

Condition 2: 21 21

with Graphic.

To obtain the admissible transformations Z and Graphic, a simple integer-sweeping algorithm such as that described by Teunissen (2019) can be employed. The algorithm involves several elementary operations, for example, adding or subtracting integer multiples from one column to or from another column, and reordering or sign-changing the columns. Careful execution of such elementary operations has been shown to deliver admissible integer transformations Z and Graphic that satisfy Graphic for a given integer design matrix R. The output of the algorithm, i.e. Z = [Z1, Z2], enables one to form integer-estimable parameterized models, thereby rendering IAR as feasible for “frequency-varying” carrier phase measurements.

If the R-matrices in Equation (17) are given to the algorithm as input, it will return admissible transformations Graphic and Graphic specified in Tables 3 and 4 as output. Let us first consider Table 3 in which the CDMA integer-estimable ambiguities Graphic follow from the entries of Graphic as

22 22

where use is made of between-receiver and satellite-satellite differencing notations (·)ik = (·)k – (·)i and (·)ik = (·)k – (·)i, respectively. As indicated in Equation (22), the integer-estimable ambiguities for GNSS CDMA are the well-known DD ambiguities or linear functions therefrom. However, this is not the case with GLONASS FDMA. According to Table 4, the FDMA integer-estimable ambiguities Graphic follow from the entries of Graphic as

23 23

View this table:
TABLE 3

Admissible Integer Transformations Graphic and Graphic that Satisfy the Two Integer-Estimability Conditions shown in Equations (19) and (21) for the CDMA R-matrix Shown in Equation (17) Corresponding to the GNSS Network Shown in Figure 1.

View this table:
TABLE 4

Admissible Integer Transformations Graphic and Graphic that Satisfy the Two Integer-Estimability Conditions Shown in Equations (19) and (21) for the FDMA R-matix in Equation (17) Corresponding to the GNSS Network shown in Figure 1.

Thus, in this example, the FDMA integer-estimable ambiguities are linear combinations of the between-receiver single-differenced (SD) ambiguities. It should also be remarked that the rather large entries of the admissible transformations in Table 4 do not harm the precision of the float ambiguity solutions, since the de-correlation step of the LAMBDA method always “conditions” the ambiguity variance matrix so that large variances are reduced, as described by Teunissen & Khodabandeh (2019).

In the following section, we will employ integer-estimable parameterization as shown in Equation (20) to generate an LTE network and show how an underlying IAR responds to several different contributing factors.

4 A SIMULATED LTE EXAMPLE

In this section we will analyze the IAR performance of a simulated network of transmitting/receiving cellular LTE signals. An illustration of this network is shown in Figure 2. The network consists of two stationary receivers r = 1, 2 (blue triangles) that track eight stationary transmitters s = 1,…, 8 (black squares). The east-north-Up local coordinates of the receivers/transmitters are also specified in the figure. As the interstation distances of the LTE network are substantially below 50 km, the network observation equations can be characterized by the full-rank model (14). For the sake of simplicity, in the following example we remove the satellite-specific parameters and work with the between-receiver SD formulation (·)1r = (·)r − (·)1. With the vector notations Graphic and Graphic, the vectorial representation of Equation (14) in its SD form is

24 24

with the m × m diagonal matrix Graphic, the m × 3 geometry matrix Graphic, and the m-vector of ones e = [1,…, 1]T. Only a single frequency-band j is assumed to be available, i.e., f = 1. Therefore, the estimable receiver code bias is absent, i.e., Graphic for j = 1. Note that, instead of receiver coordinates Graphic, the baseline Graphic is treated as unknown in the SD model as shown in Equation (24). Without loss of generality, it is thus assumed that receiver r = 2 takes the role of a rover (with unknown location), while the receiver r = 1 serves as reference (with known location).

FIGURE 2

Visualization of a simulated LTE network in which two receivers r = 1,2 (blue triangles) track eight transmitters s = 1,…, 8 (black squares). The solid and dashed lines indicate the connections between each receiver-transmitter pair. The east-north-Up local coordinates of the receivers/transmitters, together with the corresponding carrier frequencies, are specified as shown.

It should also be noted that LTE measurements can be integrated with other carrier phase-based measurements systems that experience atmospheric (ionospheric and tropospheric) delays. Inclusion of these delays as additional parameters in the model may require additional rank-deficiency removal so that the corresponding integer-estimable parameterized model will be structured appropriately. A way to do this was discussed by Teunissen (2019).

The full-rank model shown in Equation (24) can be directly utilized to position and navigate with “frequency-varying” carrier phase observables. For the single-epoch scenario, the carrier phase data Graphic are reserved for the unknown ambiguity vector Graphic; this means that only the code data Δp1r,j contribute to the estimation of the position Graphic. However, for a multi-epoch scenario in which the unknown ambiguity vector Graphic is considered time-constant (i.e., if no cycle-slip occurs), the carrier phase data can contribute to the estimation of Δxr from the second epoch and onward. This was the approach taken by Shamaei & Kassas (2019) to achieve sub-meter navigation solutions using cellular LTE signals. The disadvantage of using this formulation as shown in Equation (24) is that it ignores the integer component of Graphic, i.e. Graphic. In (ibid), it is argued that conventional IAR methods, for example, LAMBDA (Teunissen, 1995b), cannot be used to resolve the integer ambiguities of frequency-varying phase observables. However, this argument has been challenged by integer-estimable parameterization studies shown in Equation (20). Substitution with parameterization, i.e., Graphic, into Equation (24) generates an integer-estimable parameterized model

25 25

The corresponding R-vector contains integer ratios, i.e., R = [r1,…, r8]T. With reference to the frequency-relation shown in Equation (3), we have Graphic, with ωj = 1, the base frequency f0 = 1 MHz, and the integer ratios r1 = 2145, r2 = 739, r3 = 2125, r4 = 1955, r5 = 2125, r6 = 1955, r7 = 2415, and r8 = 2125 (cf. Figure 2). To obtain the admissible transformation Z = [Z1, Z2], one can insert the R-vector as input to an “integer-sweeping” algorithm (Teunissen, 2019). The corresponding results are shown in Table 5.

View this table:
TABLE 5

Admissible Integer Transformations Graphic and Graphic that Satisfy the Two Integer-Estimability Conditions Shown in Equations (19) and (21) for the R-vector in Equation (25) that Corresponds to the LTE Network shown in Figure 2.

With the integer-estimable parameterized model shown in Equation (25), we are in a position to apply IAR to the float solutions of the integer vector Graphic. To this end, we simulate undifferenced normally-distributed measurements using the MATLAB random-number generator randn, by considering the following two cases:

  • Case 1: The standard deviations of the undifferenced code and phase measurements are set to σp = 0.3 [m] and σϕ = 0.01 [cyc.], respectively.

  • Case 2: The standard deviations of the undifferenced code and phase measurements are set to σp = 1.0 [m] and σϕ = 0.03 [cyc.], respectively.

These standard deviations are typical values for code and phase measurements. However, the stated standard-deviations can take on larger values, depending on the receiver/antenna types, e.g., as described by Shamaei (2020). For each case, 5,000 samples are generated. In both cases, the measurements are assumed to be mutually uncorrelated. As shown in Figure 2, the LTE transmitters inherit a rather limited vertical diversity, meaning that the Up-component of receiver r = 2 is poorly estimable. This issue was also discussed by Shamaei & Kassas (2019). In (ibid), use of an “altimeter” that delivers extra height measurements was proposed. In the following, we therefore consider a scenario in which an auxiliary height-measurement with a standard-deviation of one meter is available. “Simulation” was used to study the ambiguity-resolution performance of the model shown in Equation (25) for different contributing factors, including the measurement precision, receiver motion, and the number of epochs involved. Simulation also enables us to conduct the required analysis and validation under controlled circumstances.

To provide insight into the IAR performance of the LTE model shown in Equation (25), we evaluated the corresponding empirical ambiguity success-rate as function of the number of epochs for both Cases 1 and 2 shown in Figure 3. The integer least-squares estimation was employed for fixing the ambiguities. The stated ambiguity success-rate was evaluated as the ratio of correctly-fixed ambiguities within a total of 5,000 samples. Given its rather precise phase and code measurements, Case 1 only required approximately 10 epochs to exhibit ambiguity success-rates >99.5% (dashed lines). By contrast, for Case 2, one needs to collect approximately 50 epochs of data to achieve large ambiguity success-rates (solid lines). Note that these results concern stationary receivers/transmitters whose geometric matrix G shown in Equation (25) does not change over time to help improve the precision of the float ambiguities (Teunissen, 1995b). To show how a change in the receiver-transmitter geometry contributes to the IAR performance, we also consider a scenario in which receiver r = 2 moves in the north-east direction with the constant velocity vector [10, 10, 0] [m/s] (east-north-Up coordinates). This is in contrast to the stationary case in which the unknown baseline Graphic is assumed time-constant. The corresponding empirical ambiguity success-rate for Case 2 is highlighted in gray. The required number of epochs has been reduced from 50 to 30 epochs.

FIGURE 3

Empirical ambiguity success-rate of the LTE network (%) as function of the number of epochs (cf. Figure 2) when receiver r = 2 is stationary (black lines) or moving (gray line). The undifferenced code and phase measurements’ standard deviations are set to σp = 0.3 [m] and σϕ = 0.01 [cyc.] (dashed line) and σp = 1.0 [m] and σϕ = 0.03 [cyc.] (solid lines), respectively.

To show how successful IAR corresponding to the model shown in Equation (25) is reflected in the positioning domain, we distinguish the following two cases:

  • Ambiguity-float: The integerness of the integer-estimable ambiguities Graphic in Equation (25) is discarded, leaving the float solution of Graphic as a real-valued vector;

  • Ambiguity-fixed: The float solutions of the integer-estimable ambiguities Graphic are resolved as integers. The fixed ambiguities are therefore maintained as known, which benefits the solutions of Graphic as documented in Equation (25).

Figure 4 presents the ambiguity-float solution’s standard deviations of the baseline Graphic (left) and the corresponding float-to-fixed standard deviation ratios (right) as function of the number of epochs. As indicated in the left panel of the figure, the precision of the Up-component before IAR is almost equal to 1 meter, i.e., equal to the precision of the auxiliary height-measurement. This implies that the ambiguity-float Up-solution is largely driven by the auxiliary height-measurement, and that the Up-component is poorly estimable in the absence of this auxiliary measurement. On the other hand, the float solutions of the horizontal components improve with a larger the number of epochs.

FIGURE 4

Standard deviations of the of the float solution baseline Graphic as shown in Figure 2 (left) and the corresponding float-to-fixed standard deviation ratios (right) as function of the number of epochs. The code and phase measurement standard deviations (undifferenced) are set to σp = 0.3 [m] and σϕ = 0.01 [cyc.] (dashed lines), and σp = 1.0 [m] and σϕ = 0.03 [cyc.] (solid lines), respectively

The right panel of Figure 4 tells us how many times the precision of the ambiguity-float solutions improves after successful IAR. The solutions of the horizontal components experience two orders of magnitude of improved precision which is commensurate with the phase-to-code precision ratio. However, the Up-solution does not generate similar precision improvement during the initial epochs. This is because of the low vertical diversity of the LTE receivers/transmitters. During the initial epochs, the Up-solution continues to rely on the auxiliary height-measurement. As the number of epochs increases, the ambiguity-resolved carrier phase data will gradually contribute to the solution of the “time-constant” Up-component, thereby increasing the corresponding float-to-fixed standard deviation ratio.

5 GALILEO INTER-FREQUENCY AMBIGUITIES

Up to now, the integer-estimable ambiguities Graphic were formed independently for each individual frequency-band j. This is because the real-valued ambiguities Graphic, detected on each frequency-band j, contain their own frequency-specific non-integer biases br,j and Graphic as documented in Equation (15). In this section, we will consider a scenario in which the SD receiver phase biases are common among the frequency bands (i.e. b1r,j = b1r,1, j ≠ 1), showing how integer-estimability theory lends itself to discovery of integer inter-frequency ambiguities. For this purpose, let us revisit the frequency relation shown in Equation (3). In a manner that is analogous to the integer ratios rs, the rational scalars ωj can be scaled up so that they also become integers. For instance, the Galileo E5a, E5b, and E5 frequencies can be characterized by Equation (3) as follows

26 26

with the base frequency f0 = 5.115 MHz. For the aforementioned Galileo triple-frequency case, the R-vector described in Equation (25) simply reads R = [1,…, 1]T (the m-vector of ones). Accordingly, the integer-sweeping algorithm (Teunissen, 2019) returns the outputs Graphic and Graphic. Thus Graphic. With this in mind, the interpretation of the estimable SD receiver phase biases follows from Graphic and the definition of br,j as described in Equation (15) as (cf. Table 2)

27 27

using between-frequency differencing notation δ1r,1j = δ1,r,jδ1r,1. Compare Equation (27) with Equation (15). If one postulates that the non-integer term δ1r,1j is absent, Equation (27) would then inherit the same structure as that of Equation (15) with the role of rs switched by ωj. Let us now follow these lines of thought and obtain an admissible transformation that satisfies the two conditions set by Equations (19) and (21) for the integer vector R = [ω1, ω2, ω3]T corresponding to the Galileo frequencies in Equation (26). This gives

28 28

Pre-multiplying the vector of estimable SD phase biases Graphic by Graphic gives the following two estimable quantities

29 29

If the conjecture δ1r,1j = 0 (j ≠ 1) is true, the above α- and β-quantities could become integer-estimable with the following interpretation

30 30

Such integer-estimable α- and β-quantities would then be of the “interfrequency” type.

To verify this conjecture, we utilized daily data sets of the Galileo triple-frequency (E5a/E5b/E5) signals collected by three zero-baselines (Table 6) and four very short baselines (Table 7). The lengths of all the short baselines involved were <21 meters and the measurement sampling-rate was 30 seconds. The data sets were collected on the 10th of April 2022. In the case of the zero-baseline CUT0-CUT2, an additional data set was also collected on the 20th of May 2022 to determine how the solutions of Equation (29) behave over time. These data are freely accessible from the Curtin GNSS Research Centre (http://saegnss2.curtin.edu.au/ldc/) and the IGS online archives of the Crustal Dynamics Data Information System (CDDIS) (https://cddis.nasa.gov/archive/gnss/data/).

View this table:
TABLE 6

The Receiver and Antenna Types of the zero-baselines Used in the Galileo Experiment Together with the Corresponding Signals.

View this table:
TABLE 7

The Receiver and Antenna types of the Short-Baselines Used in the Galileo Experiment Together with the Corresponding Signals.

Using the integer-estimable parameterized model shown in Equation (25) with known baseline (i.e. Graphic), we first computed the ambiguity-fixed solutions of the SD phase biases,i.e., Graphic. Then, we pre-multiplied the solution Graphic by the matrix Graphic in Equation (28) to evaluate the α- and β-solutions. The solutions were evaluated on a single-epoch basis, meaning that there are 2,880 individual solutions for each daily data set. These single-epoch solutions (green dots) are shown in Figures 5 (zero-baselines) and 6 (short baselines). To show how the solutions differ from those of the integer vector [0, 0]T, the corresponding integer least-squares pull-in regions (gray and blue lines), as well as their zero-centered 99.9%-confidence ellipses (in red), were also depicted in the figures.

Let us first consider the zero-baseline results. According to Table 6, each baseline is formed by the same receiver type. As indicated by the vertical axes in Figure 5, the mean values of the α-solutions are close to zero with standard deviations less than or equal to 0.0001 cycles. As the magnitudes of these means are statistically significant relative to their standard deviations, one cannot conclude that the α-quantities are integer-valued. Instead, we can conclude that the α-quantities are “close” to integers. Since the α- and β-solutions are correlated, they are LAMBDA-de-correlated by admissible transformations ZD as specified in the figure. Accordingly, the horizontal axes represent the solutions of the combinations Graphic. As shown, the solutions of β + are not zero-mean and are considerably less precise than their α-counterparts.

FIGURE 5

Zero-baseline results: single-epoch solutions (green dots) of the α- and β-quantities shown in Equation (29) using the full-rank model shown in Equation (25). The solutions are LAMBDA-de-correlated by an admissible transformation ZD. Thus, the vertical axes represent α-solutions, whereas the horizontal axes represent the solutions of the combinations β + , where k = −123 (top left), k = −114 (top right), k = −114 (bottom left), and k = −114 (bottom right). Their integer least-squares pull-in regions are indicated by gray and blue lines, whereas their zero-centered 99.9%-confidence ellipses are depicted in red.

We now turn our attention to the short baseline results. With reference to Table 7, the three baselines CUTB-CUTC, YARR-YAR3 and NYAL-NYA1 are formed by the same receiver types. By contrast, this is not the case with the baseline UNB3-UNBD. This finding may explain the rather large deviations of its α- and β-solutions from the integer vector [0, 0]T (see Figure 6). For the first three short baselines, a behavior similar to the zero-baselines is observed. The very precise α-solutions are close to zero, while the less-precise solutions of the combinations β + are not zero-mean. When receivers of identical types are used, we can therefore conclude that the phase-bias linear combination Graphic can represent an integer “inter-frequency” ambiguity that is shifted by a statistically small non-integer bias.

FIGURE 6

Short-baseline results: single-epoch solutions (green dots) of the α- and β-quantities as shown in Equation (29) using the full-rank model shown in Equation (25). The solutions are LAMBDA-de-correlated by an admissible transformation ZD. Thus, the vertical axes represent α-solutions, whereas the horizontal axes represent the solutions of the combinations β + , where k = −119 (top left), k = −15 (top right), k = −189 (bottom left), and k = −119 (bottom right). Their integer least-squares pull-in regions are indicated by gray and blue lines, whereas their zero-centered 99.9%-confidence ellipses are depicted in red.

6 CONCLUDING REMARKS

In this contribution, we introduced full-rank models to identify integer-estimable ambiguity functions and bring the observation equations of frequency-varying carrier phase measurements into a form to which IAR can be applied. We showed how one can form full-rank integer-estimable parameterized models using certain admissible transformations. We applied this concept to a number of examples, ranging from GNSS and terrestrial-based IAR to a new form of “inter-frequency” integer ambiguities that were discovered in Galileo multi-frequency carrier phase signals.

In contrast to the GNSS CDMA where carrier phase signals are broadcast on identical frequencies, the integer-estimable functions of the frequency-varying carrier phase measurements are not necessarily the well-known DD ambiguities. According to the results shown in Table 5, the first component of the integer-estimable ambiguities of the LTE network in Figure 2 reads

31 31

These ambiguities are clearly not of the DD type. By presenting the empirical ambiguity success-rate of the corresponding float solutions, we showed how IAR performance in an LTE network responds to several contributing factors such as measurement precision and the number of epochs involved. Depending on these parameters, successful ambiguity resolution was possible and beneficial and facilitated high-precision parameter estimation even when working with “frequency-varying” carrier phase data.

The applicability of integer-estimability theory was also extended to the discovery of inter-frequency ambiguities that may potentially exist among Galileo multi-frequency phase measurements. Successful resolution of these inter-frequency ambiguities may help to improve both the solutions’ precision-performance and the integrity of the underlying model.

How to cite this article

Khodabandeh, A., & Teunissen, P. J. G. (2023) Ambiguity-fixing in frequency-varying carrier phase measurements: GNSS and terrestrial examples. NAVIGATION, 70(2). https://doi.org/10.33012/navi.580

CONFLICT OF INTEREST

The authors declare no potential conflict of interests.

ACKNOWLEDGMENTS

This work was conducted as part of the first author’s Alexander von Humboldt Fellowship at the Institute for Communication and Navigation, the German Aerospace Center (DLR), Munich, with Professor Christoph Günther and Dr. Gabriele Giorgi as his hosts. The GNSS data used in this study have been made available by the Curtin GNSS Research Centre and the International GNSS Service (IGS) through the online archives of the Crustal Dynamics Data Information System (CDDIS). All support is gratefully acknowledged.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

REFERENCES

  1. Aggrey, J., & Bisnath, S. (2016). Dependence of GLONASS pseudorange inter-frequency bias on receiver–antenna combination and impact on precise point positioning. NAVIGATION, 63(4), 379391. https://doi.org/10.1002/navi.168
  2. Baarda, W. (1973). S-transformations and Criterion Matrices (Tech. Rep.). Netherlands Geodetic Commission, Publ. on Geodesy, New Series, Vol. 5(1), Delft.
  3. Banville, S., Collins, P., & Lahaye, F. (2018). Model comparison for GLONASS RTK with low-cost receivers. GPS Solutions, 22(2), 112. https://doi.org/10.1007/s10291-018-0712-3
  4. Brack, A., Männel, B., & Schuh, H. (2021). GLONASS FDMA data for RTK positioning: a five-system analysis. GPS solutions, 25(1), 113. https://doi.org/10.1007/s10291-020-01043-5
  5. Henkel, P., Mittmann, U., & Iafrancesco, M. (2016). Real-time kinematic positioning with GPS and GLONASS. In 2016 24th European Signal Processing Conference (EUSIPCO) 10631067. https://doi.org/10.1109/EUSIPCO.2016.7760411
  6. Hobiger, T., Sekido, M., Koyama, Y., & Kondo, T. (2009). Integer phase ambiguity estimation in next-generation geodetic very long baseline interferometry. Advances in Space Research, 43(1), 187192. https://doi.org/10.1016/j.asr.2008.06.004
  7. Hou, P., Zhang, B., & Liu, T. (2020). Integer-estimable GLONASS FDMA model as applied to Kalman-filter-based short-to long-baseline RTK positioning. GPS Solutions, 24(4), 114. https://doi.org/10.1007/s10291-020-01008-8
  8. Kampes, B. M., & Hanssen, R. F. (2004). Ambiguity resolution for permanent scatterer interferometry. IEEE Transactions on Geoscience and Remote Sensing, 42(11), 24462453. https://doi.org/10.1109/TGRS.2004.835222
  9. Khalife, J., Neinavaie, M., & Kassas, Z. M. (2020). Navigation with differential carrier phase measurements from megaconstellation LEO satellites. In Proc. of the 2020 IEEE/ION Position Location and Navigation Symposium (PLANS). Portland, OR. 13931404. https://doi.org/10.1109/PLANS46316.2020.9110199
  10. Khodabandeh, A., & Teunissen, P. J. G. (2015). An analytical study of PPP-RTK corrections: precision, correlation and user-impact. Journal of Geodesy, 89(11), 11091132. https://doi.org/10.1007/s00190-015-0838-9
  11. Khodabandeh, A., & Teunissen, P. J. G. (2018). On the impact of GNSS ambiguity resolution: geometry, ionosphere, time and biases. Journal of Geodesy, 92(6), 637658. https://doi.org/10.1007/s00190-017-1084-0
  12. Khodabandeh, A., & Teunissen, P. J. G. (2019). Integer estimability in GNSS networks. Journal of Geodesy, 93(9), 18051819. https://doi.org/10.1007/s00190-019-01282
  13. Koch, K.-R. (1999). Parameter estimation and hypothesis testing in linear models. Springer, Berlin. https://doi.org/10.1007/978-3-662-03976-2
  14. Leick, A., Rapoport, L., & Tatarnikov, D. (2015). GPS Satellite Surveying (Fourth ed.). John Wiley and Sons.
  15. Maróti, M., Völgyesi, P., Dóra, S., Kus`y, B., Nádas, A., Lédeczi, A.,… Molnár, K. (2005). Radio Interferometric Geolocation. In Proc. of the 3rd international conference on embedded networked sensor systems. 112. https://doi.org/10.1145/1098918.1098920
  16. Odijk, D., Zhang, B., Khodabandeh, A., Odolinski, R., & Teunissen, P. J. G. (2015). On the estimability of parameters in undifferenced, uncombined GNSS network and PPPRTK user models by means of S-system theory. Journal of Geodesy, 90(1), 1544. https://doi.org/10.1007/s00190-015-0854-9
  17. Rocken, C., Johnson, J., Hove, T. V., & Iwabuchi, T. (2005). Atmospheric water vapor and geoid measurements in the open ocean with GPS. Geophysical Research Letters, 32(L12813). https://doi.org/10.1029/2005GL022573
  18. Shamaei, K. (2020). Exploiting cellular signals for navigation: 4G to 5G. PhD Thesis, University of California, Irvine.
  19. Shamaei, K., & Kassas, Z. M. (2019). Sub-meter accurate UAV navigation and cycle slip detection with LTE carrier phase measurements. In Proc. of the 32nd International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2019), Miami, FL. 24692479. https://doi.org/10.33012/2019.17051
  20. Sleewaegen, J., Simsky, A., De Wilde, W., Boon, F., & Willems, T. (2012). Origin and compensation of GLONASS inter-frequency carrier phase biases in GNSS receivers. In Proc. of the 25th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2012), Nashville, TN. 29953001.
  21. Teunissen, P. J. G. (1985). Generalized inverses, adjustment, the datum problem and S-transformations. In E. W. Grafarend & F. Sanso (Eds.), Optimization and design of geodetic networks 1155. Springer, Berlin.
  22. Teunissen, P. J. G. (1995a). The invertible GPS ambiguity transformations. Manuscripta Geodaetia, 20, 489497.
  23. Teunissen, P. J. G. (1995b). The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation. Journal of Geodesy, 70(1–2), 6582. https://doi.org/10.1007/BF00863419
  24. Teunissen, P. J. G. (1996). Rank defect integer least-squares with applications to GPS. Bollettino di geodesia e scienze affini, 55(3), 225238.
  25. Teunissen, P. J. G. (2019). A new GLONASS FDMA model. GPS Solutions, 23, 100. https://doi.org/10.1007/s10291-019-0889-0
  26. Teunissen, P. J. G., & Khodabandeh, A. (2019). GLONASS ambiguity resolution. GPS Solutions, 23, 101. https://doi.org/10.1007/s10291-019-0890-7
  27. Teunissen, P. J. G., & Khodabandeh, A. (2021). A mean-squared-error condition for weighting ionospheric delays in GNSS baselines. Journal of Geodesy, 95(11), 112. https://doi.org/10.1007/s00190-021-01569-7
  28. Teunissen, P. J. G., & Khodabandeh, A. (2022). PPP-RTK theory for varying transmitter frequencies with satellite and terrestrial positioning applications. Journal of Geodesy, 96(11), 124. https://doi.org/10.1007/s00190-022-01665-2
  29. Teunissen, P. J. G., & Montenbruck, O. (Eds.). (2017). Springer Handbook of Global Navigation Satellite Systems. Springer. https://doi.org/10.1007/978-3-319-42928-1
  30. Teunissen, P. J. G., & Odijk, D. (2003). Rank-defect integer estimation and phase-only modernized GPS ambiguity resolution. Journal of Geodesy, 76(9), 523535. https://doi.org/10.1007/s00190-002-0285-2
  31. Viegas, D. C. d. N., & Cunha, S. R. (2007). Precise positioning by phase processing of sound waves. IEEE Transactions on Signal Processing, 55(12), 57315738. https://doi.org/10.1109/TSP.2007.900166
  32. Wanninger, L., & Wallstab-Freitag, S. (2007). Combined processing of GPS, GLONASS, and SBAS code phase and carrier phase measurements. In Proc. of the 20th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2007), Fort Worth, TX. 866875.
  33. Zaminpardaz, S., Teunissen, P. J. G., & Khodabandeh, A. (2021). GLONASS only FDMA + CDMA RTK: Performance and outlook. GPS Solutions, 25(96). https://doi.org/10.1007/s10291-021-01132-z
Loading
Loading
Loading
Loading
  • Share
  • Bookmark this Article