CN113466790B - Roland positioning calculation algorithm - Google Patents
Roland positioning calculation algorithm Download PDFInfo
- Publication number
- CN113466790B CN113466790B CN202110692526.0A CN202110692526A CN113466790B CN 113466790 B CN113466790 B CN 113466790B CN 202110692526 A CN202110692526 A CN 202110692526A CN 113466790 B CN113466790 B CN 113466790B
- Authority
- CN
- China
- Prior art keywords
- spherical
- calculating
- station
- profile
- correction term
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 38
- 238000000034 method Methods 0.000 claims abstract description 45
- 238000012937 correction Methods 0.000 claims abstract description 39
- 238000000572 ellipsometry Methods 0.000 claims abstract description 9
- 241001442234 Cosa Species 0.000 claims description 6
- 244000089409 Erythrina poeppigiana Species 0.000 claims description 6
- 235000009776 Rathbunia alamosensis Nutrition 0.000 claims description 6
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000005251 gamma ray Effects 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 239000002994 raw material Substances 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 abstract description 17
- 238000012360 testing method Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 4
- 230000006978 adaptation Effects 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
- G01S5/02—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using radio waves
- G01S5/14—Determining absolute distances from a plurality of spaced points of known location
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
- G01S5/02—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using radio waves
- G01S5/10—Position of receiver fixed by co-ordinating a plurality of position lines defined by path-difference measurements, e.g. omega or decca systems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
Abstract
The invention provides a Roland location calculation algorithm, which multiplies the measured time difference between a main station and a sub station by the propagation speed of light in the air and divides the time difference by the length of a long half axis of the earth to convert the time difference into the distance difference (radian value) between the main station and the sub stationCalculating the geodesic distance difference xi between the principal position and the primary and secondary stations by adopting an improved nested coefficient method or an improved ellipsometry i . With distance differences calculated from time differencesContinuously calculating spherical latitude correction term as referenceAnd spherical longitude correction term Deltalambda to adjust the profile so that the geodesic distance difference zeta between the profile and the primary and secondary stations i Is continuously close toWhen calculatedWhen delta lambda reaches a preset value, calculating the spherical approximate position coordinates of the receiverI.e. spherical coordinates of the receiverCompared with the existing Roland positioning calculation algorithm, the method can greatly improve the positioning accuracy, has high calculation speed and can be applied in a large range.
Description
Technical Field
The invention belongs to the field of radio navigation, and relates to a Roland positioning and resolving algorithm.
Background
The rowland location solution algorithm generally consists of two major steps, the first: the approximate position of the receiver, abbreviated as the approximate position, is calculated. When the approximate position is calculated, the earth surface is projected to the corresponding auxiliary spherical surface for calculation, and the calculated approximate position of the receiver is expressed as spherical coordinates; and a second step of: iterative receiver approximate position, iterative convergence to obtain receiver accurate position. Under the condition that the method for calculating the approximate position cannot continue to optimize, an iterative algorithm for optimizing the approximate position is selected to improve the accuracy of the Rowland positioning solution. In the approximate iteration process, the prior literature mostly adopts an Andoyer-Lambert formula to calculate the large ground wire distance. The Andoyer-Lambert formula has compact structure, symmetrical format and high real-time calculation speed, but the calculation accuracy is far smaller than that of an improved nested coefficient method and an improved macroellipsometry.
Therefore, in the approximate position iteration process, the method does not adopt an Andoyer-Lambert formula any more when calculating the large ground wire distance, and adopts an improved nested coefficient method and an improved large ellipsometry instead. The positioning accuracy of the improved Rowland positioning calculation algorithm is greatly improved, and the improved Rowland positioning calculation algorithm is fast in real-time calculation speed and can be widely applied.
Disclosure of Invention
The invention aims to provide a Roland positioning calculation algorithm which can perform high-precision positioning calculation, has high calculation speed and high calculation precision, and can be widely applied to the field of radio navigation.
The invention adopts the following technical scheme to realize the purposes:
the Roland location calculation algorithm comprises the following steps:
s1: input parameters and the geodetic coordinates of the receiver profile (B k ,L k )；
S2: calculating the large ground distance gamma from the updated position to each station by using improved nested coefficient method or improved ellipsometry i ；
S3: calculating and updating spherical azimuth angle A from approximate position to each station by adopting sine and cosine theorem of spherical triangle i Sine and cosine value sinA of (a) i And cosA i ；
S4: calculating and updating geodesic distance difference xi between main and auxiliary stations i ；
S5: calculating and updating longitude and latitude correction term of spherical surfaceAnd a correlation coefficient U of Deltalambda i And V i ；
S6: calculating and updating large ground line distance difference between main and auxiliary stations obtained by time differenceGeodesic distance difference ζ between primary and secondary stations derived from the profile i Difference of->
S7: calculating and updating spherical latitude correction term by adopting least square algorithmAnd spherical longitude correction term Deltalambda, and judging +.>And if the delta lambda reaches the preset value, executing the step 8, otherwise, returning to the step 2;
s8: calculated receiver spherical surface approximate position coordinatesNamely spherical coordinates of the receiver +.>
Further, the above-mentioned rowland location calculation algorithm is characterized in that the step S1 specifically includes:
the input parameters include: propagation time difference DeltaT between primary and secondary stations i Geodetic coordinates of stations involved in the Roland location solution (B i ,L i ) And the geodetic coordinates of the receiver profile (B k ,L k )；
Converting the geodetic coordinates of the station and the receiver profile into corresponding spherical coordinates, the spherical coordinates of the station being expressed asSpherical coordinates of the receiver profile are denoted +.>
Geodetic coordinates (B) i ,L i ) And spherical coordinatesThe conversion formula of (2) is:
λ i ＝L i
where f is the earth's flatness, a is the earth's major half-axis length, b is the earth's minor half-axis length, f=1-b/a.
Further, the update profile is calculated in the step S2 as follows:
λ k+1 ＝λ k +Δλ
wherein,spherical latitude representing the approximate position after the k+1th iteration, +.>Spherical latitude representing the k-th postamble, lambda k+1 And lambda (lambda) k Also indicated are the spherical longitudes of the profile after the k+1th and k iterations,/, respectively>Represents a spherical latitude correction term, Δλ represents a spherical longitude correction term, and +.>Δλ＝0。
Further, the method for improving the nested coefficients in the step S2 specifically includes:
s201: calculating the spherical surface warp difference delta lambda of two points of the general position and the station on the auxiliary spherical surface by iteratively calculating the geodetic warp difference correction term delta w; the geodetic coordinates of the profile and station are (B) k ,L k ) And (B) i ,L i )；
The relationship between the spherical warp difference Δλ and the geodetic warp difference Δl is:
Δλ＝ΔL+Δw
ΔL＝L i -L k
taking Δw=0 at the first iteration; if the values of the delta lambda and the delta L are not within [ -pi, pi ], adding or subtracting 2 pi to the values of the delta lambda and the delta L, and carrying out normalization processing to ensure that the values of the delta lambda and the delta L are within a normalized interval;
s202: calculating spherical angular distance sigma from position profile to station i ：
S203: calculating the critical point of the large ground wire of the overview position and the two stationsSpherical latitude>
S204: calculating the central angular distance of the large ground wire at two points of the general position and the station
S205: and calculating a geodetic correction term delta w:
K 3 ＝V[1+f+f 2 -V(3+7f-13V)]
if the earth warp difference correction term Δw does not reach the preset value, returning to execute step S202, and if the earth warp difference correction term Δw reaches the preset value, executing step S206;
s206: calculating the general position and the large ground wire distance gamma between two stations i ：
K 1 ＝1+t{1-t[3-t(5-11t)]/4}
K 2 ＝t{1-t[2-t(37-94t)/8]}
γ i ＝K 1 b(σ i -Δσ i )
Where e is the first eccentricity of the earth, a is the length of the earth's long half axis,
further, the modified macroellipsometry in the above step S2 is specifically:
s201: calculating the position of the profile and the earth's center latitude phi of the station k ，φ i ：
S202: calculating the spherical azimuth A of two points of the approximate position and the station k ，A i ：
Wherein Δl=l i -L k The method comprises the steps of carrying out a first treatment on the surface of the A refers to spherical azimuth:
s203: calculating the geocentric latitude theta of the outline position and the two stations in the cross section ellipse k ，θ i ：
Spherical surface-resolving right triangleThe method can obtain:
s204: short half shaft b for calculating cross section ellipse s ：
ρ represents the distance from a certain point P on the ellipse to the circle center O, φ represents the included angle between the OP connecting line and the x axis, and the standard equation of the meridian ellipse is:
wherein a is the length of the earth long half shaft, and b is the length of the earth short half shaft; the above is changed to polar coordinate representation, and
substituting the raw materials into the above formula, and finishing to obtain the product:
will bePoint geocentric latitude->Substituting the above, the minor half axis of the cross-section ellipse can be obtained>
Spherical surface-resolving right triangleThe method can obtain:
substituting the above intoIn the expression of (2), it is possible to obtain:
s205: calculating a cross-sectional ellipseIs of the first eccentricity e s ：
S206: calculating the length gamma of a large ellipse arc between two points i ；
The major ellipse arc length X on the cross-sectional ellipse can be expressed as:
X(θ)＝a(i 0 θ+i 2 sin 2θ+i 4 sin 4θ+i 6 sin 6θ)
wherein θ represents the geocentric latitude of a point on the cross-section ellipse; the coefficients in the formula are:
available, P 0 P k 、P 0 P i Is a length of arc: gamma ray 1 ＝X(θ k )、γ 2 ＝X(θ i ) Thus, the large elliptical arc length between two points, i.e. the large ground distance between the two points, is:
γ i ＝|γ 1 -γ 2 |＝|X(θ k )-X(θ i )|。
further, the step S3 specifically includes:
calculating a spherical azimuth A for updating the profile to the station i Sine value sinA of (a) i ：
Δλ＝λ i -λ k
Calculating a spherical azimuth A for updating the profile to the station i Cosine value cosA i ：
Further, the step S4 specifically includes:
m represents a main station, S represents a sub-station, gamma (M) represents a large ground distance from the main station to the sub-station, and gamma (S) represents a large ground distance from the main station to the sub-station; the large ground line distance difference between the main station and the auxiliary station calculated by the general position is expressed as follows:
ξ i ＝γ i (S)-γ(M)
large ground distance difference between main and auxiliary stations obtained by time differenceExpressed by Taylor formula and in general termsAnd (3) expanding:
further, the step S5 specifically includes:
computing moreNew spherical longitude and latitude correction termAnd a correlation coefficient U of Deltalambda i And V i ；
Wherein A is M Representing the spherical azimuth of the profile to the home station,representing the spherical azimuth of the profile to the secondary station.
Further, the step S7 specifically includes:
calculating and updating spherical latitude correction term by adopting least square algorithmAnd spherical longitude correction term Deltalambda, and judgeAnd if the delta lambda reaches the preset value, executing the step S8, otherwise, returning to the step S2;
order theΔλ＝λ-λ k . In combination with steps S4, S5, S6, the above formula can be written as:
if the receiver receives the propagation time difference delta T between the two main and auxiliary stations i The method can obtain:
the above formula can also be expressed as:
order theThe equation set of the above equation can be expressed as:
the least square method is adopted to solve the following problems:
further, the step S8 specifically includes:
calculated receiver approximate spherical coordinatesNamely spherical coordinates of the receiver +.>
λ＝λ k+1 。
The invention has the beneficial effects that:
in the invention, in the approximate position iteration process, the approximate position is adjusted by continuously calculating the spherical longitude and latitude correction term of the approximate position by taking the earth distance difference obtained by the known time difference as a reference, so that the earth distance difference between the main station and the auxiliary station obtained by the approximate position is continuously close to the reference line, and the calculated approximate position is the accurate position of the receiver when the spherical longitude and latitude correction term reaches a preset value. The invention greatly improves the positioning precision, has higher calculation speed and can be applied in a large range.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a schematic diagram of a modified macroellipsometry of the present invention;
FIG. 3 is a schematic view of a meridian ellipse according to the present invention;
fig. 4 is a schematic view of a spherical triangle of the present invention.
Detailed Description
Example embodiments will now be described more fully with reference to the accompanying drawings. However, the exemplary embodiments may be embodied in many forms and should not be construed as limited to the examples set forth herein; rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the concept of the example embodiments to those skilled in the art. The described features or characteristics may be combined in any suitable manner in one or more embodiments.
The invention will be described in detail with reference to the accompanying drawings, wherein the specific steps are as follows:
the principle of the invention is as follows: in the iteration process of the overlevel, the overlevel is adjusted by continuously calculating the spherical longitude and latitude correction term of the overlevel by taking the overlevel distance difference obtained by the known time difference as a reference, so that the overlevel distance difference between the main station and the auxiliary station obtained by the overlevel is continuously close to the reference line, and the calculated overlevel is the accurate position of the receiver when the spherical longitude and latitude correction term reaches a preset value.
The method is implemented according to the following steps:
step 1: input of the relevant parameters and the geodetic coordinates of the receiver profile (B k ,L k )；
The input relevant parameters include:propagation time difference DeltaT between primary and secondary stations i Geodetic coordinates of stations involved in the Roland location solution (B i ,L i ) And the geodetic coordinates of the receiver profile (B k ,L k ). Subsequent calculations require converting the geodetic coordinates of the station and receiver profile to corresponding spherical coordinates, the spherical coordinates of the station being expressed asSpherical coordinates of the receiver profile are denoted +.>
Geodetic coordinates (B) i ,L i ) And spherical coordinatesThe conversion formula of (2) is:
λ i ＝L i
where f is the earth's flatness, a is the earth's major half-axis length, b is the earth's minor half-axis length, f=1-b/a.
Step 2: calculating the large ground distance gamma from the updated position to each station by using improved nested coefficient method or improved ellipsometry i ；
First, an update profile is calculated:
λ k+1 ＝λ k +Δλ
wherein,spherical latitude representing the approximate position after the k+1th iteration, +.>Spherical latitude representing the k-th postamble, lambda k+1 And lambda (lambda) k Also indicated are the spherical longitudes of the profile after the k+1th and k iterations,/, respectively>Represents a spherical latitude correction term, Δλ represents a spherical longitude correction term, and +.>Δλ＝0。
Two algorithms for calculating the distance between large ground wires, namely an improved nested coefficient method and an improved macroellipsometry, are introduced into the invention. The present invention will specifically show the calculation step of the modified nested coefficient method in step 2.1, and the calculation step of the modified macroellipsis method in step 2.2. When the large ground wire distance is calculated, one algorithm is selected.
The specific steps of the improved nested coefficient method are as follows:
step 2.1.1: and calculating the spherical warp difference delta lambda of two points of the general position and the station on the auxiliary spherical surface by iteratively calculating the geodetic warp difference correction term delta w. The geodetic coordinates of the profile and station are (B) k ,L k ) And (B) i ,L i )。
The relationship between the spherical warp difference Δλ and the geodetic warp difference Δl is:
Δλ＝ΔL+Δw
ΔL＝L i -L k
at the first iteration, Δw=0 is taken. If the calculated values of Deltalambda and DeltaL are not within [ -pi, pi ] in the two formulas, the values of Deltalambda and DeltaL need to be added or subtracted by 2 pi, and normalization processing is carried out, so that the values of Deltalambda and DeltaL are within a normative interval.
Step 2.1.2: calculating spherical angular distance sigma from position profile to station i ：
Step 2.1.3: calculating the critical point of the large ground wire of the overview position and the two stationsSpherical latitude>
Step 2.1.4: calculating the central angular distance of the large ground wire at two points of the general position and the station
Step 2.1.5: and calculating a geodetic correction term delta w:
K 3 ＝V[1+f+f 2 -V(3+7f-13V)]
if the earth warp difference correction term Deltaw does not reach the preset value, returning to execute the step 2.1.2, and if the earth warp difference correction term Deltaw reaches the preset value, executing the step 2.1.6.
Step 2.1.6: calculating the general position and the large ground wire distance gamma between two stations i ：
K 1 ＝1+t{1-t[3-t(5-11t)]/4}
K 2 ＝t{1-t[2-t(37-94t)/8]}
γ i ＝K 1 b(σ i -Δσ i )
Where e is the first eccentricity of the earth, a is the length of the earth's long half axis,the computation of the improved nested coefficient method is completed.
The idea of the modified macroellipsometry is: the large ellipse arc length between two points is the large ground wire distance between the two points. The method comprises the following specific steps:
as shown in the modified ellipsometry schematic diagram of FIG. 2, the potential point is P k Station is P i . The right intersection point of a certain meridian on the elliptic surface and the elliptic surface of the cross section isThe geodetic latitude of the point is highest on the elliptical cross-section. P (P) 0 、Q 1 、Q 2 、Q n Is the positive intersection point of meridian and equator on ellipsoid. Phi (phi) k 、φ i 、/>Respectively represent P k 、P i 、/>The latitude of the centroid of the point A k 、A i Representing the spherical azimuth angle theta of the two points of the general position and the station k 、θ i The position of the profile, the latitude of the earth's center of the station in the cross-sectional ellipse is represented. Elliptical cross section refers to a passing ellipseIntersection of two known points on the sphere and the plane made by the center of the ellipsoid with the ellipsoid.
Step 2.2.1: calculating the position of the profile and the earth's center latitude phi of the station k ，φ i ：
Step 2.2.2: calculating the spherical azimuth A of two points of the approximate position and the station k ，A i ：
Wherein Δl=l i -L k . The quadrant determination principle of the spherical azimuth a is (a refers to any spherical azimuth):
step 2.2.3: calculating the geocentric latitude theta of the outline position and the two stations in the cross section ellipse k ，θ i ：
Spherical surface-resolving right triangleThe method can obtain:
step 2.2.4: short half shaft b for calculating cross section ellipse s ：
As shown in the schematic diagram of meridian ellipses in fig. 3, ρ represents the distance from a certain point P to the center O of the ellipse, and Φ represents the included angle between the OP connecting line and the x axis. The standard equation for a meridional ellipse is:
where a is the earth long half-shaft length and b is the earth short half-shaft length. The above is changed to polar coordinate representation, and
substituting the raw materials into the above formula, and finishing to obtain the product:
as shown in the modified macroellipsometry of fig. 2, in the meridional ellipse NOQ n In (1) willPoint geocentric latitude->Substituting the above, the minor half axis of the cross-section ellipse can be obtained>
Spherical surface-resolving right triangleThe method can obtain:
substituting the above intoIn the expression of (2), it is possible to obtain:
step 2.2.5: calculating a cross-sectional ellipseIs of the first eccentricity e s ：
Step 2.2.6: calculating the length gamma of a large ellipse arc between two points i ；
The major ellipse arc length X on the cross-sectional ellipse can be expressed as:
X(θ)＝a(i 0 θ+i 2 sin 2θ+i 4 sin 4θ+i 6 sin 6θ)
where θ represents the geocentric latitude of a point on the cross-sectional ellipse. The coefficients in the formula are:
available, P 0 P k 、P 0 P i Is a length of arc: gamma ray 1 ＝X(θ k )、γ 2 ＝X(θ i ) Thus, the large elliptical arc length between two points, i.e. the large ground distance between the two points, is:
γ i ＝|γ 1 -γ 2 |＝|X(θ k )-X(θ i )|
the calculation of the modified ellipsometry is completed and ended.
Step 3: calculating and updating spherical azimuth angle A from approximate position to each station by adopting sine and cosine theorem of spherical triangle i Sine and cosine value sinA of (a) i And cosA i ；
As shown in the spherical triangle schematic diagram of FIG. 4, the assumed point isThe station is->The north pole is N, the approximate location point P k The spherical distance to the N point is +.>Station P i The spherical distance to the N point is +.>d represents the spherical distance from the general point to the station. A is that i Representing the general site P k To station P i Is a spherical azimuth of a station, Δλ represents a spherical aberration between the potential point and the station, where Δλ=λ i -λ k 。
Calculating a spherical azimuth A for updating the profile to the station i Sine value sinA of (a) i ：
Calculating a spherical azimuth A for updating the profile to the station i Cosine value cosA i ：
Step 4: calculating and updating geodesic distance difference xi between main and auxiliary stations i ；
In the rowland location calculation, M represents a main station, S represents an auxiliary station, γ (M) represents a large ground distance from the main station to the main station, and γ (S) represents a large ground distance from the main station to the auxiliary station. Thus, the large ground distance difference between the primary and secondary stations calculated from the profile is expressed as:
ξ i ＝γ i (S)-γ(M)
large ground distance difference between main and auxiliary stations obtained by time differenceExpressed by Taylor formula and in general termsAnd (3) expanding:
step 5: calculating and updating longitude and latitude correction term of spherical surfaceAnd a correlation coefficient U of Deltalambda i And V i ；
Wherein the method comprises the steps of，A M Representing the spherical azimuth of the profile to the home station,representing the spherical azimuth of the profile to the secondary station.
Step 6: calculating and updating the geodesic distance difference xi between the master and slave stations obtained by the time difference t i Ground distance difference xi between main and auxiliary stations obtained by general position i Is a difference DeltaW of (1) i ；
Step 7: calculating and updating spherical latitude correction term by adopting least square algorithmAnd spherical longitude correction term Deltalambda, and judging +.>And if the delta lambda reaches the preset value, executing the step 8, otherwise, returning to the step 2;
order theΔλ＝λ-λ k . In connection with steps 4, 5, 6, the above formula can be written as:
if the receiver just receives the propagation time difference delta T between the two main and auxiliary stations i The method comprises the following steps of:
the above formula can also be expressed as:
order theThe equation set of the above equation can be expressed as:
the least square method is adopted to solve the following problems:
the method comprises the following steps: calculated receiver approximate spherical coordinatesNamely spherical coordinates of the receiver +.>
λ＝λ k+1
And (5) finishing calculation.
The following tables are the single-chain and double-chain cases (hyperbolic positioning is divided into single-chain positioning and multiple-chain positioning, and since multiple-chain positioning is converted into double-chain positioning after the preferred chain is performed during multiple-chain positioning, the positioning error of multiple-chain positioning solution is the positioning error of double-chain positioning solution), and the large ground wire distance algorithm is the positioning error of Andoyer-Lambert formula, improved nesting coefficient method and improved Roland positioning solution during macroellipsometry, respectively, and specific data are as follows:
single-unit chain positioning solution error with table-large ground wire distance algorithm being Andoyer-Lambert formula
/>
Single chain positioning calculating error by using table two large ground wire distance algorithm as improved nested coefficient method
/>
The three-major ground wire distance algorithm is a single chain positioning solution error for improving major ellipsometry
/>
Double-station chain positioning solution error with four-large ground wire distance algorithm of Andoyer-Lambert formula
The five-large ground wire distance algorithm is a double-station chain positioning solution error of an improved nested coefficient method
The six-major ground wire distance algorithm is a double-station chain positioning solution error for improving major ellipsometry
Examples
Under the condition of a single station chain, the geodetic coordinates of the approximate positions of test points (22 DEG 40',117 DEG 30') are input, the test station chain is an east sea (8390) station chain, the main station is Xuan Cheng, the auxiliary stations are Rongcheng and Fuping, and the propagation time difference TD between Rongcheng and Rongcheng is obtained through simulation 1 2.4030269e-03 seconds, propagation time difference TD between Xuancheng and Xuanping 2 For-2.6955252 e-03 seconds, the exact coordinates of the test points and positioning errors are calculated.
According to the implementation of the method steps of the invention, the ellipsoid parameters in the experiment are as follows: the earth long half shaft a is 6378140 m, and the earth short half shaft b is 6356755.288158 m.
Under the same condition, the accurate coordinates and positioning error data of the test point calculated by the method and the existing positioning calculation method for calculating the large ground wire distance by using the Andoyer-Lambert formula are shown in a seventh table. As can be seen from the seventh table, the positioning error of the test point calculated by the algorithm is greatly reduced compared with that of the test point calculated by the existing positioning calculation method for calculating the large ground wire distance by using the Andoyer-Lambert formula, and the positioning accuracy is greatly improved. Thus, the high-precision characteristic of the method of the invention is verified.
Positioning errors (unit: meter) of seven different positioning calculation algorithms
Other embodiments of the invention will be apparent to those skilled in the art from consideration of the specification and practice of the invention disclosed herein. This application is intended to cover any variations, uses, or adaptations of the invention following, in general, the principles of the invention and including such departures from the present disclosure as come within known or customary practice within the art to which the invention pertains. It is intended that the specification and examples be considered as exemplary only, with a true scope and spirit of the invention being indicated by the following claims.
Claims (10)
1. A rowland location solution algorithm, comprising the steps of:
s1: input parameters and the geodetic coordinates of the receiver profile (B k ,L k )；
S2: calculating the large ground distance gamma from the updated position to each station by using improved nested coefficient method or improved ellipsometry i ；
S3: calculating and updating spherical azimuth angle A from approximate position to each station by adopting sine and cosine theorem of spherical triangle i Sine and cosine value sinA of (a) i And cosA i ；
S4: calculating and updating geodesic distance difference xi between main and auxiliary stations i ；
S5: calculating and updating longitude and latitude correction term of spherical surfaceAnd a correlation coefficient U of Deltalambda i And V i ；
S6: calculating and updating large ground line distance difference between main and auxiliary stations obtained by time differenceGeodesic distance difference ζ between primary and secondary stations derived from the profile i Difference of->
S7: calculating and updating spherical latitude correction term by adopting least square algorithmAnd spherical longitude correction term Deltalambda, and judging +.>And if the delta lambda reaches the preset value, executing the step 8, otherwise, returning to the step 2;
s8: calculated receiver spherical surface approximate position coordinatesNamely spherical coordinates of the receiver +.>
2. The rowland location calculation algorithm according to claim 1, wherein the step S1 is specifically:
the input parameters include: propagation time difference DeltaT between primary and secondary stations i Geodetic coordinates of stations involved in the Roland location solution (B i ,L i ) And the geodetic coordinates of the receiver profile (B k ,L k )；
Converting the geodetic coordinates of the station and the receiver profile into corresponding spherical coordinates, the spherical coordinates of the station being expressed asSpherical coordinates of the receiver profile are denoted +.>
Geodetic coordinates (B) i ,L i ) And spherical coordinatesThe conversion formula of (2) is:
λ i ＝L i
where f is the earth's flatness, a is the earth's major half-axis length, b is the earth's minor half-axis length, f=1-b/a.
3. The rowland location calculation algorithm according to claim 2, wherein the update profile is calculated in step S2 as follows:
λ k+1 ＝λ k +Δλ
wherein,spherical latitude representing the approximate position after the k+1th iteration, +.>Spherical latitude representing the k-th postamble, lambda k+1 And lambda (lambda) k Also indicated are the spherical longitudes of the profile after the k+1th and k iterations,/, respectively>Represents a spherical latitude correction term, Δλ represents a spherical longitude correction term, and +.>Δλ＝0。
4. A rowland location calculation algorithm according to claim 3, wherein said improved nested coefficient method in step S2 is specifically:
s201: calculating the spherical surface warp difference delta lambda of two points of the general position and the station on the auxiliary spherical surface by iteratively calculating the geodetic warp difference correction term delta w; the geodetic coordinates of the profile and station are (B) k ,L k ) And (B) i ,L i )；
The relationship between the spherical warp difference Δλ and the geodetic warp difference Δl is:
Δλ＝ΔL+Δw
ΔL＝L i -L k
taking Δw=0 at the first iteration; if the values of the delta lambda and the delta L are not within [ -pi, pi ], adding or subtracting 2 pi to the values of the delta lambda and the delta L, and carrying out normalization processing to ensure that the values of the delta lambda and the delta L are within a normalized interval;
s202: calculating spherical angular distance sigma from position profile to station i ：
S203: calculating the critical point of the large ground wire of the overview position and the two stationsSpherical latitude>
S204: calculating the central angular distance of the large ground wire at two points of the general position and the station
S205: and calculating a geodetic correction term delta w:
K 3 ＝V[1+f+f 2 -V(3+7f-13V)]
if the earth warp difference correction term Δw does not reach the preset value, returning to execute step S202, and if the earth warp difference correction term Δw reaches the preset value, executing step S206;
s206: calculating the general position and the large ground wire distance gamma between two stations i ：
K 1 ＝1+t{1-t[3-t(5-11t)]/4}
K 2 ＝t{1-t[2-t(37-94t)/8]}
γ i ＝K 1 b(σ i -Δσ i )
Where e is the first eccentricity of the earth, a is the length of the earth's long half axis,
5. a rowland location calculation algorithm according to claim 3, wherein said modified macroellipses in step S2 are specifically:
s201: calculating the position of the profile and the earth's center latitude phi of the station k ，φ i ：
S202: calculating the spherical azimuth A of two points of the approximate position and the station k ，A i ：
Wherein Δl=l i -L k The method comprises the steps of carrying out a first treatment on the surface of the A refers to spherical azimuth:
s203: calculating the geocentric latitude theta of the outline position and the two stations in the cross section ellipse k ，θ i ：
Spherical surface-resolving right triangleThe method can obtain:
s204: short half shaft b for calculating cross section ellipse s ：
ρ represents the distance from a certain point P on the ellipse to the circle center O, φ represents the included angle between the OP connecting line and the x axis, and the standard equation of the meridian ellipse is:
wherein a is the length of the earth long half shaft, and b is the length of the earth short half shaft; the above is changed to polar coordinate representation, and
substituting the raw materials into the above formula, and finishing to obtain the product:
will bePoint geocentric latitude->Substituting the above, the minor half axis of the cross-section ellipse can be obtained>
Spherical surface-resolving right triangleThe method can obtain:
substituting the above intoIn the expression of (2), it is possible to obtain:
s205: calculating a cross-sectional ellipseIs of the first eccentricity e s ：
S206: calculating the length gamma of a large ellipse arc between two points i ；
The major ellipse arc length X on the cross-sectional ellipse can be expressed as:
X(θ)＝a(i 0 θ+i 2 sin2θ+i 4 sin4θ+i 6 sin6θ)
wherein θ represents the geocentric latitude of a point on the cross-section ellipse; the coefficients in the formula are:
available, P 0 P k 、P 0 P i Is a length of arc: gamma ray 1 ＝X(θ k )、γ 2 ＝X(θ i ) Thus, the large elliptical arc length between two points, i.e. the large ground distance between the two points, is:
γ i ＝|γ 1 -γ 2 |＝|X(θ k )-X(θ i )|。
6. the rowland location calculation algorithm according to claim 1, wherein the step S3 is specifically:
calculating a spherical azimuth A for updating the profile to the station i Sine value sinA of (a) i ：
Δλ＝λ i -λ k
Calculating a spherical azimuth A for updating the profile to the station i Cosine value cosA i ：
Δλ＝λ i -λ k 。
7. The rowland location calculation algorithm according to claim 1, wherein said step S4 is specifically:
m represents a main station, S represents a sub-station, gamma (M) represents a large ground distance from the main station to the sub-station, and gamma (S) represents a large ground distance from the main station to the sub-station; the large ground line distance difference between the main station and the auxiliary station calculated by the general position is expressed as follows:
ξ i ＝γ i (S)-γ(M)
large ground distance difference between main and auxiliary stations obtained by time differenceExpressed by Taylor's formula and expressed in the general position +.>And (3) expanding:
8. the rowland location calculation algorithm according to claim 1, wherein said step S5 is specifically:
calculating and updating longitude and latitude correction term of spherical surfaceAnd a correlation coefficient U of Deltalambda i And V i ；
Wherein A is M Representing the spherical azimuth of the profile to the home station,representing the spherical azimuth of the profile to the secondary station.
9. The rowland location calculation algorithm according to claim 1, wherein the step S7 is specifically:
calculating and updating spherical latitude correction term by adopting least square algorithmAnd spherical longitude correction term Deltalambda, and judging +.>And if the delta lambda reaches the preset value, executing the step S8, otherwise, returning to the step S2;
order theΔλ＝λ-λ k The method comprises the steps of carrying out a first treatment on the surface of the In combination with steps S4, S5, S6, the above formula can be written as:
if the receiver receives the propagation time difference delta T between the two main and auxiliary stations i The method can obtain:
the above formula can also be expressed as:
order theThe equation set of the above equation can be expressed as:
the least square method is adopted to solve the following problems:
10. the rowland location calculation algorithm according to claim 1, wherein said step S8 is specifically:
calculated receiver approximate spherical coordinatesNamely spherical coordinates of the receiver +.>
λ＝λ k+1 。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110692526.0A CN113466790B (en) | 2021-06-22 | 2021-06-22 | Roland positioning calculation algorithm |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110692526.0A CN113466790B (en) | 2021-06-22 | 2021-06-22 | Roland positioning calculation algorithm |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113466790A CN113466790A (en) | 2021-10-01 |
| CN113466790B true CN113466790B (en) | 2024-03-01 |
Family
ID=77869147
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110692526.0A Active CN113466790B (en) | 2021-06-22 | 2021-06-22 | Roland positioning calculation algorithm |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113466790B (en) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5017926A (en) * | 1989-12-05 | 1991-05-21 | Qualcomm, Inc. | Dual satellite navigation system |
| JP2000088585A (en) * | 1998-09-14 | 2000-03-31 | Kenwood Corp | Distance calculation device and position calculation device |
| CN106792516A (en) * | 2016-12-02 | 2017-05-31 | 武汉理工大学 | 3-D positioning method based on radio communication base station |
| CN112198537A (en) * | 2020-09-22 | 2021-01-08 | 中国科学院国家授时中心 | Rowland high-precision positioning resolving method based on difference |
-
2021
- 2021-06-22 CN CN202110692526.0A patent/CN113466790B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5017926A (en) * | 1989-12-05 | 1991-05-21 | Qualcomm, Inc. | Dual satellite navigation system |
| JP2000088585A (en) * | 1998-09-14 | 2000-03-31 | Kenwood Corp | Distance calculation device and position calculation device |
| CN106792516A (en) * | 2016-12-02 | 2017-05-31 | 武汉理工大学 | 3-D positioning method based on radio communication base station |
| CN112198537A (en) * | 2020-09-22 | 2021-01-08 | 中国科学院国家授时中心 | Rowland high-precision positioning resolving method based on difference |
Non-Patent Citations (1)
| Title |
|---|
| 燕保荣 ; 李云 ; 郭伟 ; 华宇 ; .一种基于伪距的罗兰C定位授时解算方法.时间频率学报.2020,(第02期),48-60. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113466790A (en) | 2021-10-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109708629B (en) | Aircraft cluster collaborative navigation method for performance condition of differential positioning | |
| CN105203104B (en) | A kind of gravitational field modeling method suitable for high accuracy inertial navigation system | |
| Wang et al. | Convex combination based target localization with noisy angle of arrival measurements | |
| CN109959898B (en) | Self-calibration method for base type underwater sound passive positioning array | |
| Li et al. | A profile-matching method for wireless positioning | |
| CN111397599A (en) | Improved ICCP (Integrated Circuit chip) underwater geomagnetic matching method based on triangular matching algorithm | |
| CN114061591B (en) | Contour line matching method based on sliding window data backtracking | |
| CN111380518A (en) | SINS/USBL tight combination navigation positioning method introducing radial velocity | |
| CN112540371A (en) | Near-bottom multi-beam coordinate conversion processing method | |
| CN110736484B (en) | Background magnetic field calibration method based on fusion of gyroscope and magnetic sensor | |
| Liu et al. | Tightly coupled modeling and reliable fusion strategy for polarization-based attitude and heading reference system | |
| CN113466790B (en) | Roland positioning calculation algorithm | |
| CN108917698B (en) | Azimuth angle calculation method | |
| CN113189541B (en) | Positioning method, device and equipment | |
| RU2666360C1 (en) | Target coordinates determining method and system in the “request-response” system | |
| Liu et al. | A moving source localization method for distributed passive sensor using TDOA and FDOA measurements | |
| CN107729705B (en) | Method for measuring and calculating precision of single panel of surface antenna | |
| CN116299163A (en) | Unmanned aerial vehicle track planning method, unmanned aerial vehicle track planning device, unmanned aerial vehicle track planning equipment and unmanned aerial vehicle track planning medium | |
| CN113503891B (en) | SINSDVL alignment correction method, system, medium and equipment | |
| CN111984917B (en) | Calculation method of turning center in ball great circle track turning process | |
| CN111479214B (en) | Wireless sensor network optimal target positioning method based on TOA measurement | |
| CN110687502B (en) | Short wave direction finding data set labeling method based on least square positioning | |
| CN109506645B (en) | Star sensor mounting matrix ground accurate measurement method | |
| Dong et al. | Indoor robot localization combining feature clustering with wireless sensor network | |
| CN108387897B (en) | Projectile body positioning method based on improved Gauss Newton-genetic hybrid algorithm |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |