阅读说明：本技术 一种罗兰定位解算算法 (Rowland positioning calculation algorithm ) 是由 刘江凡 张碧瑶 席晓莉 于 2021-06-22 设计创作，主要内容包括：本发明提供了一种罗兰定位解算算法,将测量到的主副台之间的时间差,乘以光在空气中传播的速度并除以地球长半轴长度转换为主副台之间的距离差(弧度值)采用改进嵌套系数法或者改进大椭圆法计算概位到主副台之间的大地线距离差ξ-i。以由时间差计算得到的距离差为基准,不断计算球面纬度修正项和球面经度修正项Δλ调整概位,使得概位到主副台之间的大地线距离差ξ-i不断靠近当计算的和Δλ达到预设值时,计算的接收机球面概位坐标即接收机的球面坐标本发明与现有的罗兰定位解算算法相比,可以大幅度提高定位精度,计算速度较快,可以大范围内进行应用。(The invention provides a Rowland positioning calculation algorithm, which converts the measured time difference between a main platform and an auxiliary platform into the distance difference (radian value) between the main platform and the auxiliary platform by multiplying the speed of light propagating in air and dividing the length of a long half shaft of the earth by the length of the long half shaft of the earth Calculating the distance difference xi between the probable position and the main and auxiliary stations by using an improved nesting coefficient method or an improved large ellipse method i . By distance differences calculated from the time differences Continuously calculating a spherical latitude correction term as a reference And the spherical longitude correction term delta lambda adjusts the approximate position to make the distance difference xi between the approximate position and the main station and the secondary station i Is constantly close to When calculated When the sum delta lambda reaches a preset value, the calculated receiver spherical approximate position coordinates I.e. spherical coordinates of the receiver Compared with the existing Rowland positioning calculation algorithm, the method can greatly improve the positioning accuracy, has higher calculation speed, and can be applied in a large range.)

1. A roland positioning calculation algorithm, characterized by comprising the steps of:

s1: geodetic coordinates (B) of input parameters and receiver approximate bits_k,L_k)；

S2: calculating the geodesic distance gamma from the updated probable bit to each station by using an improved nesting coefficient method or an improved large ellipse method_i；

S3: spherical azimuth angles A from the approximate position to each station are calculated and updated by adopting sine and cosine theorem of spherical triangle_iSine and cosine value sinA of_iAnd cosA_i；

S4: calculating and updating the earth between the main station and the secondary stationDifference xi of line distance_i；

S5: calculating and updating spherical longitude and latitude correction termCorrelation coefficient U of sum Delta lambda_iAnd V_i；

S6: calculating and updating large ground distance difference between primary and secondary stations obtained by time differenceLarge earth distance difference xi between primary and secondary stations obtained from the bit-of-probability_iDifference of (2)

S7: calculating and updating spherical latitude correction term by adopting least square algorithmAnd a spherical longitude correction term Delta lambda, and judgingAnd whether the delta lambda reaches a preset value, if so, executing the step 8, otherwise, returning to the step 2;

s8: calculated receiver spherical approximate position coordinatesI.e. spherical coordinates of the receiver

2. The rowland positioning 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_iGeodetic coordinates of stations participating in the Rowland positioning solution (B_i,L_i) And geodetic coordinates (B) of the receiver profile_k,L_k)；

Converting geodetic coordinates of the station and receiver profile into corresponding spherical coordinates, the spherical coordinates of the station being expressed asSpherical coordinates of receiver profile are expressed as

Geodetic coordinates (B)_i,L_i) Coordinates with spherical surfaceThe conversion formula of (1) is as follows:

λ_i＝L_i

wherein f is the oblateness of the earth, a is the length of the longer half axis of the earth, b is the length of the shorter half axis of the earth, and f is 1-b/a.

3. The roland positioning calculation algorithm according to claim 2 wherein the calculation of the update profile in step S2 is as follows:

λ_k+1＝λ_k+Δλ

wherein the content of the first and second substances,represents the spherical latitude of the probable location after the (k + 1) th iteration,spherical latitude representing the k-th after-probable position, and, similarly, λ_k+1And λ_kAlso denoted spherical longitudes for the approximate after the k +1 th and k-th iterations respectively,representing a spherical latitude correction term, and delta lambda representing a spherical longitude correction term, and taking the correction term when iterating for the first timeΔλ＝0。

4. The Rowland positioning calculation algorithm according to claim 3, wherein the improved nesting coefficient method in step S2 is specifically:

s201: calculating the spherical warp difference delta lambda of the approximate position and the station on the auxiliary circular spherical surface by iteratively calculating the earth warp difference correction item delta w; the geodetic coordinates of the approximate and station are respectively (B)_k,L_k) And (B)_i,L_i)；

The relationship between the spherical warp difference Δ λ and the earth warp difference Δ L is:

Δλ＝ΔL+Δw

ΔL＝L_i-L_k

during the first iteration, taking delta w as 0; if the values of the delta lambda and the delta L are not in [ -pi, pi ], adding or subtracting 2 pi to the values of the delta lambda and the delta L, and carrying out normalized processing to make the values of the delta lambda and the delta L in a normalized interval;

s202: calculating the spherical angular distance sigma from the probable position to the station_i：

S203: calculating the critical point of the large ground wire passing through the general position and the stationSpherical latitude of

S204: calculating the central angular distance between the earth wires of the general position and the station

S205: calculating an earth warp difference correction term Δ w:

K₃＝V[1+f+f²-V(3+7f-13V)]

if the earth difference correction item delta w does not reach the preset value, returning to execute the step S202, and if the earth difference correction item delta w reaches the preset value, executing the step S206;

s206: calculating the large ground wire distance gamma between the two points of the approximate position and the station_i：

K₁＝1+t{1-t[3-t(5-11t)]/4}

K₂＝t{1-t[2-t(37-94t)/8]}

γ_i＝K₁b(σ_i-Δσ_i)

Wherein e is the first eccentricity of the earth, a is the length of the longer half axis of the earth,

5. the Rowland positioning calculation algorithm according to claim 3, wherein the improved ellipse method in the step S2 is specifically:

s201: calculating geocentric latitude phi of the probable location and the station_k，φ_i：

S202: calculating the spherical azimuth angle A of the approximate station and the station_k，A_i：

Wherein Δ L ═ L_i-L_k(ii) a A is the spherical azimuth:

s203: calculating geocentric latitude theta of the approximate position and the station in the section ellipse_k，θ_i：

Spherical relief surface right-angled triangleThe following can be obtained:

s204: short semi-axis b for calculating section ellipse_s：

Rho represents the distance from a certain point P to the center O of the ellipse, phi represents the included angle between the OP connecting line and the x axis, and the standard equation of the meridian ellipse is as follows:

wherein, a is the length of the earth's major semi-axis, b is the length of the earth's minor semi-axis; the above formula is expressed by polar coordinates, and

substituting the formula into the formula and finishing to obtain:

will be provided withGeocentric latitude of pointSubstituting the above formula to obtain a semi-minor axis with an oval cross section

Spherical relief surface right-angled triangleThe following can be obtained:

substituting the above formula intoIn the expression of (a), it can be:

s205: calculating the cross-sectional ellipseFirst eccentricity e_s：

S206: calculating the large ellipse arc length gamma between two points_i；

The large ellipse arc length X on the cross-sectional ellipse can be expressed as:

X(θ)＝a(i₀θ+i₂ sin2θ+i₄ sin4θ+i₆ sin6θ)

wherein theta represents the geocentric latitude of a certain point on the section ellipse; the coefficients in the formula are:

can be obtained, P₀P_k、P₀P_iArc length of (d): gamma ray₁＝X(θ_k)、γ₂＝X(θ_i) Thus, the large elliptical arc length between two points, i.e., the geodetic distance between two points, is:

γ_i＝|γ₁-γ₂|＝|X(θ_k)-X(θ_i)|。

6. the rowland positioning calculation algorithm according to claim 1, wherein the step S3 is specifically:

calculating and updating spherical azimuth A of probable position to station_iSine value sinA of_i：

Δλ＝λ_i-λ_k

Calculating and updating spherical azimuth A of probable position to station_iCosine value of (cosA)_i：

Δλ＝λ_i-λ_k。

7. The rowland positioning calculation algorithm according to claim 1, wherein the step S4 is specifically:

using M to represent the main station and S to represent the secondary station, then gamma (M) represents the distance between the probable bit and the main station, and gamma (S) represents the distance between the probable bit and the secondary station; the large difference in ground distance between the primary and secondary stations calculated from the bit-profile is represented as:

ξ_i＝γ_i(S)-γ(M)

large difference in distance of ground lines between primary and secondary stations to be obtained from time differenceExpressed by Taylor's formula and in the approximate positionUnfolding:

8. the rowland positioning calculation algorithm according to claim 1, wherein the step S5 is specifically:

calculating and updating spherical longitude and latitude correction termCorrelation coefficient U of sum Delta lambda_iAnd V_i；

Wherein A is_MIndicating the spherical azimuth of the overview to the main stage,indicating the spherical azimuth of the overview to the secondary stage.

9. The rowland positioning calculation algorithm according to claim 1, wherein the step S7 is specifically:

calculating and updating spherical latitude correction term by adopting least square algorithmAnd a spherical longitude correction term Delta lambda, and judgingAnd whether the delta lambda reaches a preset value, if so, executing the step S8, otherwise, returning to the step S2;

order toΔλ＝λ-λ_k(ii) a In conjunction with steps S4, S5, S6, the above equation may be written as:

if the receiver receives the propagation time difference Delta T between the two main stations and the secondary station_iThe following can be obtained:

the above formula can also be expressed as:

order toThe system of equations of the above equation can be expressed as:

solving by adopting a least square method to obtain:

10. the rowland positioning calculation algorithm according to claim 1, wherein the step S8 is specifically:

calculated receiver approximate position spherical coordinateI.e. spherical coordinates of the receiver

λ＝λ_k+1。

Technical Field

The invention belongs to the field of radio navigation, and relates to a Rowland positioning calculation algorithm.

Background

The rowland positioning solution algorithm is generally divided into two major steps, the first step: the approximate position of the receiver, abbreviated as approximate bits, is calculated. When calculating the approximate position, projecting the earth surface to the corresponding auxiliary spherical surface for calculation, and expressing the calculated receiver approximate position by spherical coordinates; the second step is that: and iterating the approximate bit of the receiver, and iteratively converging to obtain the accurate position of the receiver. And under the condition that the method for calculating the probable bits cannot be optimized continuously, selecting an iterative algorithm for optimizing the probable bits to improve the accuracy of the Rowland positioning calculation. In the general bit iteration process, the prior documents mostly adopt an Andoyer-Lambert formula to calculate the geodesic distance. The Andoyer-Lambert formula has the advantages of compact structure, symmetrical format and high real-time calculation speed, but the calculation precision is far smaller than that of an improved nesting coefficient method and an improved large ellipse method.

Therefore, in the general bit iteration process, the method of the invention does not adopt the Andoyer-Lambert formula any more when calculating the geodetic distance, and adopts an improved nesting coefficient method and an improved ellipse method instead. The improved Rowland positioning calculation algorithm has greatly improved positioning accuracy, high real-time calculation speed and wide application.

Disclosure of Invention

The invention aims to provide a Rowland 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.

In order to achieve the purpose, the invention adopts the following technical scheme:

the Rowland positioning calculation algorithm comprises the following steps:

s1: geodetic coordinates (B) of input parameters and receiver approximate bits_k,L_k)；

S2: calculating the geodesic distance gamma from the updated probable bit to each station by using an improved nesting coefficient method or an improved large ellipse method_i；

S3: spherical azimuth angles A from the approximate position to each station are calculated and updated by adopting sine and cosine theorem of spherical triangle_iSine and cosine value sinA of_iAnd cosA_i；

S4: calculating and updating large earth wire distance difference xi between main station and secondary station_i；

S5: calculating and updating spherical longitude and latitude correction termCorrelation coefficient U of sum Delta lambda_iAnd V_i；

S6: calculating and updating large ground distance difference between primary and secondary stations obtained by time differenceLarge earth distance difference xi between primary and secondary stations obtained from the bit-of-probability_iDifference of (2)

S7: calculating and updating spherical latitude correction term by adopting least square algorithmAnd a spherical longitude correction term Delta lambda, and judgingAnd whether the delta lambda reaches a preset value, if so, executing the step 8, otherwise, returning to the step 2;

s8: calculated receiver spherical approximate position coordinatesI.e. spherical coordinates of the receiver

Further, the above-mentioned rowland positioning calculation algorithm is characterized in that the step S1 specifically includes:

the input parameters include: propagation time difference DeltaT between primary and secondary stations_iEarth coordinates (B) of a station participating in Rowland positioning solution_i,L_i) And geodetic coordinates (B) of the receiver profile_k,L_k)；

Converting geodetic coordinates of the station and receiver profile into corresponding spherical coordinates, the spherical coordinates of the station being expressed asSpherical coordinates of receiver profile are expressed as

Geodetic coordinates (B)_i,L_i) Coordinates with spherical surfaceThe conversion formula of (1) is as follows:

λ_i＝L_i

wherein f is the oblateness of the earth, a is the length of the longer half axis of the earth, b is the length of the shorter half axis of the earth, and f is 1-b/a.

Further, the update approximate bits are calculated in step S2 as follows:

λ_k+1＝λ_k+Δλ

wherein the content of the first and second substances,represents the spherical latitude of the probable location after the (k + 1) th iteration,spherical latitude representing the k-th after-probable position, and, similarly, λ_k+1And λ_kAlso denoted spherical longitudes for the approximate after the k +1 th and k-th iterations respectively,representing a spherical latitude correction term, and delta lambda representing a spherical longitude correction term, and taking the correction term when iterating for the first timeΔλ＝0。

Further, the method for improving the nesting coefficient in step S2 specifically includes:

s201: calculating the spherical warp difference delta lambda of the approximate position and the station on the auxiliary circular spherical surface by iteratively calculating the earth warp difference correction item delta w; the geodetic coordinates of the approximate and station are respectively (B)_k,L_k) And (B)_i,L_i)；

The relationship between the spherical warp difference Δ λ and the earth warp difference Δ L is:

Δλ＝ΔL+Δw

ΔL＝L_i-L_k

during the first iteration, taking delta w as 0; if the values of the delta lambda and the delta L are not in [ -pi, pi ], adding or subtracting 2 pi to the values of the delta lambda and the delta L, and carrying out normalized processing to make the values of the delta lambda and the delta L in a normalized interval;

s202: calculating the spherical angular distance sigma from the probable position to the station_i：

S203: calculated ruleCritical point of large earth wire between station and stationSpherical latitude of

S204: calculating the central angular distance between the earth wires of the general position and the station

S205: calculating an earth warp difference correction term Δ w:

K₃＝V[1+f+f²-V(3+7f-13V)]

if the earth difference correction item delta w does not reach the preset value, returning to execute the step S202, and if the earth difference correction item delta w reaches the preset value, executing the step S206;

s206: calculating the large ground wire distance gamma between the two points of the approximate position and the station_i：

K₁＝1+t{1-t[3-t(5-11t)]/4}

K₂＝t{1-t[2-t(37-94t)/8]}

γ_i＝K₁b(σ_i-Δσ_i)

Wherein e is the first eccentricity of the earth, a is the length of the longer half axis of the earth,

further, the modified ellipse method in step S2 specifically includes:

s201: calculating geocentric latitude phi of the probable location and the station_k，φ_i：

S202: calculating the spherical azimuth angle A of the approximate station and the station_k，A_i：

Wherein Δ L ═ L_i-L_k(ii) a A is the spherical azimuth:

s203: calculating geocentric latitude theta of the approximate position and the station in the section ellipse_k，θ_i：

Spherical relief surface right-angled triangleThe following can be obtained:

s204: short semi-axis b for calculating section ellipse_s：

Rho represents the distance from a certain point P to the center O of the ellipse, phi represents the included angle between the OP connecting line and the x axis, and the standard equation of the meridian ellipse is as follows:

wherein, a is the length of the earth's major semi-axis, b is the length of the earth's minor semi-axis; the above formula is expressed by polar coordinates, and

substituting the formula into the formula and finishing to obtain:

will be provided withGeocentric latitude of pointSubstituting the above formula to obtain a semi-minor axis with an oval cross section

Spherical relief surface right-angled triangleThe following can be obtained:

substituting the above formula intoIn the expression of (a), it can be:

s205: calculating the cross-sectional ellipseFirst eccentricity e_s：

S206: calculating the large ellipse arc length gamma between two points_i；

The large ellipse arc length X on the cross-sectional ellipse can be expressed as:

X(θ)＝a(i₀θ+i₂ sin 2θ+i₄ sin 4θ+i₆ sin 6θ)

wherein theta represents the geocentric latitude of a certain point on the section ellipse; the coefficients in the formula are:

can be obtained, P₀P_k、P₀P_iArc length of (d): gamma ray₁＝X(θ_k)、γ₂＝X(θ_i) Thus, the large elliptical arc length between two points, i.e., the geodetic distance between two points, is:

γ_i＝|γ₁-γ₂|＝|X(θ_k)-X(θ_i)|。

further, the step S3 is specifically:

calculating and updating spherical azimuth A of probable position to station_iSine value sinA of_i：

Δλ＝λ_i-λ_k

Calculating and updating spherical azimuth A of probable position to station_iCosine value of (cosA)_i：

Further, the step S4 is specifically:

using M to represent the main station and S to represent the secondary station, then gamma (M) represents the distance between the probable bit and the main station, and gamma (S) represents the distance between the probable bit and the secondary station; the large difference in ground distance between the primary and secondary stations calculated from the bit-profile is represented as:

ξ_i＝γ_i(S)-γ(M)

large difference in distance of ground lines between primary and secondary stations to be obtained from time differenceExpressed by Taylor's formula and in the approximate positionUnfolding:

further, the step S5 is specifically:

calculating and updating spherical longitude and latitude correction termCorrelation coefficient U of sum Delta lambda_iAnd V_i；

Wherein A is_MIndicating the spherical azimuth of the overview to the main stage,indicating the spherical azimuth of the overview to the secondary stage.

Further, the step S7 is specifically:

calculating and updating spherical latitude correction term by adopting least square algorithmAnd a spherical longitude correction term Delta lambda, and judgingAnd whether the delta lambda reaches a preset value, if so, executing the step S8, otherwise, returning to the step S2;

order toΔλ＝λ-λ_k. In conjunction with steps S4, S5, S6, the above equation may be written as:

if the receiver receives the propagation time difference Delta T between the two main stations and the secondary station_iThe following can be obtained:

the above formula can also be expressed as:

order toThe system of equations of the above equation can be expressed as:

solving by adopting a least square method to obtain:

further, the step S8 is specifically:

calculated receiver approximate position spherical coordinateI.e. spherical coordinates of the receiver

λ＝λ_k+1。

The invention has the beneficial effects that:

in the invention, in the process of probability iteration, the spherical longitude and latitude correction item of the probability is continuously calculated by taking the geodetic distance difference obtained by the known time difference as a reference to adjust the probability, so that the geodetic distance difference between the main station and the auxiliary station obtained by the probability is continuously close to the reference line, and when the spherical longitude and latitude correction item reaches a preset value, the calculated probability is the accurate position of the receiver. 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 of the modified macroellipse process of the present invention;

FIG. 3 is a schematic meridional ellipse of the present invention;

fig. 4 is a schematic view of the spherical triangle of the present invention.

Detailed Description

Example embodiments will now be described more fully with reference to the accompanying drawings. Example embodiments may, however, be embodied in many different 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 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 is described in detail with reference to the accompanying drawings, and the specific steps are as follows:

the principle of the invention is as follows: in the process of the approximate location iteration, the spherical longitude and latitude correction item of the approximate location is continuously calculated to adjust the approximate location 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 location is continuously close to the reference line, and when the spherical longitude and latitude correction item reaches a preset value, the approximate location obtained by calculation is the accurate position of the receiver.

The method is implemented according to the following steps:

step 1: geodetic coordinates (B) of input-related parameters and of receiver approximate bits_k,L_k)；

The input relevant parameters include: propagation time difference DeltaT between primary and secondary stations_iEarth coordinates (B) of a station participating in Rowland positioning solution_i,L_i) And geodetic coordinates (B) of the receiver profile_k,L_k). Subsequent calculations require conversion of the geodetic coordinates of the station and receiver profiles into corresponding spherical coordinates, which are expressed as stationSpherical coordinates of receiver profile are expressed as

Geodetic coordinates (B)_i,L_i) Coordinates with spherical surfaceThe conversion formula of (1) is as follows:

λ_i＝L_i

wherein f is the oblateness of the earth, a is the length of the longer half axis of the earth, b is the length of the shorter half axis of the earth, and f is 1-b/a.

Step 2: calculating the geodesic distance gamma from the updated probable bit to each station by adopting an improved nesting coefficient method or an improved large ellipse method_i；

First, the update approximate bit is calculated:

λ_k+1＝λ_k+Δλ

wherein the content of the first and second substances,represents the spherical latitude of the probable location after the (k + 1) th iteration,spherical latitude representing the k-th after-probable position, and, similarly, λ_k+1And λ_kAlso denoted spherical longitudes for the approximate after the k +1 th and k-th iterations respectively,representing a spherical latitude correction term, and delta lambda representing a spherical longitude correction term, and taking the correction term when iterating for the first timeΔλ＝0。

The invention introduces two algorithms for calculating the geodetic distance, namely an improved nesting coefficient method and an improved large ellipse method. The invention will specifically show the calculation step of the improved nesting coefficient method in step 2.1, and will specifically show the calculation step of the improved big ellipse method in step 2.2. When the geodetic distance is calculated, one algorithm is selected.

The specific steps of the improved nesting coefficient method are as follows:

step 2.1.1: and (3) calculating the spherical longitude difference delta lambda of the two points of the approximate position and the station on the auxiliary circular spherical surface by iteratively calculating the earth longitude difference correction item delta w. The geodetic coordinates of the approximate and station are respectively (B)_k,L_k) And (B)_i,L_i)。

The relationship between the spherical warp difference Δ λ and the earth warp difference Δ L is:

Δλ＝ΔL+Δw

ΔL＝L_i-L_k

for the first iteration, Δ w is taken to be 0. In the two formulas, if the values of the delta lambda and the delta L are not in [ -pi, pi ] after calculation, the values of the delta lambda and the delta L need to be subjected to addition and subtraction of 2 pi for normalization processing, so that the values of the delta lambda and the delta L are in a standard interval.

Step 2.1.2: calculating the spherical angular distance sigma from the probable position to the station_i：

Step 2.1.3: calculating the critical point of the large ground wire passing through the general position and the stationSpherical latitude of

Step 2.1.4: calculating the central angular distance between the earth wires of the general position and the station

Step 2.1.5: calculating an earth warp difference correction term Δ w:

K₃＝V[1+f+f²-V(3+7f-13V)]

and if the large-scale deviation correction item delta w does not reach the preset value, returning to execute the step 2.1.2, and if the large-scale deviation correction item delta w reaches the preset value, executing the step 2.1.6.

Step 2.1.6: calculating the large ground wire distance gamma between the two points of the approximate position and the station_i：

K₁＝1+t{1-t[3-t(5-11t)]/4}

K₂＝t{1-t[2-t(37-94t)/8]}

γ_i＝K₁b(σ_i-Δσ_i)

Wherein e is the first eccentricity of the earth, a is the length of the longer half axis of the earth,and completing the calculation of the improved nesting coefficient method.

The idea of the improved big ellipse method is as follows: the large ellipse arc length between two points is the large earth wire distance between two points. The method comprises the following specific steps:

as shown in the schematic diagram of the modified ellipsometry of FIG. 2, the approximate point is P_kThe station is P_i. The cross section ellipse and the orthogonal point of a certain meridian on the ellipsoid areThe geodetic latitude of the point is highest on the elliptical cross-section. P₀、Q₁、Q₂、Q_nIs the orthogonal point of the meridian and the equator on the ellipsoid. Phi is a_k、φ_i、Respectively represent P_k、P_i、Geocentric latitude of the point, A_k、A_iIndicating the spherical azimuth angle theta of the two points of the profile and the station_k、θ_iAnd the geocentric latitude of the station in the section ellipse is shown. A cross-sectional ellipse refers to the intersection of a plane made through two known points on the surface of an ellipsoid and the center of the ellipsoid with the ellipsoid.

Step 2.2.1: calculating geocentric latitude phi of the probable location and the station_k，φ_i：

Step 2.2.2: calculating the spherical azimuth angle A of the approximate station and the station_k，A_i：

Wherein Δ L ═ L_i-L_k. The quadrant determination principle of the spherical azimuth A is as follows (A refers to any spherical azimuth):

step 2.2.3: calculating geocentric latitude theta of the approximate position and the station in the section ellipse_k，θ_i：

Spherical relief surface right-angled triangleThe following can be obtained:

step 2.2.4: short semi-axis b for calculating section ellipse_s：

As shown in the meridian ellipse diagram of FIG. 3, ρ represents the distance from a point P on the ellipse to the center O, and φ represents the angle between the OP line and the x-axis. The standard equation for a meridional ellipse is:

wherein, a is the length of the longer half axis of the earth, and b is the length of the shorter half axis of the earth. The above formula is expressed by polar coordinates, and

substituting the formula into the formula and finishing to obtain:

as shown in the modified macrocephalic method of figure 2, at meridional ellipse NOQ_nIn the process, theGeocentric latitude of pointSubstituting the above formula to obtain a semi-minor axis with an oval cross section

Spherical relief surface right-angled triangleThe following can be obtained:

substituting the above formula intoIn the expression of (a), it can be:

step 2.2.5: calculating the cross-sectional ellipseFirst eccentricity e_s：

Step 2.2.6: calculating the large ellipse arc length gamma between two points_i；

The large ellipse arc length X on the cross-sectional ellipse can be expressed as:

X(θ)＝a(i₀θ+i₂ sin 2θ+i₄ sin 4θ+i₆ sin 6θ)

where θ represents the geocentric latitude of a point on the cross-sectional ellipse. The coefficients in the formula are:

can be obtained, P₀P_k、P₀P_iArc length of (d): gamma ray₁＝X(θ_k)、γ₂＝X(θ_i) Thus, the large elliptical arc length between two points, i.e., the geodetic distance between two points, is:

γ_i＝|γ₁-γ₂|＝|X(θ_k)-X(θ_i)|

and finishing the calculation of the improved big ellipse method.

And step 3: spherical azimuth angles A from the approximate position to each station are calculated and updated by adopting sine and cosine theorem of spherical triangle_iSine and cosine value sinA of_iAnd cosA_i；

As shown in the schematic diagram of spherical triangle in FIG. 4, the approximate point isThe station isIf the north pole is N, then the approximate bit point P_kSpherical distance to N point ofStation P_iSpherical distance to N point ofd representing probability point to stationSpherical distance. A. the_iRepresenting a probable point P_kTo station P_iΔ λ represents the spherical aberration between the isocenter and the station, where Δ λ ═ λ_i-λ_k。

Calculating and updating spherical azimuth A of probable position to station_iSine value sinA of_i：

Calculating and updating spherical azimuth A of probable position to station_iCosine value of (cosA)_i：

And 4, step 4: calculating and updating large earth wire distance difference xi between main station and secondary station_i；

In the Rowland positioning calculation, M represents the main station, S represents the auxiliary station, and then gamma (M) represents the geodesic distance from the probable bit to the main station, and gamma (S) represents the geodesic distance from the probable bit to the auxiliary station. Therefore, the large difference in ground distance between the primary and secondary stations calculated from the bit-profile is represented as:

ξ_i＝γ_i(S)-γ(M)

large difference in distance of ground lines between primary and secondary stations to be obtained from time differenceExpressed by Taylor's formula and in the approximate positionUnfolding:

and 5: calculating and updating spherical longitude and latitude correction termCorrelation coefficient U of sum Delta lambda_iAnd V_i；

Wherein A is_MIndicating the spherical azimuth of the overview to the main stage,indicating the spherical azimuth of the overview to the secondary stage.

Step 6: calculating and updating large earth wire distance difference xi between main station and secondary station obtained by time difference_t ⁱDifference ξ between the large distance and the main station and the secondary station derived from the position of the probability_iDifference Δ W of_i；

And 7: calculating and updating spherical latitude correction term by adopting least square algorithmAnd a spherical longitude correction term Delta lambda, and judgingAnd whether the delta lambda reaches a preset value, if so, executing the step 8, otherwise, returning to the step 2;

order toΔλ＝λ-λ_k. In connection with steps 4, 5, 6, the above equation can be written as:

if the receiver receives the propagation time difference Delta T between two main and auxiliary stations_iAnd finishing to obtain:

the above formula can also be expressed as:

order toThe system of equations of the above equation can be expressed as:

solving by adopting a least square method to obtain:

step 8: calculated receiver approximate position spherical coordinateI.e. spherical coordinates of the receiver

λ＝λ_k+1

And finishing the calculation.

The following tables are respectively under the conditions of a single chain and a double chain (hyperbolic positioning is divided into single chain positioning and multi-chain positioning, and multi-chain positioning is converted into double chain positioning after optimized station chains are carried out during multi-chain positioning, so that positioning errors solved by multi-chain positioning are positioning errors solved by double chain positioning), the geodesic distance algorithm is respectively positioning errors solved by a Rowland positioning method during an Andoyer-Lambert formula, an improved nesting coefficient method and an improved large ellipse method, and specific data are as follows:

single station chain positioning resolving error with apparent-big ground wire distance algorithm as Andoyer-Lambert formula

Single-chain positioning resolving error for improving nesting coefficient method by using table two geodesic distance algorithm

Single-chain positioning resolving error of table three-large ground wire distance algorithm for improving large ellipse method

Table four geodetic distance algorithm is double-station chain positioning resolving error of Andoyer-Lambert formula

Double-chain positioning resolving error of table five-earth wire distance algorithm for improving nesting coefficient method

Double-chain positioning resolving error of table six-large ground wire distance algorithm for improving large ellipse method

Examples

In the case of a single station chain, earth coordinates of the approximate position of the test points (22 degrees 40 ', 117 degrees 30') are input, the test station chain is an east China sea (8390) station chain, the main station of the test station chain is Xuan Cheng, the secondary stations of the test station chain are Rong and Luo, and the propagation time difference TD between the Xuan Cheng and the Rong is simulated₁2.4030269e-03 seconds, propagation time difference TD between Xuan Cheng and roughu₂For-2.6955252 e-03 seconds, the exact coordinates and positioning errors of the test points are calculated.

The method steps are implemented, and the ellipsoid parameters in the experiment are as follows: the earth's major semi-axis a is 6378140 meters, and the earth's minor semi-axis b is 6356755.288158 meters.

Under the same condition, the accurate coordinates and the positioning error data of the test points calculated by the method and the positioning calculation method for calculating the geodetic distance by using the conventional Andoyer-Lambert formula are shown in the seventh table. It can be seen from table seven that the positioning error of the test point calculated by the algorithm of the present invention is greatly reduced and the positioning accuracy is greatly improved compared with the positioning error of the test point calculated by the existing positioning calculation method for calculating the geodetic distance by using the Andoyer-Lambert formula. Thereby verifying the high-precision characteristic of the method.

Positioning error (unit: meter) of seven different positioning calculation algorithms of watch

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.