Using SDP Relaxation for Optimization of the Satellites Set Chosen for Positioning

Abstract

A method for determining attitude of an object having multiple GNSS antennas, the method including receiving GNSS signals from at least five satellites, wherein at least 2 of the five belong to a different satellite constellation than the other satellites; processing each of the GNSS signals to generate pseudorange code and carrier phase measurements; resolving carrier phase ambiguities for all the received GNSS signals; generating unbiased carrier phase measurements based on the resolving; determining the attitude, including heading, pitch, and roll angles ,,, respectively, by solving a quadratically constrained quadratic minimization problem through finding a minimum of a linear function subject to a linear matrix inequality constraint; and outputting the attitude.

Claims

1. A method for determining attitude of an object having multiple GNSS antennas, the method comprising: receiving GNSS signals from at least five satellites, wherein at least two of the five belong to a different satellite constellation than the other satellites; processing each of the GNSS signals to generate pseudorange code and carrier phase measurements; resolving carrier phase ambiguities for all the received GNSS signals; generating unbiased carrier phase measurements based on the resolving; determining the attitude, including heading, pitch, and roll angles , , , respectively, by solving a quadratically constrained quadratic minimization problem through finding a minimum of a linear function subject to a linear matrix inequality constraint; and outputting the attitude.

2. The method of claim 1, wherein the finding of the minimum uses Semi-Definite Programming.

3. The method of claim 1, wherein the quadratically constrained quadratic minimization problem solves for the matrix Q in $\underset{Q^T.Math.Q=I}{\min }.Math.\left[\mathrm{tr}.Math.X^{0.Math.T}.Math.Q^T.Math.{\stackrel{\_}{J}}^T.Math.W.Math.\stackrel{\_}{J}.Math.{\mathrm{QX}}^0-2.Math.\mathrm{tr}.Math.{\stackrel{\_}{B}}^T.Math.W.Math.\stackrel{\_}{J}.Math.{\mathrm{QX}}^0\right]$ where Q is an orthogonal matrix $Q=\left(\begin{array}{ccc}\cos \left(\right).Math.\cos \left(\right)& \cos \left(\right).Math.\sin \left(\right).Math.\sin \left(\right)-\cos \left(\right).Math.\sin \left(\right)& \sin \left(\right).Math.\sin \left(\right)+\cos \left(\right).Math.\cos \left(\right).Math.\sin \left(\right)\\ \cos \left(\right).Math.\sin \left(\right)& \cos \left(\right).Math.\cos \left(\right)+\sin \left(\right).Math.\sin \left(\right).Math.\sin \left(\right)& \cos \left(\right).Math.\sin \left(\right).Math.\sin \left(\right)-\cos \left(\right).Math.\sin \left(\right)\\ -\sin \left(\right)& \cos \left(\right).Math.\sin \left(\right)& \cos \left(\right).Math.\cos \left(\right)\end{array}\right),$ W is a weighting matrix, B is a matrix of unbiased carrier phase across-antenna differenced measurements collected for all baselines and all signals, J is a matrix of directional cosines corresponding to all signals, and X.sup.0 is a matrix of all the baselines expressed in a body frame of the object.

4. The method of claim 3, wherein the matrix W is a diagonal weighting matrix representing an inverse of variance for each satellite.

5. The method of claim 3, wherein the attitude determination uses baselines defined by vector distances between the antennas, and wherein the baselines are selected using a spanning tree to minimize a condition number of the matrix X.sup.0X.sup.0T.

6. The method of claim 3, wherein the matrix Q must satisfy LMI (linear matrix inequality) of $Z0,.Math.\left(\begin{array}{cc}I& Q\\ Q^T& I\end{array}\right)0,$ where $Z=\left(\begin{array}{cc}Y& {\mathrm{QX}}^0\\ X^{0.Math.T}.Math.Q^T& I\end{array}\right)$ is a matrix linearly dependent on the matrix Q Y is auxiliary symmetric 33 matrix, and I is a 33 identity matrix.

7. A system for determining attitude of an object, the system comprising: a plurality of antennas receiving GNSS signals from at least five satellites, wherein at least two of the five belong to a different satellite constellation than the other satellites; a signal processor processing each of the GNSS signals to generate pseudorange code and carrier phase measurements; wherein the signal processor resolves carrier phase ambiguities for all the pseudorange code and carrier phase measurements and generates unbiased carrier phase measurements based on the resolving; an attitude determination block determining and outputting the attitude, including heading, pitch, and roll angles ,,, respectively, by applying Semi-Definite Programming to the unbiased carrier phase measurements.

8. The system of claim 7, wherein the linear function is G,Z, where $G=\left(\begin{array}{cc}{\stackrel{\_}{J}}^T.Math.W.Math.\stackrel{\_}{J}& -{\stackrel{\_}{J}}^T.Math.W.Math.\stackrel{\_}{B}\\ -{\stackrel{\_}{B}}^T.Math.W.Math.\stackrel{\_}{J}& 0\end{array}\right),$ Q is an orthogonal matrix $Q=\left(\begin{array}{ccc}\cos \left(\right).Math.\cos \left(\right)& \cos \left(\right).Math.\sin \left(\right).Math.\sin \left(\right)-\cos \left(\right).Math.\sin \left(\right)& \sin \left(\right).Math.\sin \left(\right)+\cos \left(\right).Math.\cos \left(\right).Math.\sin \left(\right)\\ \cos \left(\right).Math.\sin \left(\right)& \cos \left(\right).Math.\cos \left(\right)+\sin \left(\right).Math.\sin \left(\right).Math.\sin \left(\right)& \cos \left(\right).Math.\sin \left(\right).Math.\sin \left(\right)-\cos \left(\right).Math.\sin \left(\right)\\ -\sin \left(\right)& \cos \left(\right).Math.\sin \left(\right)& \cos \left(\right).Math.\cos \left(\right)\end{array}\right),$ W is a weighting matrix, B is a matrix of unbiased carrier phase across-antenna differenced measurements collected for all baselines and all signals, J is a matrix of directional cosines corresponding to all signals, and X.sup.0 is a matrix of all the baselines expressed in a body frame of the object; where $Z=\left(\begin{array}{cc}Y& {\mathrm{QX}}^0\\ X^{0.Math.T}.Math.Q^T& I\end{array}\right)$ is a matrix linearly dependent on the matrix Q, Y is auxiliary symmetric 33 matrix, and I is a 33 identity matrix.

9. The system of claim 7, wherein the matrix W is a diagonal weight matrix representing an inverse of variance for each satellite.

10. The system of claim 7, wherein the attitude determination uses baselines defined by vector distances between the antennas, and wherein the baselines are selected using a spanning tree to minimize a condition number of the matrix X.sup.0X.sup.0T.

11. The system of claim 7, wherein the matrix Q must satisfy LMI (linear matrix inequality) of $Z0,.Math.\left(\begin{array}{cc}I& Q\\ Q^T& I\end{array}\right)0,$ where $Z=\left(\begin{array}{cc}Y& {\mathrm{QX}}^0\\ X^{0.Math.T}.Math.Q^T& I\end{array}\right)$ is a matrix linearly dependent on the matrix Q, Y is auxiliary symmetric 33 matrix, and I is a 33 identity matrix.

Description

BRIEF DESCRIPTION OF THE ATTACHED FIGURES

[0089] The accompanying drawings, which are included to provide a further understanding of the invention and are incorporated in and constitute a part of this specification, illustrate embodiments of the invention and together with the description serve to explain the principles of the invention.

[0090] In the drawings:

[0091] FIG. 1 illustrates different reference coordinate frames.

[0092] FIG. 2 illustrates full antenna location graph.

[0093] FIG. 3 illustrates examples of different spanning trees.

[0094] FIG. 4 shows a structure of an exemplary multi-antenna navigation receiver.

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

[0095] Reference will now be made in detail to the embodiments of the present invention.

[0096] With the spanning tree determined above and consisting of (N1) edges (unknown vectors) d x.sub.kl, let us combine all vectors b*.sub.,kl(t) to the M(N1) matrix B (the notation (t) is omitted for brevity). All (N1) vectors can be similarly combined into the 3(N1) matrix (unknown) X, all values .sub.kl(t) can be combined into the 1(N1) row . Let also define the matrix J=.sup.1A and vector E=.sup.1e. Then collection of (N1) equations (12) can be rewritten as the matrix equation where all dimensions match according to definitions determined in this paragraph:

B=JX+E(21)

[0097] Substitution of the relationship X=QX.sup.0 into the equation (21) and applying the least squares method leads to the following formulation:

[00013]$\begin{array}{cc}{.Math.{\mathrm{JQX}}_0+E.Math..Math.-B.Math.}_F^2.fwdarw.\underset{Q,}{\min }& \left(22\right)\end{array}$

[0098] over all matrices subject to Q.sup.TQ=I and real valued , where I stands for the identity matrix of the appropriate dimension.

[0099] The minimization problem (22) can now be rewritten as

[00014]$\begin{array}{cc}\underset{Q^T.Math.Q=I}{\min }.Math.\left(\underset{}{\min }.Math.{.Math.{\mathrm{JQX}}^0+E.Math..Math.-B.Math.}_F^2\right)& \left(23\right)\end{array}$

[0100] Taking the internal minimum in (23) and using the trace notation for the Frobenius norm, (23) can be rewritten as

[00015]$\begin{array}{cc}\underset{Q^T.Math.Q=I}{\min }.Math.{\mathrm{tr}\left(\stackrel{\_}{J}.Math.{\mathrm{QX}}^0-\stackrel{\_}{B}\right)}^T.Math.\left(\stackrel{\_}{J}.Math.{\mathrm{QX}}^0-\stackrel{\_}{B}\right)& \left(24\right)\end{array}$

where J and B are matrices defined as

[00016]$J=\left(I-\frac{1}{E^T.Math.E}.Math.{\mathrm{EE}}^T\right).Math.J,\stackrel{\_}{B}=\left(I-\frac{1}{E^T.Math.E}.Math.{\mathrm{EE}}^T\right).Math.B.$

[0101] Usually optimal processing assumes weighting of measurements. The maximum likelihood method uses an inverse variance-covariance matrix for weighting: W=C.sup.1. Knowledge of the variance-covariance matrix is assumed because we know the physical nature of errors. Therefore, the problem of equation (24) takes the form

[00017]$\begin{array}{cc}\underset{Q^T.Math.Q=I}{\min }.Math.{\mathrm{tr}\left(\stackrel{\_}{J}.Math.{\mathrm{QX}}^0-\stackrel{\_}{B}\right)}^T.Math.W\left(\stackrel{\_}{J}.Math.{\mathrm{QX}}^0-\stackrel{\_}{B}\right)& \left(25\right)\end{array}$

with appropriate weighting matrix dimension.

[0102] Using properties of the matrix trace, (24) can be rewritten as

[00018]$\begin{array}{cc}\underset{Q^T.Math.Q=I}{\min }.Math.\left[\mathrm{tr}.Math.X^{0.Math.T}.Math.Q^T.Math.{\stackrel{\_}{J}}^T.Math.W.Math.\stackrel{\_}{J}.Math.{\mathrm{QX}}^0-2.Math.\mathrm{tr}.Math.{\stackrel{\_}{B}}^T.Math.W.Math.\stackrel{\_}{J}.Math.{\mathrm{QX}}^0\right]& \left(26\right)\end{array}$

[0103] The problem with solving expression (26) subject to minimization, is that it is a nonlinear programming problem, since it assumes a minimization of the scalar nonlinear function (the expression in the square brackets) of 9 variables (entries of the matrix Q) over 6 nonlinear quadratic constraints (diagonal and upper diagonal entries of the matrix Q.sup.TQ). Generally, these problems are hard to solve reliably and in a predictable time. In many practical applications (such as a small aircraft performing surveying), expression (26) has to be solved many times per second (sometimes several hundred times per second). This is computationally impractical.

[0104] Summing up, in the general case, the expression (26) is hard to solve because of quadratic dependence of the cost function on the matrix Q and the quadratic constraints.

[0105] The problem (26) is formulated after the optimal spanning tree has been selected among all possible spanning trees in the full antenna graph. The natural way to optimize selection is to look for such a tree that provides minimum value of the condition number of the matrix X.sup.0X.sup.0T. The smaller the condition number, the better the entire problem (26) is numerically conditioned. Large values mean ill-conditioning. Small values mean well-conditioning.

[0106] The problem (26) can be simplified if we chose another, not optimal weighting matrix in such a way that J.sup.TWJ=I (for example W=J(JJ.sup.T).sup.2J.sup.T=J.sup.+TJ.sup.+, J.sup.+ is pseudoinverse). Then the problem (26) takes the form

[00019]$\begin{array}{cc}\underset{Q^T.Math.Q=I}{\max }.Math.\mathrm{tr}.Math.{\stackrel{\_}{B}}^T.Math.W.Math.\stackrel{\_}{J}.Math.{\mathrm{QX}}^0=\underset{Q^T=Q=I}{\max }.Math.\mathrm{tr}.Math.{\tilde{X}}^T.Math.{\mathrm{QX}}^0,.Math.{\tilde{X}}^T={\stackrel{\_}{B}}^T.Math.W.Math.\stackrel{\_}{J}& \left(27\right)\end{array}$

which is similar to (20) and can be efficiently solved via SVD as shown above. Therefore, the methods proposed in [1] (and commonly used conventionally) can be derived from (26) with an artificial and non-optimal choice of the weighting matrix. Phrased another wayin the conventional approach of [1], a computationally efficient solution of (27) for matrix Q uses a non-optimal matrix W. The desired solution is therefore a solution of (26), rather than a solution of (27). To do this, the SDP relaxation approach is used, to make the solution of (26) tractable.

[0107] Non-convexity of the optimization problem (26) makes it very inefficient for computationally efficient solution. To make the problem numerically tractable we use semi-definite programming (an SDP approach to relax non-convex constraints) and the cost function replacing them by their convex analogs.

[0108] The convex problem allows for a numerical solution with polynomial computational complexity (see [4]). The SDP relaxation applied to the problem of solving equation (26) for attitude determination using GNSS signals and the algorithm for solution of the attitude determination problem is therefore a novel step of the present invention.

[0109] After the matrix Q is determined, it allows for calculation of orientation of the body frame with respect to ECEF. To translate it to the orientation of the body frame relative to the local horizon it must be multiplied by appropriate rotation matrix, rotating the reference frame from ECEF to the local horizon, see for example [3]. After that the three attitude/orientation angles ,, (which are the quantities that are ultimately of interest in this exercise) can be determined according to (1).

[0110] The expression subject to minimization in (26) can be equivalently transformed using concept of inner product of matrices. Let R.sub.1 and R.sub.2 be square matrices of the same size. Define the inner product of them as

R.sub.1,R.sub.2=trR.sub.1R.sub.2 (28)

so that R,R=R.sub.F.sup.2. Let us denote

[00020]$\begin{array}{cc}G=\left(\begin{array}{cc}{\stackrel{\_}{J}}^T.Math.W.Math.\stackrel{\_}{J}& -{\stackrel{\_}{J}}^T.Math.W.Math.\stackrel{\_}{B}\\ -{\stackrel{\_}{B}}^T.Math.W.Math.\stackrel{\_}{J}& 0\end{array}\right)& \left(29\right)\end{array}$

[0111] and keeping in mind the relation X=QX.sup.0 to denote

[00021]$\begin{array}{cc}Z=\left(\begin{array}{cc}Y& X\\ X^T& I\end{array}\right)& \left(30\right)\end{array}$

where Y=XX.sup.T.

[0112] Then the problem (26) is equivalent to the following problem

G,Z.fwdarw.min (31)

[0113] subject to constraints

[00022]$\begin{array}{cc}Y-{\mathrm{QX}}^0.Math.X^{0.Math.T}.Math.Q^T=0,& \left(32\right)\\ Z=\left(\begin{array}{cc}Y& {\mathrm{QX}}^0\\ X^{0.Math.T}.Math.Q^T& I\end{array}\right),& \left(33\right)\\ Q^T.Math.Q-I=0,& \left(34\right)\end{array}$

[0114] Conditions (32) and (34) are non-convex. Substituting for them less constraining semi-definite inequality constraints

YQX.sup.0X.sup.0TQ.sup.T0 (35)

and

Q.sup.TQI0 (36)

means semidefinite relaxation. Here the symbol means positive semi-definiteness of matrices. Condition (35) is equivalent (due to the so called Schur Lemma, see [4]) to the condition

Z0. (37)

[0115] Condition (36) is equivalent to

[00023]$\begin{array}{cc}\left(\begin{array}{cc}I& Q\\ Q^T& I\end{array}\right)0& \left(38\right)\end{array}$

[0116] Note that the cost function subject to minimization in (31) is linearly dependent on 9 entries of the matrix Q due to the equation (33). Left sides of matrix inequalities (37) and (38) are also linearly dependent on entries of Q. These inequalities are known as linear matrix inequalities (LMI) and establish convex constraints on Q. Thus, the attitude determination problem is reduced to the convex optimization problem for which there are very efficient numerical algorithms, see [4].

[0117] Possible algorithms include: [0118] Ellipsoid method, and [0119] Path-following method for -log det(.) barrier functions.

[0120] These and other similar algorithms are known in the art, and many of them are readily available in Linux.

[0121] Numerical experiments performed by the inventor shows that in practice relaxed and initial problems in many examples give the same result. Therefore, in practice, relaxation does not necessarily lead to a need to improve the result applying additional orthogonalization of the matrix

[0122] If the number of antennas is 3 then N1=2 and number of columns in the matrix X.sub.0 is 2. In this case (and only in this case) the number of columns must be increased by adding third column as an outer product (vector product) of first two columns

d x.sub.3.sup.0=d x.sub.1.sup.0d x.sub.2.sup.0 (39)

[0123] The choice of the spanning tree in the full graphs (FIG. 2) should be performed in such a way that the condition number of the matrix X.sup.0TX.sup.0 is minimized. This problem is solved only once, not in real time and therefore its computational impact is not a limiting factor.

[0124] The structure of the multi-antenna navigation receiver is shown in FIG. 4. Multiple antennas are connected to the multi-antenna receiver consisting among others of multi-channel radio frequency (RF) part. RF signals are convolved with reference signals generated by the numerically controlled oscillator and digitized by A/D converter. Result of convolution is used for multiple signal tracking. Multiple signals corresponding to plurality of satellites, multiple frequency bands, and multiple antennas are used for attitude calculations as described in the patent. The processor runs signal tracking, positioning, and attitude determination algorithms.

[0125] Having thus described the different embodiments of a system and method, it should be apparent to those skilled in the art that certain advantages of the described method and apparatus have been achieved. It should also be appreciated that various modifications, adaptations, and alternative embodiments thereof may be made within the scope and spirit of the present invention. The invention is further defined by the following claims.

REFERENCES (ALL OF WHICH ARE INCORPORATED HEREIN BY REFERENCE IN THEIR ENTIRETY)

[0126] 1. C. E. Cohen, Attitude Determination, Global Positioning System: Theory and Applications, Volume 2, edited by B. W. Parkinson and J. J. Spilker, Vol. 164, Progress in Astronautics and Aeronautics, AIAA, Reston, Va., 1994. [0127] 2. C. Berge, The Theory of Graphs, Courier Corp., 1962. [0128] 3. A. Leick, L. Rapoport, D. Tatarnikov, GPS Satellite Surveying, Wiley & Sons, 2015 [0129] 4. Ben-Tal, A., Nemirovski, A., Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications, MOS-SIAM Series on Optimization, 2001.