Bias-Bounded Estimation of AmbiguiTy: A Method for Radio Interferometric Positioning

Carrier phase signals form the basis of various interferometric measurement models and estimation techniques that have parameters of which some are integer. In this paper, integer least-squares estimation theory is extended and applied to bias-bounded mixed-integer models. This extension accommodates the presence of bounded real-valued parameters in mixed-integer models through incorporating prior knowledge of a set, in which the parameters reside, into the estimation process. This enables one to jointly estimate the ambiguous phase cycles and the parameters of interest. To compute such mixed-integer estimates, a fast search strategy is developed that makes use of a dual-ellipsoid encompassing region. The volume of the stated region is quantified and its links to existing ellipsoidal search spaces are highlighted. Simulated and real-world data are employed to illustrate the theory. It is then, for the first time, shown that the proposed method makes single-epoch, phase-only positioning feasible with Global Navigation Satellite Systems (GNSS).


I. INTRODUCTION
B EAT patterns of the measurable carrier waves form the basis of various high-precision measurement systems and remote-sensing applications. Next to Global Navigation Satellite Systems (GNSS) in which radio-frequency carrier phase signals are processed to deliver ultra precise parameter solutions [1]- [3], several other interferometric techniques, such as Very Long Baseline Interferometry (VLBI) [4], Interferometric Synthetic Aperture Radar (InSAR) [5], and underwater acoustic carrier phase positioning [6], employ carrier beat phase measurements in their estimation process. In recent years, carrier phase-based parameter estimation has also received attention in wireless sensor localization techniques like the Radio Interferometric Positioning Systems (RIPS) [7]- [9]. The phase measurement model of the applications stated above can be cast in the following system of observation equations [10]  where the phase observation vector y ∈ R m is linked to the unknown integer-valued ambiguity vector z and the unknown real-valued parameter vector x through the known full-column rank design matrix A ∈ R m×n . The mixed-integer model (1) is clearly underdetermined as y is fully reserved for the ambiguity vector z. Therefore, the parameters of interest x cannot be determined with the sole use of the information content in (1). To have (1) solvable for x, one therefore has to exploit the integer constraints z ∈ Z m , resolving the phase ambiguities z to their correct integers. The approach often taken is to extend the system (1) by additional observation equations [3]. The provision of additional measurements leads to an initial solution for x, saŷ x ∈ R n . The real-valued float ambiguity solutionẑ = y−Ax would then serve as input to methods of integer ambiguity resolution [10] to carry out the integer-mappingẑ →ž (ž ∈ Z m ), thereby producing the 'ambiguity-resolved' phase observation vector y −ž. Provided that the integer-valued solutionž represents the correct ambiguity vector z, the updated measurement model E(y−ž) = Ax can then be solved for x.
Although the above approach has been routinely used, it is restrictive for the cases where additional measurements and, therefore, the initial solutionx is not available. For instance, while GNSS users enjoy the support of pseudo-range data with time-constant phase ambiguities [3], RIPS and InSAR users are often left with phase measurements only [5], [7]. To compensate the lack of such information, the common procedure followed in practice is to include 'pseudo-observations' from external sources or to impose constraints on the phase ambiguities z. The important consequence of this practice is that the quality of the integer solutionž (and that of x) is driven by the assumption made about the quality of such pseudo-observations. In the event that access to external sources is not provided nor would a proper quality description of such pseudo-observations be specified, the information content in the observables y cannot be fully exploited.
In the absence of such additional information, the question that comes to the fore is whether one can take recourse to prior knowledge about the geometry of the measurement setup to determine the unknown ambiguities z, thus leveraging the ambiguity-resolved measurements y −ž in the estimation of x. In this contribution we therefore aim to develop an estimation method upon which the phase observation vector y is directly mapped tož ∈ Z m on the basis of the information content in (1) and prior knowledge of a set in which the parameter vector x lies. To this end, the existing theory of integer least-squares estimation [11]- [13] is extended so as to accommodate the presence This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ of 'bounded' real-valued parameters in the mixed-integer model (1). Next to a strategy for computing the ambiguity solutionž, necessary and sufficiency conditions are derived to guarantee that every observation vector gets mapped to a unique ambiguity solution.
The organization of this paper is as follows. In Section II, we introduce the bias-bounded mixed-integer model on which our discussion is based. The system of GNSS carrier phase observation equations is briefly reviewed as leading example. Section III is devoted to the new integer estimation method. Its admissibility conditions are discussed and its links to existing integer estimators are highlighted. To address how the proposed integer estimationž ∈ Z m can be searched and computed, we discuss the geometry of its associated search space in Section IV. After quantifying the volume of the stated search space, a fast search strategy for the computation ofž is then presented in Section V. The strategy benefits from the 'canonical' form of the mixed-integer model (1). Finally, concluding remarks are given in Section VI.
A number of examples are presented throughout the text to illustrate the theory, while the proofs of the main results are provided in the appendix. We make use of the following notation: The expectation and dispersion operators are shown by E(·) and D(·), respectively. The m-dimensional spaces of real, rational and integer numbers are, respectively, denoted by R m , Q m , and Z m . Subsets are indicated by calligraphic fonts, e.g., M ⊂ R n or S ⊂ R m . Vectors and matrices are indicated by bold-italic lowercase and uppercase letters, respectively. Thus a ∈ R is a scalar, a ∈ R m is an m-vector, and A ∈ R m×n is an m × n matrix. The transpose of a matrix is indicated by the superscript T , i.e. (·) T . The identity matrix of order m is denoted as I m , and A ⊗ B denotes the Kronecker product of A and B. The notation || · || Q is the weighted norm whose weight matrix is given by the inverse of the positive definite matrix Q. Thus ||a|| 2 Q = a T Q −1 a. When Q = I m , the weighted norm reduces to the standard Euclidean norm || · ||, i.e. ||a|| = ||a|| I m . The gamma and beta functions are given by Γ(x) = ∞ 0 v x−1 exp(−v)dv and β(x, y) = Γ(x)Γ(y)/Γ(x + y), respectively.

II. BIAS-BOUNDED MIXED-INTEGER MODEL
Consider two GNSS receivers that simultaneously track the satellites q and s on frequency-band j. One receiver serves as pivot (with known location), and the other as rover (with unknown location). It is assumed that the distance between the receivers is short enough (e.g., less than 10-20 km) so as to neglect the presence of between-receiver atmospheric delays. With this in mind, the corresponding GNSS differential phase observation equations read [10] where φ qs j , expressed in cycles, denotes the phase measurement of the satellite pair qs on frequency-band j. The unknown phase integer ambiguity is denoted by z qs j ∈ Z. The term ρ qs (u) denotes the unknown double-differenced range which is linked to the measurement φ qs j through the known carrier wavelength λ j . In the context of GNSS observation modelling, the term 'double-differenced' (DD) is used to indicate that ρ qs (u) is formed by two successive differences of the receiver-to-satellite distances: 1) between-satellite and 2) between-receiver differences. Fig. 1 shows a simple visualization of GNSS carrier phase-based positioning, together with the definition of ρ qs (u).
As presented in the figure, ρ qs (u) is a nonlinear function of the rover state-vector u ∈ R ν (ν = 3 if 3D-localization is considered). Given the known state-vectors of the reference receiver and GNSS satellites, the goal of carrier phase-based positioning is to determine the unknown vector u through a phase-driven solution of DD ranges. Note that the system of observation equations (2) is not GNSS-specific, but also forms the basis of other interferometric techniques [10]. For instance, (2) also underlies the measurement setup of the Radio Interferometric Positioning System (RIPS) in wireless sensor networks [7], be it that a different naming is used to refer to the DD range ρ qs (u), the so-called Q-range [8]. Now consider the provision of multiple GNSS carrier phase signals that are transmitted by (n+1) satellites on f frequencies. Thus n differential phase measurements φ qs j (s = 1, . . . , n+1; s = q) per frequency can be formed. To express (2) in its vector form, we define the phase observa- With a likewise definition for the ambiguity vector z ∈ Z m , the vectorial form of the observation equations (2) can then be expressed as The unknown vector ρ(u) ∈ R n , as a vector-function of the state-vector u ∈ R ν (ν ≤ n), contains the DD ranges ρ qs (u) (s = q) that are scaled by the wavelength of the first frequency λ 1 . The design matrix A ∈ R m×n is structured by the f -vector a containing the wavelength-ratios λ 1 /λ j (j = 1, . . . , f). As the number of unknown ambiguities z is as many as the number of observations y, the system of observation equations (3) is not solvable for z and u. To properly cope with this underdetermined system, we first consider the estimability of u as function of an arbitrary integer vector z o ∈ Z m . Letǔ z o be a solution of u under E(y−z o ) = Ax. By allowing z o to take on any integer vectors in Z m , one indeed arrives at a class of solutions for u. Such a class has infinitely many membersǔ z o (z o ∈ Z m ), out of which the vectorǔ z o =z represents the soughtfor solution. The idea is now to incorporate possible values, that the state-vector u can take, into the model, with the intention to rule out the redundant solutionsǔ z o =z . Such possible values follow from the geometry of the measurement setup. For instance, in order to neglect between-receiver atmospheric delays, the distance between the GNSS pivot and rover receivers should not exceed certain values [10]. Moreover, the rover receiver is often known to be within a certain area. Likewise, the inter-node distances of the RIPS-based transceivers cannot be longer than a threshold (e.g. 170 meters) so as to be able to detect the interference phase signals [7]. For the sake of argument, let such values be characterized by the set {u ∈ R ν | ||u − u o || ≤ r o }, in which the state-vector u o and the radius r o are given. If the operator ρ(u) is bounded, the parameter vector x is then forced to reside in a 'proper' subset of R n . This notion is made precise in the following.
Definition 1 (Bias-bounded mixed-integer model): Let A ∈ R m×n be a given matrix of full-column rank, and let Q yy ∈ R m×m be a given positive definite matrix. With M being a given nonempty and bounded set in R n , the following model will be referred to as a bias-bounded mixed-integer model. The nonempty set M ⊂ R n is said to be bounded if for all x ∈ M, there exist x o ∈ M and h ≥ 0 such that ||x−x o || ≤ h. The term 'mixed-integer' indicates that the model (5) is composed of both real-(x) and integer-valued (z) unknown parameters, whereas the term 'bias-bounded' is given to highlight that the bias Ax = E(y) − z is bounded through x. When the bias is absent, i.e. x = 0, y represents an integer-mean random vector, meaning that y can directly serve as input for integer ambiguity resolution. In the presence of bias however, one would need to devise an integer estimator y →ž which can, next to the observation vector y and its variance matrix Q yy , also take the role of A and M into account. Such integer estimator will be discussed in Section III.
The model (5) is versatile and can be applicable to any carrier phase-based interferometric techniques. For the InSAR permanent scatterer technique, the real-valued vector x may (among others) contain 'displacement' parameters of terrestrial points [5], whereas for the VLBI technique, x may contain Earth's tropospheric parameters [4]. Depending on the application at hand, the real-valued vector x can be assumed to reside in an a-priori known bounded set M ⊂ R n . From a geometrical point of view, the bounded set M describes a 'segment' of a (nonlinear) manifold in R n that lies inside an n-ball {x ∈ R n |||x − x o || ≤ h; x o ∈ M}. A GNSS example is presented below to provide further insight into the set M.
Example 1: If the distance between the rover receiver and a certain location u o is known to be not longer than a positive scalar r o , i.e. ||u − u o || ≤ r o , then the set M can be characterized as The set (6) describes a 'segment' of a ν-dimensional manifold in R n . To show that the set is bounded, we make use of the definition of the DD range ρ qs (u) in (4) as follows ( Fig. 1) where matrix D n+1 ∈ R (n+1)×n forms between-entry differences [14]. Thus vector D T n+1 (·) ∈ R n contains the second-tolast entries of vector (·) ∈ R (n+1) , from which the first entry of (·) is subtracted. From (7) and the inequality ||u − u o || ≤ r o follows that (see Appendix) . Thus, the set (6) lies entirely inside an n-ball, representing a bounded set in R n .

III. BIAS-BOUNDED INTEGER ESTIMATION
As stated earlier, the integer ambiguity resolution procedures, that are often followed in practice, rely on the float ambiguity solutionẑ = y − Ax, thereby requiring the initial solutionx. In this section we present an integer estimator upon which y ∈ R m is directly mapped to the integer solutionž ∈ Z m on the basis of the three pieces of information in (5) only, i.e. Q yy ∈ R m×m , A ∈ R m×n and M ⊂ R n .

A. The Class of Admissible Integer Estimators
Before presenting the integer estimation sought, we first briefly review the existing theory of integer estimation that is based on the seminal contributions of Teunissen [11]- [13]. We discuss basic properties which an integer estimator is expected to possess. To uniquely specify an integer estimator of the unknown ambiguity vector z ∈ Z m , one needs to characterize its mapping y →ž from R m to Z m . Due to the discrete nature of Z m , such mapping is not one-to-one, but many-to-one. This means that for every integer-valued vector z o ∈ Z m , there will be a set of real-valued vectors w ∈ R m whose image is z o . Such set is referred to as the pull-in region of z o and can be expressed as [11] Thus the pull-in regions are subsets of R m . The integer estimatoř z, and therefore the mapping y →ž, can be specified by its corresponding pull-in regions as follows [12] By defining the pull-in regions (9) in an arbitrary fashion, one can therefore devise many integer estimators. However, such estimators and thus their pull-in regions must possess properties in order to be considered admissible in some sense. For instance, the union of the pull-in regions S z o (z o ∈ Z m ) must cover the whole m-dimensional space of real numbers, i.e. z o ∈Z m S z o = R m . Otherwise, there would be gaps in which not every y ∈ R m could be mapped to an integer-valued vector. Also, the pull-in regions should not have an overlap, where Int stands for the 'interior' of a set. Otherwise, y could be mapped to more than one integer-valued vectors. And finally, the pull-in regions S z o (z o ∈ Z m ) must be integer-translated copies of one another, i.e.
Accordingly, if the input vector y is shifted by an integer amount, the output integer solutionž must be shifted by the same integer amount. Any integer estimator, having the three properties stated above, belongs to the class of 'admissible integer estimators' [11].

Definition 2 (Admissible integer estimators):
The integermapping y →ž whose pull-in regions, given by (9), possess the following three properties is said to be an admissible integer estimator [11].
Two well-known examples of such admissible integer estimators, that are often used in GNSS, are discussed below.
Example 2 (Rounding and integer least-squares estimators): Suppose that the observation vector y ∈ R m is integer-mean. The corresponding model follows then from (5) by setting the bias Ax = E(y) − z to zero, i.e. x = 0. Thus E(y) = z. The goal is to estimate the unknown integer ambiguity vector z ∈ Z m . The simplest integer estimator of z follows by rounding the entries of y to their nearest integers. The rationale behind this is that the integer rounding delivers the estimatež R ∈ Z m as the closest integer vector to y with respect to the standard Euclidean norm || · ||, that isž withž R = y , and · being the entry-wise integer rounding operator. The integer rounding has been shown to be an admissible integer estimator with the following pull-in regions [12] with c i ∈ R m being a unit vector that has a '1' as its ith entry and zeros otherwise. The rounding pull-in regions (13) are m-dimensional cubes, centred at z 0 ∈ Z m , all having sides of length one.  These samples (red dots) will therefore be mapped to wrong integer vectors. As an estimate of the corresponding ambiguity success-rate, the ratio of the number of samples inside the blue square (green dots) to the total number of samples is computed as 98.7%. This estimator does not take the variance matrix Q yy , and therefore, the dispersion of the random vector y, into account. To simulate the samples in Fig. 2 which clearly differs from a scaled identity matrix, on which the weight of the standard Euclidean norm in (12) is based. To take Q yy into account, the standard norm || · || should be replaced by its weighted version || · || Q yy . This giveš The estimatorž ILS is known as the integer least-square (ILS) estimator [1]. In 1999, Teunissen [12] proved, for integer-mean normally distributed input vectors y, that the ILS estimator has the largest success rate among all admissible ambiguity estimators. The ILS pull-in regions are given as These z o -centred pull-in regions are formed by intersecting halfspaces in R m . Their two-dimensional examples are hexagons that are depicted in the bottom panel of Fig. 2. The shape and orientation of the ILS hexagons are governed by the variance matrix Q yy , thus following the dispersion of the input vector y. As shown in the figure, the z-centred blue hexagon captures the majority of the samples. Therefore, the ambiguity success rate considerably increases (from 98.7% to 99.9%) by switching from integer rounding to the ILS method. The price to pay for this spectacular performance is that the ILS solution (15) has to be obtained through an integer search. This is in contrast to the integer rounding with which the solution (12) directly follows from y . A computationally efficient way of finding the ILS solutionž ILS is to use the method of LAMBDA [1]. For a detailed discussion on the LAMBDA method, we refer to, e.g., [15] or [16].

B. A New Member
So far we have learned that the ILS estimation (15) is the optimal way of resolving the unknown integer ambiguity vector z of the bias-free model E(y) = z. It is optimal in the sense of achieving the highest probability of correct integer estimation, i.e. the largest ambiguity success-rate. The presence of the unaccounted bias Ax = E(y) − z can, however, significantly reduce the ILS success-rate and thus increase the chance of unsuccessful ambiguity resolution [17].
Instead of minimizing the squared-norm ||y − z|| 2 Q yy over z ∈ Z m , the idea now is to make use of the information content in the bias-bounded mixed-integer model (5) and minimize the squared-norm ||y −Ax −z|| 2 Q yy over both z ∈ Z m and x ∈ M. Application of such mixed-integer least-squares principle to solvable models is in fact not new and has been widely studied in several contributions, see e.g. [18]- [20]. To the author's knowledge however, its application to the model (5)-which without the constraints z ∈ Z m and x ∈ M is not solvable-is new. We treat the mixed-integer least-squares problem as a two-step minimization problem as follows where Thus in the first step, the squared-norm is to be minimized over x ∈ M, while the unknown vector z is kept fixed. This gives the objective function (18). In the second step, F y (z) is to be minimized over z ∈ Z m . We now define our new integer estimation on the basis of the objective function F y (z).

Definition 3 (Bias-bounded Estimation of AmbiguiTy):
With respect to the model (5) and the objective function (18), the integer estimation y →ž BEAT witȟ will be referred to as the method of BEAT. The corresponding pull-in regions are characterized as Compare the BEAT solution (19) with its ILS counterpart (15). For the trivial case M = {x o }, i.e. when M has only one member, BEAT coincides with the ILS estimation, be it that the input vector y is replaced by y For the general case M ⊂ R n however, these two integer estimators (and their performances) are different. Before discussing the BEAT performance, one first has to identify the conditions under which BEAT becomes an admissible integer estimator of z. Theorem 1 (Admissibility of BEAT): Consider the following minimizerš The method of BEAT (19) is an admissible integer estimator of the ambiguity vector z, in model (5), if and only if The necessary and sufficient condition (23) is needed to guarantee that every observation vector y ∈ R m gets mapped to a unique integer-valued vector. For the trivial case M = {x o } when BEAT reduces to the ILS method, all the minimizers (22) Thus the condition (23) reduces to z o − z l = 0 which always holds by definition. This makes sense as the ILS estimator is already an admissible integer estimator. Now recall from Section II that we imposed the constraint x ∈ M with the intention to rule out redundant solutions for u (thus for x). The natural question, with reference to the condition (23), is whether the constraint x ∈ M does what it is supposed to do. The answer follows as a direct consequence of Theorem 1.
Corollary 1 (Admissibility of BEAT): For the unconstrained case of model (5) in which the bounded set M ⊂ R n is replaced by M = R n , the necessary and sufficient condition (23) reduces to with the projector P ⊥ A = I m −A A + , and the least-squares inverse of A given by For the unconstrained case M = R n , the condition (24) states that BEAT is an admissible integer estimator, if no nonzero integer vector z o in Z m can be spanned by the columns of the design matrix A ∈ R m×n . Otherwise, there exists some vector In the next section we show that the stringent condition (24) almost never holds, as matrix A is often formed by rational numbers in practice, i.e. A ∈ Q m×n . The restriction upon which the unknown vector x is forced to be inside the bounded set M ⊂ R n therefore seems to be legitimate. But one still needs to quantify the extent to which the set M should be bounded. As finding the minimizers (22) for every z o ∈ Z m is not always a trivial task, the following corollary presents an alternative easy-to-verify sufficiency condition.
Corollary 2 (Admissibility of BEAT): Let x o ∈ M, and W ∈ R m×m be a given positive definite matrix. BEAT (19) is then an admissible integer estimator, if the bounded set M in (5) possesses the following property The upper-bound (25) enables one to check whether BEAT is admissible for a given bias-bounded mixed-integer model. It does, however, not enable one to conclude that BEAT is not admissible. It is thus possible for models to have BEAT admissible, while the sufficiency condition (25) does not hold. The workings of (25) can be seen as follows. Let the values, which the vector u can take on, imply the set ||x − x o || Q ≤ h. The task is to verify if the scalar h is smaller than the minimum of (1/2)||z|| W or not. With matrix W , the smallest weighted norm ||z|| W over the nonzero integer vectors z ∈ Z m can be readily computed by, e.g., the LAMBDA method [1], [15]. If 2 h is smaller than the computed minimum weighted norm, the method of BEAT is then declared to be admissible for that mixed-integer model. The following example provides an initial insight into the BEAT performance.

Example 3 (BEAT with different bounded sets):
Consider the two-dimensional observation vector y ∈ R 2 of Example 2 with variance matrix (14), but now with a noninteger mean given by Thus A = [0, 1] T (n = 1), and the bounded set M ⊂ R is given by the interval |x| ≤ h (thus x o = 0). The goal is to use BEAT for resolving the unknown integer vector z ∈ Z 2 . The choice W = I 2 also leads to Q = 1. As the smallest norm ||z|| over z ∈ Z 2 \{0} is equal to 1, the condition (25) is evaluated as |x| < 0.5. Accordingly, if h < 0.5, BEAT can be considered admissible for model (26).
To assess the BEAT success-rate, 100,000 normallydistributed samples of y are simulated. We choose the true values of the unknown parameters to be z = [0, 0] T and x = 0.14. The following three cases are considered: In comparison to the results of Example 2, the presence of the nonzero bias x lowers the ILS success-rate from 99.9% to 97.5%. Such reduction in the success-rate can be largely mitigated by taking the constraint |x| ≤ h (h > 0) into consideration. The success-rate increases to 99.3% for Case 2, i.e. when |x| ≤ 0.25. Upon choosing a less strict constraint |x| ≤ 0.35 (i.e. Case 3), the success-rate drops to 99.0% which is still larger than the bias-affected ILS success-rate.

IV. THE NONELLIPSOIDAL SEARCH SPACE OF BEAT
In the previous section, the method of BEAT as an admissible integer estimator was introduced and its performance for a simple two-dimensional example (Example 3) was briefly discussed. However, it has not yet been addressed how the BEAT estimators of z ∈ Z m and x ∈ M, given in model (5), can be obtained. Considering the two-step minimization problem (17), the BEAT estimators of z and x are expressed aš As previously shown in (6), the constraint x ∈ M can be parametrized in terms of a (nonlinear) function of the parameter vector u. This means, for a given vectorž ∈ Z m , that the second expression of (27) is just a constrained (nonlinear) least-squares problem. The minimizerx ∈ M can therefore be computed using standard trust-region techniques such as the 'Levenberg-Marquardt' method, see e.g. [21]. The computation of the first expression, i.e. the minimizerž, does however contain more involved steps.
As an indication of the complexity involved in thežcomputation, consider the definition of the objective function F y (z) in (18). In order to evaluate function F y (z) for every z ∈ Z m , one has to compute the corresponding minimizerx z given in (22). Therefore, one may be required to repeatedly solve the second minimization problem of (27) forx z over different integer vectors z ∈ Z m to eventually find the minimizerž. Note also the discrete domain of z ∈ Z m . As with the ILS solution (15), the minimizerž has to be searched inside a subset of Z m , known as the search space. The BEAT search space is characterized as where χ is a positive scalar which governs the size of the search space F. Given the search space (28), one may be inclined to take a naive approach and, for a given χ 2 , try to collect all integer vectors (members) that lie inside F, thereby evaluating the objective function F y (z) for all members. The member that returns the smallest value for F y (z) is then declared as the solutionž. In this regard, two questions arise: 1) how can one enumerate (identify) all the members of F? and 2) how can F be guaranteed not to have many members so that the evaluation of F y (z) is carried out in a timely manner? To address these questions, in the following we discuss the geometry of F so as to develop a computationally efficient integer search strategy for the computation ofž.

A. Encompassing Regions
To better appreciate the intricacies of the search space (28), we first express its defining objective function F y (z) as a sum of two squared-norms.
Lemma 1 (Decomposition of F y ): The BEAT objective function F y (z), given in (18), can be expressed as Clearly, the conditionx z ∈ M always holds for the unconstrained case M = R n . Thus when no restriction is placed on x, the BEAT search space F(χ 2 ) in (28) Such search space can therefore have infinitely many integer vectors (members). Now let us switch from the unconstrained case M = R n to a constrained case in which the bounded set M represents an n-ball, i.e. M = {x ∈ R n |||x − x o || ≤ h}. For this case, the BEAT objective function F y (z) can be shown to take the following form For the above special form of F y (z), the search space (28) follows as a union of two parts. The first part is the intersection of two degenerate m-dimensional hyper- The second part is formed by a union of the hyper-ellipsoids ||y −w −z|| 2 Q yy ≤ χ 2 , where w = Ax (x ∈ R n : ||x −x o || = h). This search space is visualized by the blue area in the right panel of Fig. 4. In contrast to the unconstrained case M = R n (left panel), the search space is now bounded, thereby containing finite members (integer vectors). For the twodimensional case shown in the figure, the intersection of the two stated degenerate hyper-ellipsoids is a 'parallelogram' of which the two sides, parallel to R (Q yy B), are extended by the two ellipses ||y−w −z|| 2 Q yy ≤ χ 2 , with w = A x (x = x o ± h). As shown in Fig. 4, both of the special forms 1) M = R n and 2) M = {x ∈ R n |||x − x o || ≤ h} lead to search spaces that are of a nonellipsoidal type, meaning that the search space cannot be described by a single hyper-ellipsoid. The complexity of such nonellipsoidal search space is driven by the bounded set M. This is in marked contrast to the search space of the ILS method. As indicated by the ILS objective function ||y x o −z|| 2 Q yy in (21), the ILS search space is ellipsoidal as it is given by the following hyper-ellipsoid (compare with 28) While all the members of an ellipsoidal search space can effectively be identified using the LAMBDA or MLAMBDA methods [1], [16], the identification of members inside a nonellipsoidal search space is often a challenging task [22]. This can potentially make the BEAT integer search a time-consuming process.
To have the LAMBDA ellipsoidal search strategy also applicable to a nonellipsoidal search space, Teunissen [22] proposed the usage of an ellipsoidal region that can encompass the search space. This would then allow one to identify all the members of the search space, plus those integer vectors lying outside the search space (but inside the region). By evaluating the objective function F y (z) for all the integer vectors inside the region, one identifies the minimizerž as the vector delivering the smallest value for F y (z). Clearly, the smaller the volume of such region, the fewer the number of to-be-searched integer vectors, thus the faster the integer search becomes. As illustrated in Fig. 4, the ILS search space E(χ 2 ) (red ellipse) always lies inside the BEAT search space F(χ 2 ) given in (28), i.e. E(χ 2 ) ⊂ F(χ 2 ), irrespective of the nonempty set M (see Appendix). Therefore, such ellipsoidal region cannot be used as it cannot encompass the entire BEAT search space. The following lemma presents two alternative encompassing regions. and in which the positive scalars χ 2 h and χ 2 h (B T z) are, respectively, given by with χ 2 (B T z) = χ 2 − ||t − B T z|| 2 Q tt and γ being the smallest eigenvalue of matrix Qxx.
If the BEAT search space F(χ 2 ) is nonempty, then the encompassing regions (33) and (34) are both guaranteed to contain the minimizerž in (27). The ILS and BEAT search spaces, given in (32) and (28), are related to the above two encompassing regions as (Appendix) together with the following property Thus for the special case h = 0, i.e. for the single-member set M = {x o }, the regions reduce to the ILS search space E(χ 2 ). The ellipsoidal region E(χ 2 h ) covers all the other regions. This region is visualized in the right panel of Fig. 4 (green ellipse). Given the ellipsoidal region (33), one can directly employ the LAMBDA method to enumerate all the integer vectors inside the region and pick the minimizer of F y (z). In the following example, we discuss how the members of such region are identified.
Example 4: (Enumeration inside an ellipsoidal search space) The goal is to find all the integer vectors inside the region ||y − z|| 2 Q yy ≤ χ 2 . Let the LDL decomposition of the variance matrix Q yy be given as Q yy = LDL T , with L a unit lower triangular matrix and D = diag(d 1 , . . . , d m ) a diagonal matrix. With y = [y i ] and L −1 = [l ij ] (i, j = 1, . . . , m), the following m sequential search intervals can be formulated [15] . . . , m). The scalars δ i and χ 2 i are, respectively, computed as with χ 2 1 = d 1 χ 2 . The procedure of finding the members of ||y − z|| 2 Q yy ≤ χ 2 can now be described as follows. In the first step, all the integers inside the interval I 1 are collected. For each of these integers, say z 1 , the second interval I 2 (z 1 ) is formed. For each of the integers inside I 2 (z 1 ), say z 2 , the third interval I 3 (z [2] ) is formed to give the corresponding integers z 3 . This sequential procedure is repeated for the remaining intervals 4, . . . , m). For some integer z i , the subsequent interval I i+1 (z [i] ) may happen to be empty. In this case, the integer z i is discarded and its next integer inside I i (z [i−1] ) is chosen to form I i+1 (z [i] ). In this way, all the integer vectors, say z = [z 1 , . . . , z m ] T , inside the hyper-ellipsoid ||y − z|| 2 Q yy ≤ χ 2 are identified.

B. The Volume of the Encompassing Regions
With the aid of the sequential search intervals (38) for the region (33), all the members of the BEAT search space (28) can be enumerated. This addresses the first of the two questions previously raised in the beginning of this section. However, it does not address the second question, namely, whether the region is guaranteed not to have many members. By avoiding having an abundance of unnecessary integer vectors in the search space, one can speed up the search for the minimizerž in (27). To address the second question, one needs to measure the volume of such region. As shown in [15], the volume of a search space can approximate the number of its members. To have a timely integer search, one therefore prefers the volume to be reasonably small. Here and in the following, the volume of a search space is referred to as the volume of the geometrical body that contains all real-valued vectors inside the search space. Accordingly, the volume of the ILS search space E(χ 2 ) is given by [15] Vol(E(χ 2 )) = χ m V m |Q yy | with | · | denoting the determinant of a matrix. The volume of the unit m-ball is given as V m = π m/2 /Γ(m/2 + 1). Thus the volume (40) monotonically increases with the scalar χ, the number of ambiguities m, and the noise-level of the measurements y which is characterized by |Q yy |.
With (40), one can infer the volume-ratio of the two ellipsoidal regions E(χ 2 h ) and E(χ 2 ), that is The above volume-ratio tells us how many times the encompassing region E(χ 2 h ) is larger than the ILS search space E(χ 2 ). This ratio is driven by the quantity (h/χ √ γ  'hybrid' structure (see Appendix) Thus the shape of F h (χ 2 ) depends on the dimension of the unknown vector x ∈ R n (i.e., n) and the scalar b = m − n. An expression for the volume of F h (χ 2 ) is presented below. Theorem 2 (Volume of F h ): The ratio of the volume of the encompassing region F h (χ 2 ), in (34), to the volume of the ILS search space E(χ 2 ), in (32), can be expressed as where n i , i = 0, . . . , n, are the 'binomial' coefficients. The volume-ratio (43) tells us how many times the encompassing region F h (χ 2 ) is larger than the ILS search space E(χ 2 ). For n = 0 and b = m, this ratio reduces to 1 as F h (χ 2 ) reduces to E(χ 2 ) (cf. 42). For n = m and b = 0 on the other hand, the ratio becomes equal to (41), since F h (χ 2 ) becomes E(χ 2 h ). Using (43), one can therefore measure to what extent F h (χ 2 ) is smaller than E(χ 2 h ) when b = 0. The top panel of Fig. 5 shows values of the volume-ratio (43) as a function of b and n. In the color map small values are indicated by green, while large values are indicated by red. As shown, the ratio increases with n, but remains almost bounded as b increases. The bottom panel of Fig. 5 shows this ratio as increasing functions of the 'scaled' bias-to-noise ratio (h/χ √ γ). The increase in the ratio is shown to be significantly smaller when b > n. For m = 50, the volumeratio (41), shown in red (n = 50, b = 0), can become as large as 10 15 . For the same dimension however, the volume-ratio (43) shown in green (n = 5, b = 45), reduces to about 10 3 , i.e. 10 12 times smaller. By switching from E(χ 2 h ) to F h (χ 2 ), one can thus considerably speed up the BEAT integer search when n is not too large and b > n. Fortunately, this is the case with applications like carrier phase-based positioning. As previously discussed in (3), the number of transmitters (satellites) are bounded (i.e. n is bounded), whereas each transmitter sends signals on multiple frequencies (i.e. m = fn). Thus b = (f − 1)n, meaning that the more the number of frequencies, the larger the scalar b becomes as compared to n. The encompassing region F h (χ 2 ), in (34), is therefore more suited to fast carrier phase-based positioning as it contains fewer unnecessary integer vectors.

V. A DUAL-ELLIPSOID SEARCH STRATEGY
In this section, a fast search strategy is presented for the BEAT solutionž, in (27), on the basis of the encompassing region F h (χ 2 ). As the region F h (χ 2 ) is 'nonellipsoidal,' it cannot be directly expressed into the m sequential search intervals (38). We therefore first present a canonical form for the bias-bounded mixed-integer model (5). It is then shown, under the canonical form, how the integer search inside F h (χ 2 ) can be conducted via two sequential LAMBDA ellipsoidal searches.

A. Bias-Bounded Model in Canonical Form
The idea is to transform the ambiguity vector z ∈ Z m , and therefore, the observation vector y ∈ R m into forms such that the defining terms of F h (χ 2 ), i.e. ||t − B T z|| 2 Q tt ≤ χ 2 and ||x z − x 0 || 2 Qxx ≤ χ 2 h (B T z), represent two non-degenerate hyper-ellipsoids in R b and R n , respectively. As it will be shown, this is done by using the change of variables [z T 1 , z T 2 ] T = z T z such that z T A = [0, H T ] T , where H ∈ R n×n is a nonsingular matrix. The transformation matrix z should preserve the integerness of the ambiguities upon switching from the 'original' vector z ∈ Z m to the 'transformed' vectors z 1 ∈ Z b and z 2 ∈ Z n , meaning that both z and its inverse z −1 must be formed by integer entries. Matrix z is thus required to be unimodular, i.e. z should be integer with determinant |z| = ±1. In the GNSS context, such matrix is referred to as the Z-transformation [23].
For matrices A ∈ R m×n whose column-spaces cannot be spanned by 'rational' vectors, Z-transformations with property z T A = [0, H T ] T , do not exist (see Appendix). In the following, the column-space of the design matrix A ∈ R m×n is therefore assumed to be spanned by a rational matrix N ∈ Q m×n , i.e. A = N x for some x ∈ R n×n . This is a plausible assumption, since in many, if not all, applications the coefficients of the underlying observation equations are formed by rational numbers. The existence of such Z-transformations is then guaranteed by virtue of the Hermite normal form for rational matrices [14], [24]. An immediate consequence of the assumption A = N x is that the condition (24) never holds when the set M is unbounded, i.e. when M = R n . This follows from the equality showing that the columns of the design matrix A span part of the columns of z −T . Thus where z o ∈ Z m can be a column of z −T . Therefore, the set M must be bounded in order to have an 'admissible' BEAT integer estimator. Another consequence of the assumption A = N x is the following.

Theorem 3 (Bias-bounded model in canonical form):
Let the Z-transformation z = [z 1 , z 2 ] ∈ Z m×m be given such that Then, the canonical form of the bias-bounded mixed-integer model (5) reads In the canonical form of the bias-bounded model (5), the observation vector y ∈ R m is replaced by the transformed observation vector [y T 1 , y T 2 ] T , with y 1 ∈ R b being 'integer-mean' (i.e. E(y 1 ) = z 1 ∈ Z b ), leaving the remaining n observations y 2 ∈ R n 'biased' asx = E(y 2 ) − z 2 . As the canonical model itself is a bias-bounded mixed-integer model, the results of the previous sections carry over to (44). Accordingly, one would only need to replace the design matrix A by [0, I n ] T , the bounded set M byM, and the variance matrix Q yy by D([y T 1 , y T 2 ] T ). Working with the canonical model (44) has the advantage that the BEAT solution of the transformed ambiguities [z T 1 , z T 2 ] T can be computed using a rather straightforward 'ellipsoidallike' search (see the next subsection). Applying the backtransformation to such solution, one can then recover the BEAT solution of the original ambiguities z.
Corollary 3 (BEAT solutions in canonical form): Given the canonical model (44), the BEAT solutions of [z T 1 , z T 2 ] T ∈ Z m andx ∈M are given as (ž 1 ,ž 2 ) = arg min in which the objective function F y (z 1 , z 2 ) reads The solutions of the original parameters (27) follow from the back-transformationž = z −T [ž T 1 ,ž T 2 ] T andx = H −1x . Compare the BEAT solutions (45) with their original counterparts (27). The objective function F y (z), in (29), is replaced by the objective function F y (z 1 , z 2 ) in (46). Before presenting the search strategy offered by the canonical model, we first give examples of the Z-transformation z ∈ Z m×m for the carrier phase observation equations (3).
Example 5: (Z-transformations for RIPS and GNSS) Teunissen [25] proposed a general algorithm that delivers Ztransformations for the design matrix A to fulfil the equality z T A = [0, H T ] T , see also [26]. For the 'carrier phase observation equations' (3) in which the design matrix takes the simple form of A = a ⊗ I n (a ∈ R f ), one can, however, compute such Z-transformations in a straightforward way. To see this, let matrix z ∈ Z m×m take the form of z = U ⊗ I n , where U ∈ Z f ×f is unimodular. Thus U must be integer with determinant |U | = ±1. The task is to find U such that where κ is a scalar. The last equality holds when U T a = [0, 1/κ] T . To find matrix U , let us now take a closer look at the f -vector a = [λ 1 /λ j ] (j = 1, . . . , f). As the wavelengthratios λ 1 /λ j are assumed to be rational numbers, they can be expressed as λ 1 /λ j = p j /p 1 , where the integer p j ∈ Z can, for instance, be the carrier phase frequency on the frequency-band j (j = 1, . . . , f). In case of RIPS, the f frequencies p j are often chosen to be equally spaced as [7] with the positive integer δ ∈ Z being the frequency channel separation. For this case, the unimodular matrix U can be chosen as The triple (g, c, d) = gcd(p 1 , δ) can be computed by the 'Extended Euclidean algorithm' [26], where scalar g ∈ Z is the greatest common divisor (gcd) of p 1 and δ, with the scalars c ∈ Z and d ∈ Z satisfying c p 1 + d δ = g. With reference to (48) and (50), the Z-transformation sought is therefore given by z = [z 1 , z 2 ], where z 1 = U 1 ⊗ I n and z 2 = U 2 ⊗ I n . It is now not difficult to verify, for A = a ⊗ I n , that z T 1 A = 0 and z T 2 A = (1/κ)I n , where κ = p 1 /g. With an application of the sufficiency condition (25), the scalar κ can be shown to play a decisive role in checking whether BEAT is admissible. Accordingly, in order for BEAT to be an admissible integer estimator for model (3), it is sufficient that the following condition holds (Appendix) are considered, it is sufficient to have u inside a sphere with a radius of r o ≈ 59.75m (as λ 1 ≈ 69cm if p 1 = 433MHz).
Let us now turn our attention to the GNSS carrier phase model. In case of GNSS, the f frequencies p j are not necessarily equally spaced. In that case, the unimodular matrix U , in z = U ⊗ I n , can be recursively computed through the following algorithmic steps [25] • j ← 2, (g, c, d) ← gcd(p j−1 , p j ), • • end while (52) The

B. A Sequential Dual-Ellipsoid Search
Given the canonical model (44), we are in a position to develop a search strategy to enumerate all the members of the encompassing region F h (χ 2 ) in (34). The stated region should encompass the BEAT search space corresponding to the canonical model (44), that is (cf. 46) Note that volume of the above search space is identical to that of (28), i.e. the BEAT search space corresponding to the orignal model (5). This is because the transformation [z T 1 , z T 2 ] T = z T z, with determinant |z| = ±1, is volume-preserving [1]. The following result shows how the region F h (χ 2 ), containing (54), can be expressed in terms of two non-degenerate hyper-ellipsoids.
Corollary 4 (Dual-ellipsoid search region): Upon the volume-preserving transformation [z T 1 , z T 2 ] T = z T z in (44), the region F h (χ 2 ), in (34), takes the following form in which two sequential hyper-ellipsoids are given by withx o = Hx o . The BEAT search space F(χ 2 ), given in (54), lies entirely inside the dual-ellipsoid region (55). Given the dual-ellipsoid region (55), a strategy to find the BEAT solutions (45) is as follows. To have a small search region, the search is initialized by a small nonnegative scalar χ 2 . Then, as with Example 4, the LAMBDA enumeration is conducted for the two sequential hyper-ellipsoids (56). The member with smallest objective value F y (z 1 , z 2 ) is thus declared as the BEAT solution. To speed up the search, the 'decorrelation' property of the LAMBDA method is employed in each of the hyperellipsoids E 1 (χ 2 ) and E 2 (χ 2 h (z 1 )). In the event that the choice of χ 2 leads to an empty search space, χ 2 is to be increased incrementally so as to have at least one member inside the search space. The following example is aimed to illustrate both the computational and quality performances of the BEAT method in carrier phase-based positioning.
Example 6: (BEAT performance for RIPS and GNSS) Two data-sets are considered: 1) simulated RIPS data (Fig. 6 ), and 2) real-world GNSS data (Fig. 7 ). We first discuss the RIPS results. As shown in the top-panel of Fig. 6, two radio receivers Rx-1 and Rx-2 collect beat phase measurements from 5 transmitters Tx-i (i = 1, . . . , 5). The task is to estimate the unknown state-vector of Rx-2 (u), given the known state-vectors of Rx-1 and the transmitters. As with Example 5, the first carrier frequency and channel separation are set to p 1 = 433MHz and δ = 5MHz, respectively (cf. 49). To analyze the execution time required for the BEAT integer search over different radii r o (||u − u o || ≤ r o ) when signals are transmitted on f = 11 frequencies, 10,000 normally-distributed beat-phase samples with the standard-deviation of 0.05 [cycles] are generated. The boxplots of the computed execution times are given in the middlepanel of Fig. 6, with medians ranging from 9 [ms] (r o = 10m) to 432 [ms] (r o = 30m). This demonstrates, for this example, that the search for the BEAT ambiguity solutionž ∈ Z m (m = fn = 44) can be conducted within fractions of a second. The bottom-panel of Fig. 6 shows BEAT ambiguity success-rates [%] for different radii r o and numbers of frequencies f . As shown, the success-rate increases as f increases, but decreases as r o increases. To show the role taken by the measurements' precision, the results when the beat-phase standard-deviation is    8. BEAT (black) and RTK (blue) performance compared: Commutative histograms of the wrongly-fixed solutions' radial error ||û − u||, corresponding to the short-baseline YARR-YAR3 (Fig. 7), when the number of visible satellites is limited to four. 0.03 [cycles] are presented inside the brackets. In this case, when f = 8 frequencies are considered, the ambiguity success-rate is not lower than 99.3%. Interestingly, the success-rate becomes 100% when the number of frequencies increases to f = 11. Now consider the GNSS results shown in Fig. 7. A daily Galileo quadruple-frequency dataset (cf. 53), with a 30-second sampling-rate, was collected from two IGS (International GNSS Service) stations YARR and YAR3 on 2 April 2021. Although the pseudo-range data could have also been utilized, we consider a single-epoch, phase-only scenario to illustrate the BEAT performance. However, one should bear in mind that the rover can make use of pseudo-range (code) data when they are available. Instead, the rover state-vector u is assumed to be a-priori known to lie inside the sphere ||u − u o || ≤ 8 metres. When the phase ambiguities are correctly fixed (green dots), the positioning errors are shown to be at the cm-level. It can also be observed that wrong ambiguity-fixing occurs (red dots) when only 4 to 5 Galileo satellites are visible, leading to positioning errors at the dm-level. To see how both the BEAT and standard GNSS realtime kinematic (RTK) methods respond to wrong ambiguityfixing when the number of visible satellites is limited to its minimum value of four, we present the commutative histograms of such solutions' radial error ||û − u|| in Fig. 8 . The figure illustrates, in case of a limited number of visible satellites, that BEAT delivers bounded positioning errors as compared to RTK if the phase ambiguities are wrongly fixed. This is because BEAT incorporates the constraint ||u − u o || ≤ 8m into the estimation process. When a modest number of satellites/transmitters are present, BEAT can thus be used for high-precision phase-only positioning. BEAT can therefore serve as a standard estimation method for applications like RIPS [7] or opportunistic navigation using Low-Earth-Orbiting (LEO) communication satellites [27] where the phase data are the only observations available.

VI. CONCLUSION
In this contribution, we developed and presented a new integer estimation method, the method of BEAT. BEAT employs the mixed-integer least-squares principle as optimality criterion, and extends the applicability of current integer estimation theory to bias-bounded mixed-integer models. This extension allows one to jointly estimate the integer vector z ∈ Z m and bounded real vector x ∈ M of the model (5), which would otherwise be unsolvable if prior knowledge of the 'bounded' set M ⊂ R n is not taken into consideration.
It was shown why and to what extent the set M should be bounded in order for BEAT to be an 'admissible' integer estimator. The admissibility conditions (11) are needed to ensure that every observation phase vector y ∈ R m gets mapped to a unique integer vectorž ∈ Z m . As with the integer least-squares estimation, the BEAT ambiguity estimationž ∈ Z m has to be searched, but inside a nonellipsoidal search-space. We therefore proposed a fast search strategy using a dual-ellipsoid region which encompasses the BEAT search-space.
Examples were presented to illustrate the BEAT positioning performance where, for the first time, the feasibility of GNSS single-epoch, phase-only positioning was shown. We believe that the proposed method opens a new path for the processing of phase-only measurements that are dedicated to ranging and navigation techniques, including those using the carrier phase signals of recent mega-constellation LEO satellites. This would, of course, demand supporting numerical studies that can illustrate the applicability of BEAT.

Supplementary Proofs
Proof of Equation (8) λ 1 ||D n+1 || ||u − u o || ≤ ||D n+1 || √ n + 1 (r o /λ 1 ) (57) The first inequality follows from the fact that the norm of an n-vector is never larger than √ n times its (in absolute value) largest entry. The second inequality follows from the triangleinequality while the third (last) inequality follows from the constraint ||u − u o || ≤ r o . The inequality (8) follows then from (57) if we show that ||D n+1 || = √ n + 1, or that the largest eigenvalue of D T n+1 D n+1 is n + 1. Let e n be an n-vector of ones. Then D n+1 = [−e n , I n ] T . Thus D T n+1 D n+1 = I n + e n e T n is an n × n symmetric matrix having n − 1 eigenvalues equal to 1, and one eigenvalue equal to n + 1.
Proof of Theorem 1: BEAT is admissible if its pull-in regions (20) satisfy the three conditions (11). The first condition of (11) already holds since for every y ∈ R m , the objective function F y (z) = min x∈M ||y − Ax − z|| 2 Q yy has bounded minimum value(s) for some integer vector z o . The third condition of (11) also holds as for every v ∈ S z o , BEAT , there exists w = v − z o such that ||t − B T z|| 2 Q tt (when n = 0) and ||y x o − z|| 2 Q yy = ||x z − x o || 2 Qxx (when n = m). Proof of Theorem 2: Let w ∈ R m represent all real vectors inside the region F h (χ 2 ). The corresponding geometrical body is given by The volume of the region F h (χ 2 ) can then be obtained by evaluating the following integral Vol(F h (χ 2 )) = The first equality is due to the definition of F l , while the second equality follows as V nχ n (v 1 ) is the volume of the n-ball ||v 2 || ≤χ(v 1 ). The third equality is due to the definition of χ(v 1 ), while the fourth equality follows from an application of the n-dimensional spherical coordinates rule for integrals [28], with radius α = ||v 1 || and angles θ i (i = 1, . . . , b − 1). The fifth (last) equality follows, as the product of the trigonometric integrals in (71) is identical to bV b . Using the change of variable α = sin(η), along with an application of the binomial theorem, the integral In the equality expression, we define the projector P D n+1 = D n+1 (D T n+1 D n+1 ) −1 D T n+1 . The first inequality follows from the fact that the eigenvalues of a projector are zeros or ones. The second inequality follows from the fact that the norm of an n-vector is never larger than √ n times its (in absolute value) largest entry, and the triangle-inequality (58).