A dual-system weighting precision factor direct calculation method

Abstract

The application discloses a double-system weighted precision factor direct calculation method, which is applied to the technical field of navigation and aims at the following problems of the prior art: a large amount of calculation is needed for matrix inversion for multiple times in the calculation of WDOP value; and numerical stability problems are prone to occurring in the matrix inversion process when a geometric observation matrix is in a state of illness. * The expression of the weighted precision factor (WDOP) is transformed, det(N) and tr(N * ) are calculated respectively, and thus the WDOP value is solved, so that the matrix inversion operation is not needed for multiple times, and the calculation efficiency is improved and the calculation stability is effectively ensured by the method.

Description

Technical Field

[0001] This invention belongs to the field of navigation technology, and specifically relates to a direct calculation technique for the weighted accuracy factor (WDOP) of a multi-mode satellite navigation system. Background Technology

[0002] Multi-mode satellite navigation systems, which combine the BeiDou system with other satellite navigation systems (such as the Global Positioning System, GPS), will gradually become the mainstream development direction and inevitable trend in the fields of navigation, positioning, velocity measurement, and timing. In multi-mode satellite navigation systems, the number of satellites that terminal equipment can receive will increase significantly, thereby improving the system's positioning accuracy, availability, continuity, and integrity.

[0003] In positioning calculations, the Geometric Dilution of Precision (GDOP) is crucial for satellite selection by terminal devices, positioning accuracy assessment, and system ill-conditioning diagnosis. However, GDOP typically assumes that all currently visible satellites have the same ranging accuracy, failing to consider the impact of ranging accuracy from different satellites (especially those from different systems) on the positioning calculation process. In reality, for multi-system terminal device positioning calculations, due to factors such as satellite ephemeris errors, satellite elevation angles, and carrier-to-noise ratios, the ranging accuracy of different satellites from various satellite navigation systems (such as BeiDou and GPS satellites) varies significantly. Therefore, in multi-system positioning calculations, a weighted matrix is ​​usually introduced to measure the ranging accuracy of different satellites. Accordingly, the Geometric Dilution of Precision (WDOP) under the influence of the weighted matrix is ​​used for satellite selection by multi-system terminal devices, positioning accuracy assessment, and system ill-conditioning diagnosis.

[0004] In the process of positioning and solving multi-system terminal equipment, calculating the WDOP value using traditional methods requires multiple matrix inversions, which significantly increases the computational load. Furthermore, when the geometric observation matrix is ​​ill-conditioned, the matrix inversion process is prone to numerical stability problems. Therefore, avoiding the matrix inversion process and realizing direct calculation of multi-system WDOP is of great significance. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention proposes a direct calculation method for the multi-system weighted precision factor (WDOP). This method eliminates the need for multiple matrix inversion operations and directly provides the mathematical expression for WDOP.

[0006] The technical solution adopted in this invention is: a method for direct calculation of multi-system weighted precision factor (WDOP), comprising:

[0007] S1. When the BD / GPS terminal device observes 3 BD satellites and 2 GPS satellites, the geometric observation matrix H is obtained;

[0008] S2. Based on the geometric observation matrix H and the weighting matrix W, the expression for the weighted precision factor (WDOP) is obtained:

[0009]

[0010] Where tr(·) denotes the trace of a matrix, (·) T This indicates the matrix transpose; the weighted matrix W is a diagonal matrix, and its diagonal elements represent the weights of the i-th satellite;

[0011] S3, Let N = H T WH, transform the weighted precision factor expression in step S2 into:

[0012]

[0013] Where det(·) denotes the calculation of the determinant of a matrix, (·) * This indicates the calculation of the adjoint matrix;

[0014] S4. By calculating det(N) and tr(N) respectively * ), to obtain the WDOP value.

[0015] The process of calculating det(N) in step S4 is as follows:

[0016] A1. When 3 BD satellites and 2 GPS satellites are observed, since both the geometric observation matrix H and the weighting matrix W are square matrices, the determinant det(N) can be simplified to:

[0017] det(N) = det(H) T WH)=[det(H)] 2 det(W)

[0018] A2. Since the weighting matrix W is a diagonal matrix, its determinant can be simplified to:

[0019]

[0020] A3. After Schur decomposition and corresponding algebraic operations, the determinant det(H) has the following properties:

[0021]

[0022] A4. Based on the expressions in steps A2 and A3, transform det(N) into:

[0023]

[0024] Step S4 describes the calculation of tr(N) * The process is as follows:

[0025] B1, Order Represents the adjoint matrix N * The i-th diagonal element then has

[0026]

[0027] B2. Since N = H T WH=(n ij ) 5×5 Given a 5th-order real symmetric matrix, its adjoint matrix N is derived. * The diagonal elements are as follows:

[0028]

[0029]

[0030]

[0031]

[0032]

[0033] B3. According to step B2, tr(N) * Transformed into:

[0034] tr(N * ) = n 11 A+n 12 B+n 13 C+n 14 D+n 15 E+n 22 F+n 23 G+n 24 H+n 25 I

[0035] Where A, B, C, D, E, F, G, H, and I represent correlation coefficients.

[0036] The formula for calculating the correlation coefficient in step B3 is as follows:

[0037] A = n 33 (n 44 n 55 -n 45 n 54 )+n 34 (n 45 n 53 -n 43 n 55 )+n35 (n 43 n 54 -n 44 n 53 )+n 22 (n 44 n 55 -n 45 n 54 )+n 24 (n 45 n 52 -n 42 n 55 )+n 25 (n 42 n 54 -n 44 n 52 )+n 22 (n 33 n 55 -n 35 n 53 )+n 23 (n 35 n 52 -n 32 n 55 )+n 25 (n 32 n 53 -n 33 n 52 )+n 22 (n 33 n 44 -n 34 n 43 )+n 23 (n 34 n 42 -n 32 n 44 )+n 24 (n 32 n 43 -n 33 n 42 )

[0038] B=n 21 (n 45 n 54 -n 44 n 55 )+n 24 (n 41 n 55 -n 45 n 51 )+n 25 (n 44 n 51 -n 41 n 54 )+n 21(n 35 n 53 -n 33 n 55 )+n 23 (n 31 n 55 -n 35 n 51 )+n 25 (n 33 n 51 -n 31 n 53 )+n 21 (n 34 n 43 -n 33 n 44 )+n 23 (n 31 n 44 -n 34 n 41 )+n 24 (n 33 n 41 -n 31 n 43 )

[0039] C=n 31 (n 45 n 54 -n 44 n 55 )+n 34 (n 41 n 55 -n 45 n 51 )+n 35 (n 44 n 51 -n 41 n 54 )+n 21 (n 32 n 55 -n 35 n 52 )+n 22 (n 35 n 51 -n 31 n 55 )+n 25 (n 31 n 52 -n 32 n 51 )+n 21 (n 32 n 44 -n 34 n 42 )+n 22 (n34 n 41 -n 31 n 44 )+n 24 (n 31 n 42 -n 32 n 41 )

[0040] D=n 31 (n 43 n 55 -n 45 n 53 )+n 33 (n 45 n 51 -n 41 n 55 )+n 35 (n 41 n 53 -n 43 n 51 )+n 21 (n 42 n 55 -n 45 n 52 )+n 22 (n 45 n 51 -n 41 n 55 )+n 25 (n 41 n 52 -n 42 n 51 )+n 21 (n 33 n 42 -n 32 n 43 )+n 22 (n 31 n 43 -n 33 n 41 )+n 23 (n 32 n 41 -n 31 n 42 )

[0041] E=n 31 (n 44 n 53 -n 43 n 54 )+n 33 (n 41 n 54 -n 44 n 51 )+n34 (n 43 n 51 -n 41 n 53 )+n 21 (n 44 n 52 -n 42 n 54 )+n 22 (n 41 n 54 -n 44 n 51 )+n 24 (n 42 n 51 -n 41 n 52 )+n 21 (n 33 n 52 -n 32 n 53 )+n 22 (n 31 n 53 -n 33 n 51 )+n 23 (n 32 n 51 -n 31 n 52 )

[0042] F=n 33 (n 44 n 55 -n 45 n 54 )+n 34 (n 45 n 53 -n 43 n 55 )+n 35 (n 43 n 54 -n 44 n 53 )

[0043] G=n 32 (n 45 n 54 -n 44 n 55 )+n 34 (n 42 n 55 -n 45 n 52 )+n 35 (n 44 n 52 -n 42 n54 )

[0044] H = n 32 (n 42 n 55 -n 45 n 53 )+n 33 (n 45 n 52 -n 42 n 55 )+n 35 (n 42 n 53 -n 43 n 52 )

[0045] I = n 32 (n 43 n 54 +n 44 n 53 )+n 33 (n 42 n 54 -n 44 n 52 )+n 34 (n 43 n 52 -n 42 n 53 ).

[0046] The beneficial effects of this invention are as follows: Traditional methods for calculating the WDOP value require multiple matrix inversions, significantly increasing the computational load; furthermore, when the geometric observation matrix is ​​ill-conditioned, the matrix inversion process is prone to numerical stability issues. This invention transforms the expression for the weighted precision factor (WDOP) by calculating det(N) and tr(N) separately. * This method solves the WDOP value without requiring multiple matrix inversion operations. The method of this invention not only improves computational efficiency but also effectively ensures computational stability. Attached Figure Description

[0047] Figure 1 This is a flowchart of the present invention. Detailed Implementation

[0048] To facilitate understanding of the technical content of this invention by those skilled in the art, the following description, in conjunction with the accompanying drawings, further illustrates the invention.

[0049] Taking a BD / GPS terminal device as an example, the positioning calculation process requires solving for five unknown parameters (including 3D position information and two time parameters), thus requiring at least five satellites for positioning calculation. This invention primarily considers the direct calculation method for WDOP under the condition of five satellites. In BD / GPS terminal devices, when five satellites are observed, this invention mainly discusses the case of three BD (or GPS) satellites and two GPS (or BD) satellites.

[0050] like Figure 1 As shown, the implementation process of this invention is as follows:

[0051] S1. Without loss of generality, when the BD / GPS terminal device observes 3 BD satellites and 2 GPS satellites, the geometric observation matrix H is as follows:

[0052]

[0053] Where: h i Let represent the direction cosine vector between the i-th satellite and the BD / GPS terminal device, which can be calculated based on the initial position of the BD / GPS terminal device and the position of the i-th satellite. Let With r i =(x i ,y i ,z i Let represent the initial position coordinates of the terminal device and the position coordinates of the i-th satellite, respectively. In addition, h i It is a unit vector, i.e., ||h i ||=1.

[0054] S2. In BD / GPS terminal equipment, after introducing a weighting matrix, the expression for the weighted accuracy factor (WDOP) is as follows:

[0055]

[0056] Where tr(·) denotes the trace of a matrix, (·) T This indicates the matrix transpose; the weighted matrix W is a diagonal matrix, and its diagonal elements represent the weights of the i-th satellite.

[0057] S3. Calculating the WDOP value according to equation (2) requires multiple matrix inversions, which significantly increases the computational load. In addition, when the geometric observation matrix is ​​ill-conditioned, the matrix inversion process is prone to numerical instability problems.

[0058] Let N = H T If WH, then equation (2) can be transformed into:

[0059]

[0060] Where det(·) denotes the calculation of the determinant of a matrix, (·) * This indicates the calculation of the adjoint matrix.

[0061] S4. By calculating det(N) and tr(N) respectively * The WDOP value is obtained by calculating det(N) and tr(N) respectively. * ).

[0062] (1) Calculate the determinant det(N)

[0063] When three BD satellites and two GPS satellites are observed, since both the geometric observation matrix H and the weighting matrix W are square matrices, the determinant det(N) can be simplified to:

[0064] det(N) = det(H) T WH)=[det(H)] 2 det(W) (4)

[0065] Since the weighting matrix W is a diagonal matrix, its determinant can be simplified to:

[0066]

[0067] For the determinant det(H), after Schur decomposition and corresponding algebraic operations, we have:

[0068]

[0069] Combining equations (5) and (6), det(N) can be transformed into:

[0070]

[0071] (2) Calculate tr(N) * )

[0072] make Represents the adjoint matrix N * The i-th diagonal element then has

[0073]

[0074] Since N = H T WH=(n ij ) 5×5 Given a 5th-order real symmetric matrix, its adjoint matrix N is derived. * The diagonal elements are as follows:

[0075]

[0076]

[0077]

[0078]

[0079]

[0080] Substituting equations (9a)-(9e) into equation (8), then tr(N) * This can be transformed into:

[0081] tr(N * ) = n 11 A+n 12 B+n 13 C+n 14 D+n 15 E+n 22 F+n 23 G+n 24 H+n 25 I (10)

[0082] The correlation coefficients are as follows:

[0083]

[0084]

[0085]

[0086]

[0087]

[0088] as well as

[0089]

[0090] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Various modifications and variations can be made to the invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the invention should be included within the scope of the claims of the invention.

Claims

1. A method for directly calculating the weighted precision factor of a multi-system system, characterized in that, include: S1. When the BD / GPS terminal device observes 3 BD satellites and 2 GPS satellites, the geometric observation matrix H is obtained; S2. Based on the geometric observation matrix H and the weighting matrix W, the expression for the weighted accuracy factor is obtained: Where tr(·) denotes the trace of a matrix, (·) T This indicates the matrix transpose; the weighted matrix W is a diagonal matrix, and its diagonal elements represent the weights of the i-th satellite; S3, Let N = H T WH, transform the weighted precision factor expression in step S2 into: Where det(·) denotes the calculation of the determinant of a matrix, (·) * This indicates the calculation of the adjoint matrix; S4. By calculating det(N) and tr(N) respectively * ), to obtain the WDOP value; Step S4 describes the calculation of tr(N) * The process is as follows: B1, Order Represents the adjoint matrix N * The i-th diagonal element then has B2. Since N = H T WH=(n ij ) 5×5 Given a 5th-order real symmetric matrix, its adjoint matrix N is derived. * The diagonal elements are as follows: B3. According to step B2, tr(N) * Transformed into: tr(N * )=n 11 A+n 12 B+n 13 C+n 14 D+n 15 E+n 22 F+n 23 G+n 24 H+n 25 I; Where A, B, C, D, E, F, G, H, and I represent correlation coefficients; The formula for calculating the correlation coefficient in step B3 is: A=n 33 (n 44 n 55 -n 45 n 54 )+n 34 (n 45 n 53 -n 43 n 55 )+n 35 (n 43 n 54 -n 44 n 53 )+n 22 (n 44 n 55 -n 45 n 54 )+n 24 (n 45 n 52 -n 42 n 55 )+n 25 (n 42 n 54 -n 44 n 52 )+n 22 (n 33 n 55 -n 35 n 53 )+n 23 (n 35 n 52 -n 32 n 55 )+n 25 (n 32 n 53 -n 33 n 52 )+n 22 (n 33 n 44 -n 34 n 43 )+n 23 (n 34 n 42 -n 32 n 44 )+n 24 (n 32 n 43 -n 33 n 42 ) B=n 21 (n 45 n 54 -n 44 n 55 )+n 24 (n 41 n 55 -n 45 n 51 )+n 25 (n 44 n 51 -n 41 n 54 )+n 21 (n 35 n 53 -n 33 n 55 )+n 23 (n 31 n 55 -n 35 n 51 )+n 25 (n 33 n 51 -n 31 n 53 )+n 21 (n 34 n 43 -n 33 n 44 )+n 23 (n 31 n 44 -n 34 n 41 )+n 24 (n 33 n 41 -n 31 n 43 ) C=n 31 (n 45 n 54 -n 44 n 55 )+n 34 (n 41 n 55 -n 45 n 51 )+n 35 (n 44 n 51 -n 41 n 54 )+n 21 (n 32 n 55 -n 35 n 52 )+n 22 (n 35 n 51 -n 31 n 55 )+n 25 (n 31 n 52 -n 32 n 51 )+n 21 (n 32 n 44 -n 34 n 42 )+n 22 (n 34 n 41 -n 31 n 44 )+n 24 (n 31 n 42 -n 32 n 41 ) D=n 31 (n 43 n 55 -n 45 n 53 )+n 33 (n 45 n 51 -n 41 n 55 )+n 35 (n 41 n 53 -n 43 n 51 )+n 21 (n 42 n 55 -n 45 n 52 )+n 22 (n 45 n 51 -n 41 n 55 )+n 25 (n 41 n 52 -n 42 n 51 )+n 21 (n 33 n 42 -n 32 n 43 )+n 22 (n 31 n 43 -n 33 n 41 )+n 23 (n 32 n 41 -n 31 n 42 ) E=n 31 (n 44 n 53 -n 43 n 54 )+n 33 (n 41 n 54 -n 44 n 51 )+n 34 (n 43 n 51 -n 41 n 53 )+n 21 (n 44 n 52 -n 42 n 54 )+n 22 (n 41 n 54 -n 44 n 51 )+n 24 (n 42 n 51 -n 41 n 52 )+n 21 (n 33 n 52 -n 32 n 53 )+n 22 (n 31 n 53 -n 33 n 51 )+n 23 (n 32 n 51 -n 31 n 52 ) F=n 33 (n 44 n 55 -n 45 n 54 )+n 34 (n 45 n 53 -n 43 n 55 )+n 35 (n 43 n 54 -n 44 n 53 ) G=n 32 (n 45 n 54 -n 44 n 55 )+n 34 (n 42 n 55 -n 45 n 52 )+n 35 (n 44 n 52 -n 42 n 54 ) H=n 32 (n 42 n 55 -n 45 n 53 )+n 33 (n 45 n 52 -n 42 n 55 )+n 35 (n 42 n 53 -n 43 n 52 ) I=n 32 (n 43 n 54 +n 44 n 53 )+n 33 (n 42 n 54 -n 44 n 52 )+n 34 (n 43 n 52 -n 42 n 53 )。 2. The method for directly calculating the weighted precision factor of a multi-system system according to claim 1, characterized in that, The process of calculating det(N) in step S4 is as follows: A1. When 3 BD satellites and 2 GPS satellites are observed, since both the geometric observation matrix H and the weighting matrix W are square matrices, the determinant det(N) can be simplified to: it(N)=it(H T WH)=[that(H)] 2 it(W) A2. Since the weighting matrix W is a diagonal matrix, its determinant can be simplified to: A3. After Schur decomposition and corresponding algebraic operations, the determinant det(H) has the following properties: A4. Based on the expressions in steps A2 and A3, transform det(N) into:

Images