Method for improving TDOA positioning precision of MEOSAR system

Abstract

The invention discloses a method for improving TDOA positioning precision of an MEOSAR system, and relates to the technical field of satellite positioning. Comprising the following steps: analyzing a precise ephemeris, obtaining a satellite position and a satellite speed, and constructing a satellite position matrix and a time difference of arrival matrix of signals forwarded by a plurality of satellites; initializing a particle swarm under a space geodetic coordinate system, and setting position optimization parameters of particles and an initialization range corresponding to the parameters; evaluating the advantages and disadvantages of the particle positions according to the fitness function, and obtaining the optimal positions of particle individuals and the global optimal position; updating the positions of the particles in the constraint space, using boundary conditions to limit the updating positions of the particles to meet the solution space, and carrying out iteration; and updating parameters and judging whether the maximum number of iterations is reached, if not, continuing iteration, and if the algorithm is finished, outputting a globally optimal solution. Compared with a traditional PSO algorithm, the method has the advantages that the calculation precision is improved, the positioning error is remarkably reduced, and the positioning effect is better.

Description

Technical Field

[0001] The present invention relates to the field of satellite positioning technology, and more specifically to a method for improving TDOA positioning accuracy of a MEOSAR system. Background Art

[0002] As the main direction of the future development of the global satellite search and rescue system (COSPAS-SARSAT), the Medium Earth Orbit Search and Rescue Satellite System (MEOSAR) has become a key technology for responding to maritime, aviation and land rescue missions worldwide with its efficient, real-time and high-precision positioning capabilities. In the MEOSAR system, TDOA positioning technology is widely used in the positioning of distress signals, and can provide high-efficiency and high-rescue positioning services. However, the TDOA positioning algorithm itself is a typical nonlinear optimization problem, and it often faces great computational challenges when facing complex geographical environments and satellite coverage conditions. The traditional particle swarm optimization (PSO) algorithm is widely used in solving nonlinear optimization problems such as TDOA positioning due to its global search capability, strong adaptability and gradient independence. In the TDOA positioning problem of the MEOSAR system, the PSO algorithm can find the positioning solution through an iterative optimization process and achieve a good balance between computational efficiency and accuracy.

[0003] However, the PSO algorithm cannot be directly used for satellite passive positioning because it is a nonlinear optimization problem subject to the constraints of the earth's surface. The update mechanism of PSO depends on the position, velocity, individual optimal solution and global optimal solution of the particles. In the absence of constraints, the movement of particles depends entirely on these factors, which will generate solutions that do not meet the constraints, resulting in low positioning efficiency and large positioning errors. Therefore, PSO is not directly applicable to constrained optimization problems.

[0004] In recent years, with the continuous development and application of MEOSAR systems, how to improve positioning accuracy while ensuring computational efficiency has become an important research direction in the field of satellite search and rescue. To this end, many researchers have tried to further improve the accuracy and reliability of TDOA positioning by improving optimization algorithms and introducing constraints. In summary, improving the existing PSO algorithm so that it can adapt to the special constraints on the earth's surface has important significance and application value for improving the performance of the MEOSAR system.

[0005] Therefore, it is an urgent problem for those skilled in the art to propose a method for improving the TDOA positioning accuracy of the MEOSAR system to solve the difficulties existing in the prior art. Summary of the invention

[0006] In view of this, the present invention provides a method for improving the TDOA positioning accuracy of a MEOSAR system, which is used to solve the technical problems existing in the prior art.

[0007] In order to achieve the above object, the present invention provides the following technical solutions:

[0008] A method for improving TDOA positioning accuracy of a MEOSAR system, comprising the following steps:

[0009] Analyze precise ephemeris, obtain satellite position and satellite speed, and construct satellite position matrix and arrival time difference matrix of signals forwarded via multiple satellites;

[0010] Initialize the particle swarm in the spatial geodetic coordinate system, set the particle position optimization parameters and the initialization range corresponding to the parameters;

[0011] Evaluate the quality of particle positions according to the fitness function to obtain the best individual particle position and the global best position;

[0012] Update the particle position in the constraint space, use boundary conditions to limit the particle update position to satisfy the solution space, and iterate;

[0013] Update the parameters and determine whether the maximum number of iterations has been reached. If not, continue iterating. When the algorithm ends, output the global optimal solution.

[0014] Optionally, the specific content of obtaining the satellite position is:

[0015] The steps to obtain satellite positions using orbital parameters are as follows:

[0016] Calculate the average angular velocity of the satellite:

[0017]

[0018] Where n0 is the reference time t oe The average angular velocity of the earth is GM, which is the product of the gravitational constant G and the total mass of the earth M. GM = 3.986005 × 10 14 m 3 / s 2 ;

[0019] n=n0+Δn

[0020] Where Δn is the perturbation correction given by the broadcast ephemeris, and n is the average angular velocity of the satellite at the time of observation;

[0021] Calculate the satellite's mean anomaly when the signal is transmitted:

[0022] Δt=a0+a1(t 1 -t oc )+a2(t 1-t oc ) 2

[0023] t=t 1 -△t

[0024] Among them, t 1 is the observation time, t oc is the reference time, t is the time corrected by the satellite clock,

[0025] Δt is the satellite clock error, a0 is the satellite clock bias, a1 is the satellite clock drift, and a2 is the satellite clock drift rate;

[0026] M k =M0+n(tt oc )

[0027] Among them, M0 is the mean anomaly angle at the reference time, M k is the mean anomaly angle of the satellite when the signal is transmitted;

[0028] Calculate the eccentric anomaly angle:

[0029] E k =M k +esinE k

[0030] Calculate the true anomaly angle:

[0031]

[0032] Where, e is the eccentricity of the satellite orbit;

[0033] Calculate the ascending angle:

[0034]

[0035] Where ω is the periapsis angular distance;

[0036] Compute the perturbation correction:

[0037]

[0038]

[0039] Among them, C uc , C us , C rc , C rs , C ic , C is are 6 perturbation correction parameters, δ uk is the perturbation correction term of the ascending distance angle, δ rk is the satellite radial vector perturbation correction term, δ ikis the satellite orbit inclination perturbation correction term;

[0040] Ascending intersection angle, satellite radius vector, and orbital inclination after perturbation correction:

[0041]

[0042] Where i0 is t oe The orbital inclination at the moment is given by the Kepler six parameters of the broadcast ephemeris; is the rate of change of orbital inclination i, given by the nine perturbation parameters in the broadcast ephemeris; A is the major radius of the satellite orbit; is the satellite radial vector perturbation correction term;

[0043] Calculate the coordinates of the satellite in the orbital plane coordinate system:

[0044] x′ k =r k cosu k

[0045] y′ k =r k Sinu k

[0046] Among them, x′ k is the projection coordinate of the satellite on the orbital plane along the horizontal axis of the orbital plane, y′ k is the projection coordinate of the satellite on the orbital plane along the longitudinal axis of the orbital plane;

[0047] Calculate the right ascension of the ascending node Ω:

[0048]

[0049] in, is the rate of change of the ascending node with respect to time, Ω0 is the right ascension of the ascending node at the reference time, Ω0 and Provided by satellite ephemeris, Ω k To calculate the time t k The right ascension of the ascending node, is the angular velocity of the Earth's rotation, t k is the calculation time of the satellite position, t oe is the reference time of the satellite ephemeris;

[0050] Calculate the coordinates of the satellite in the Earth-fixed coordinate system:

[0051]

[0052] Optionally, the specific content of obtaining satellite speed is:

[0053] The steps to obtain satellite speed using orbital parameters are as follows:

[0054] Calculate the first-order derivatives of the mean anomaly angle and the eccentric anomaly angle with respect to time respectively:

[0055]

[0056] in, is the first-order derivative of the mean anomaly with respect to time, is the first-order derivative of the true anomaly with respect to time;

[0057] The first-order derivative of the true anomaly is:

[0058]

[0059] The first derivative of the angular distance of the ascending node:

[0060]

[0061] in, It is the first-order derivative of the true anomaly with respect to time, indicating the rate of change of the true position angle of the satellite in orbit with time. is the first derivative of the angular distance of the ascending node with respect to time;

[0062] The first derivative of the perturbation correction term:

[0063]

[0064] in, is the first-order derivative of the perturbation correction term of the ascending distance angle with respect to time, is the first-order derivative of the satellite radial vector perturbation correction term with respect to time, is the first-order derivative of the satellite orbit inclination correction term with respect to time;

[0065] The first-order derivatives of the ascending node angular distance, satellite radius vector, and orbital inclination after perturbation correction are:

[0066]

[0067] in, is the first-order derivative of the angular distance of the ascending node with respect to time, indicating the rate of change of the angular distance of the ascending node with time. is the first-order derivative of the satellite radius vector, which indicates the rate of change of the distance from the satellite to the center of the earth over time. is the first-order derivative of the orbital inclination with respect to time, indicating the rate of change of the orbital inclination with time, is the orbital inclination, that is, the angle between the orbital plane and the Earth's equatorial plane, is the first derivative with respect to time, indicating the rate at which the orbital inclination changes with time;

[0068] Calculate the satellite's velocity in the orbital plane:

[0069]

[0070] in, is the first-order derivative of the satellite's projection coordinates on the orbital plane along the horizontal axis of the orbital plane with respect to time, that is, the velocity component of the satellite along the horizontal axis of the orbital plane. It is the first-order derivative of the projection coordinates of the satellite on the orbital plane along the longitudinal axis of the orbital plane with respect to time, that is, the velocity component of the satellite along the longitudinal axis of the orbital plane;

[0071] The first derivative of the right ascension of the ascending node:

[0072]

[0073] in, is the first derivative of the right ascension of the ascending node, is the angular velocity of the ascending node right ascension, which refers to the precession velocity of the satellite orbit. is the rate of change of right ascension caused by the rotation of the Earth;

[0074] Calculate the satellite's velocity in the ECEF coordinate system:

[0075]

[0076] Optionally, the particle swarm is initialized in the spatial geodetic coordinate system, and the specific contents of setting the particle position optimization parameters and the initialization range corresponding to the parameters are as follows:

[0077] Set the particle position optimization parameters to B, L, H, which are latitude, longitude, and altitude respectively;

[0078] x(i, 1) = 180*rand()-90; % Initialize the latitude in the range of [-90, 90];

[0079] x(i, 2) = 360*rand()-180; % Initialize longitude in the range of [-180, 180];

[0080] x(i, 3) = 0.2*rand()-0.1; % Initialization height is within the range of [-100, 100] meters;

[0081] Among them, i represents the i-th particle, x(i, j) represents the j-th position optimization parameter of the i-th particle, x(i, 1) represents latitude B, x(i, 2) represents longitude L, and x(i, 3) represents altitude H.

[0082] Optionally, the fitness function is:

[0083]

[0084] Among them, s iis the position coordinate of the ith satellite, r is the position coordinate of the ground station, g(y) is the position coordinate of the distress beacon, s j is the position coordinate of the jth satellite, c is the speed of light, TDOA ij is the time difference between the distress signal passing through the ith satellite and the ith satellite and arriving at the ground station.

[0085] Optionally, update the particle position in the constraint space, use boundary conditions to limit the particle update position to satisfy the solution space, and perform the iteration as follows:

[0086] The particle velocity is updated as follows:

[0087] v i (t+1)=wv i (t)+c1r1(pBest i -x i (t))+c2r2(gBest-x i (t))

[0088] Among them, v i (t) is the velocity vector of particle i at the tth iteration, w is the inertia weight, which controls the inertia of the particle, that is, the influence of the velocity of the previous step on the current velocity, c1 is the cognitive coefficient, that is, the individual learning factor, which controls the degree to which the particle approaches its historical optimal position, c2 is the social coefficient, that is, the group learning factor, which controls the degree to which the particle approaches the global optimal position, r1 and r2 are random numbers between [0, 1], which increase the randomness and diversity of the search, and pBest i is the best historical position of particle i in the tth iteration, gBest is the global best historical position, which means the best position of all particles in the group in the tth iteration;

[0089] The particle position update formula is:

[0090] x i (t+1)=x i (t)+v i (t+1)

[0091] Among them, x i (t) is the position vector of particle i in the tth iteration, v i (t+1) is the velocity vector of particle i in the t+1th iteration.

[0092] It can be seen from the above technical solutions that, compared with the prior art, the present invention discloses a method for improving the TDOA positioning accuracy of a MEOSAR system, and its beneficial effects are:

[0093] 1) Compared with the traditional PSO algorithm, the calculation accuracy is improved, the positioning error is significantly reduced, and the positioning effect is better;

[0094] 2) The search space of particles is constrained within the range of restrictions, which greatly reduces the ineffective exploration of the particle swarm, saves resources and enables particles to find the optimal solution faster, greatly speeding up the convergence of the algorithm and achieving higher efficiency in practical applications;

[0095] 3) It has stronger anti-noise ability and adaptability. It can significantly improve positioning accuracy when the signal-to-noise ratio is low, and continuously optimize the results when the signal quality gradually improves. It shows stronger robustness and stability than the standard PSO algorithm, especially in weak signal environments, showing obvious advantages. BRIEF DESCRIPTION OF THE DRAWINGS

[0096] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying creative work.

[0097] Figure 1 A flow chart of a method for improving TDOA positioning accuracy of a MEOSAR system provided by the present invention;

[0098] Figure 2 A schematic diagram of the particle confinement space provided by the present invention;

[0099] Figure 3 Schematic diagram of positioning error under medium signal-to-noise ratio conditions provided by the present invention; wherein 3a is an optimization algorithm and 3b is a traditional algorithm;

[0100] Figure 4 Schematic diagram of the convergence curve provided by the present invention; wherein 4a is the optimization algorithm, and 4b is the traditional algorithm;

[0101] Figure 5 This is a comparison chart of positioning errors under different signal-to-noise ratio conditions provided by the present invention. DETAILED DESCRIPTION

[0102] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0103] See also Figure 1As shown, the present invention discloses a method for improving the TDOA positioning accuracy of a MEOSAR system, comprising the following steps:

[0104] Analyze precise ephemeris, obtain satellite position and satellite speed, and construct satellite position matrix and arrival time difference matrix of signals forwarded via multiple satellites;

[0105] Initialize the particle swarm in the spatial geodetic coordinate system, set the particle position optimization parameters and the initialization range corresponding to the parameters;

[0106] Evaluate the quality of particle positions according to the fitness function to obtain the best individual particle position and the global best position;

[0107] Update the particle position in the constraint space, use boundary conditions to limit the particle update position to satisfy the solution space, and iterate;

[0108] Update the parameters and determine whether the maximum number of iterations has been reached. If not, continue iterating. When the algorithm ends, output the global optimal solution.

[0109] Specifically, the input of this method is the objective function to be optimized and various parameters including the number of particles, learning factor, maximum number of iterations, optimization target dimension, satellite position matrix, ground station position, and arrival time difference matrix. The output is the global optimal position (B, L, H) of the particle swarm, that is, the distress beacon solution position in the spatial geodetic coordinate system.

[0110] In the MEOSAR system, the signals sent by the distress beacon pass through different medium-orbit satellites to reach the ground receiving station. The ground station antenna receives these signals, records and compares the time of signal arrival, and then obtains the time difference, which reflects the difference in distance between the signal from the beacon to different satellites. A set of time differences defines a hyperbolic surface, and the location of the signal source is located on this hyperbolic surface. Multiple sets of time differences can define multiple hyperbolic surface equations. Therefore, more than three time difference equations can uniquely determine the location of the signal source.

[0111] Considering the equation described in the Earth-centered Earth-fixed coordinate system, the observation equation for the arrival time of the signal forwarded by satellite i can be described as:

[0112]

[0113] Among them, (x i ,y i , z i ) is the satellite S i The location coordinates of the beacon signal, (x, y, z) are the location coordinates of the beacon signal, (x g ,y g , z g) is the position coordinate of the ground receiving station, T0 is the transmission time of the beacon signal, which is an unknown quantity, Δt is the advance of the beacon clock relative to the ground station clock, which is an unknown quantity, δL i are various error parameters, including ionospheric error, tropospheric error, multipath effect error, etc.

[0114] Subtracting the arrival time observation equations of two different satellites can obtain a set of arrival time difference equations, which can eliminate the unknowns Δt and T0, and offset a large part of the time error and the error caused by the multipath effect:

[0115]

[0116] TDOA ij =TOA i -TOA j

[0117] Combine different equations for time difference and earth ellipse: The location of the distress beacon can be determined.

[0118] In the MEOSAR problem, we want to solve the beacon position, so the construction idea of ​​the fitness function should be centered around how to minimize the positioning error, and is set as:

[0119]

[0120] Among them, S i , S j are the coordinates of the i-th and j-th satellites respectively, r is the position coordinate of the ground station, and p is the position coordinate of the beacon to be solved. Since the distance is calculated, the above coordinates are all in the Cartesian coordinate system.

[0121] The earth is an approximate ellipsoid, and limiting the latitude, longitude, and altitude within a specific range is intuitive in the geodetic coordinate system. There is a height difference between the earth's ellipsoid reference model and the real earth's surface due to terrain undulations such as mountains, oceans, and other terrain features. Excluding extreme terrain such as high mountains and deep valleys, the average height difference between most areas of the world's land and ocean and the ellipsoid surface is usually around plus or minus 100 meters, so the range of latitude, longitude, and altitude is set to:

[0122] -90 ° ≤B≤90 °

[0123] -180 ° ≤L≤180 °

[0124] -100≤H≤100

[0125] The constraint condition that the particle swarm flies on or near the surface of the earth can be met.

[0126] When evaluating the particle position, it is necessary to convert the spatial geodetic coordinate system into the Cartesian coordinate system and then substitute it into the fitness function:

[0127]

[0128] According to the fitness value, the individual optimal position and the global optimal position under the current number of iterations are determined and the update continues. The particle velocity update formula and position update formula are as follows:

[0129] v i (t+1)=wv i (t)+c1r1(pBest i -x i (t))+c2r2(gBest-x i (t))

[0130] x i (t+1)=x i (t)+v i (t+1)

[0131] Among them, the velocity of particle i at generation t is v i (t), position x i (t), the individual optimal position is pBest i , the global optimal position is gBest.

[0132] w is called the inertia weight, which controls the particle's ability to maintain its current velocity as it moves forward (a larger inertia weight allows the particle to search extensively in the search space, while a smaller inertia weight allows the particle to search locally).

[0133] The following dynamic adjustment method settings are used in the iteration, which gradually decreases as the number of iterations increases:

[0134]

[0135] Among them, t is the current iteration number, M is the maximum iteration number, ω max ,ω min These are the maximum and minimum inertia weights set, generally 0.9 and 0.4.

[0136] c1 and c2 are learning factors that determine the speed at which particles approach the individual optimal position and the group optimal position, and are generally set to the same constant 2.

[0137] r1 and r2 are random numbers between [0, 1], which increase the randomness of the search.

[0138] The following strategy is used during the update process to limit the particles to update within the constraints, ensuring that the algorithm can use reasonable latitude, longitude, and altitude values ​​throughout the optimization process:

[0139] B = mod (B + 90, 180) - 90

[0140] L = mod (L + 180, 360) - 180

[0141] H=0.2*rand()-0.1

[0142] When the algorithm ends, the currently found global best position and its corresponding fitness value are output as the optimal solution to the problem.

[0143] For details, see Figure 2 The figure shows a schematic diagram of particle confinement space.

[0144] Furthermore, the specific content of obtaining the satellite position is:

[0145] The steps to obtain satellite positions using orbital parameters are as follows:

[0146] Calculate the average angular velocity of the satellite:

[0147]

[0148] Where n0 is the reference time t oe The average angular velocity of the earth is GM, which is the product of the gravitational constant G and the total mass of the earth M. GM = 3.986005 × 10 14 m 3 / s 2 ;

[0149] n=n0+Δn

[0150] Where Δn is the perturbation correction given by the broadcast ephemeris, and n is the average angular velocity of the satellite at the time of observation;

[0151] Calculate the satellite's mean anomaly when the signal is transmitted:

[0152] Δt=a0+a1(t 1 -t oc )+a2(t 1 -t oc ) 2

[0153] t=t 1 -△t

[0154] Among them, t 1 is the observation time, t ocis the reference time, t is the time after the satellite clock correction, Δt is the satellite clock error, which indicates the total deviation of the satellite clock relative to the standard time (such as GPS system time or other reference time), a0 is the satellite clock bias, a1 is the satellite clock drift (indicating the rate of change of the satellite in unit time), and a2 is the satellite clock drift rate (indicating the acceleration of the clock drift over time, usually a second-order correction term used to compensate for nonlinear drift);

[0155] M k =M0+n(tt oc )

[0156] Among them, M0 is the mean anomaly angle at the reference time, M k is the mean anomaly angle of the satellite when the signal is transmitted;

[0157] Calculate the eccentric anomaly angle:

[0158] E k =M k +esinE k

[0159] Calculate the true anomaly angle:

[0160]

[0161] Where, e is the eccentricity of the satellite orbit;

[0162] Calculate the ascending angle:

[0163]

[0164] Where ω is the periapsis angular distance;

[0165] Compute the perturbation correction:

[0166]

[0167] Among them, C uc , C us , C rc , C rs , C ic , C is are 6 perturbation correction parameters, δ uk is the perturbation correction term of the ascending distance angle, δ rk is the satellite radial vector perturbation correction term, δ ik is the satellite orbit inclination perturbation correction term;

[0168] Ascending intersection angle, satellite radius vector, and orbital inclination after perturbation correction:

[0169]

[0170]

[0171]

[0172] Where i0 is t oe The orbital inclination at the moment is given by the Kepler six parameters of the broadcast ephemeris; is the rate of change of orbital inclination i, given by the nine perturbation parameters in the broadcast ephemeris; A is the major radius of the satellite orbit; is the satellite radial vector perturbation correction term;

[0173] Calculate the coordinates of the satellite in the orbital plane coordinate system:

[0174] x′ k =r k cosu k

[0175] y′ k =r k Sinu k

[0176] Among them, x′ k is the projection coordinate of the satellite on the orbital plane along the horizontal axis of the orbital plane, y′ k is the projection coordinate of the satellite on the orbital plane along the longitudinal axis of the orbital plane;

[0177] Specifically, the orbital plane refers to the orbital plane of the satellite, which is uniquely determined by the satellite's orbital elements (six Kepler parameters). The two coordinates describe the two-dimensional projection position of the satellite in the orbital plane. Its physical meaning is the components of the satellite in two directions along the orbital plane.

[0178] Calculate the right ascension of the ascending node Ω:

[0179]

[0180] in, is the rate of change of the ascending node with respect to time, Ω0 is the right ascension of the ascending node at the reference time, Ω0 and Provided by satellite ephemeris, Ω k To calculate the time t k The right ascension of the ascending node, is the angular velocity of the Earth's rotation (about 7.2921151467×10- 5 rad / s), t k is the calculation time of the satellite position, t oe is the reference time of the satellite ephemeris;

[0181] Calculate the coordinates of the satellite in the Earth-fixed coordinate system:

[0182]

[0183] Furthermore, the specific content of obtaining the satellite speed is:

[0184] The steps to obtain satellite speed using orbital parameters are as follows:

[0185] Calculate the first-order derivatives of the mean anomaly angle and the eccentric anomaly angle with respect to time respectively:

[0186]

[0187] in, is the first-order derivative of the mean anomaly with respect to time, is the first-order derivative of the true anomaly with respect to time;

[0188] The first-order derivative of the true anomaly is:

[0189]

[0190] The first derivative of the angular distance of the ascending node:

[0191]

[0192] in, It is the first-order derivative of the true anomaly with respect to time, indicating the rate of change of the true position angle of the satellite in orbit with time. is the first derivative of the angular distance of the ascending node with respect to time;

[0193] The first derivative of the perturbation correction term:

[0194]

[0195] in, is the first-order derivative of the perturbation correction term of the ascending distance angle with respect to time, is the first-order derivative of the satellite radial vector perturbation correction term with respect to time, is the first-order derivative of the satellite orbit inclination correction term with respect to time;

[0196] The first-order derivatives of the ascending node angular distance, satellite radius vector, and orbital inclination after perturbation correction are:

[0197]

[0198] in, is the first-order derivative of the angular distance of the ascending node with respect to time, indicating the rate of change of the angular distance of the ascending node with time. is the first-order derivative of the satellite radius vector, which indicates the rate of change of the distance from the satellite to the center of the earth over time. is the first-order derivative of the orbital inclination with respect to time, indicating the rate of change of the orbital inclination with time, is the orbital inclination, the angle between the orbital plane and the Earth's equatorial plane, is the first derivative with respect to time, indicating the rate at which the orbital inclination changes with time;

[0199] Calculate the satellite's velocity in the orbital plane:

[0200]

[0201] in, is the first-order derivative of the satellite's projection coordinates on the orbital plane along the horizontal axis of the orbital plane with respect to time, that is, the velocity component of the satellite along the horizontal axis of the orbital plane. It is the first-order derivative of the projection coordinates of the satellite on the orbital plane along the longitudinal axis of the orbital plane with respect to time, that is, the velocity component of the satellite along the longitudinal axis of the orbital plane;

[0202] The first derivative of the right ascension of the ascending node:

[0203]

[0204] in, is the first derivative of the right ascension of the ascending node, is the angular velocity of the ascending node right ascension, which refers to the precession velocity of the satellite orbit. is the rate of change of right ascension caused by the rotation of the Earth;

[0205] Calculate the satellite's velocity in the ECEF coordinate system:

[0206]

[0207] Furthermore, the particle swarm is initialized in the spatial geodetic coordinate system, and the specific contents of setting the particle position optimization parameters and the initialization range corresponding to the parameters are as follows:

[0208] Set the particle position optimization parameters to B, L, H, which are latitude, longitude, and altitude respectively;

[0209] x(i, 1) = 180*rand()-90; % Initialize the latitude in the range of [-90, 90];

[0210] x(i, 2) = 360*rand()-180; % Initialize longitude in the range of [-180, 180];

[0211] x(i, 3) = 0.2*rand()-0.1; % Initialization height is within the range of [-100, 100] meters;

[0212] Among them, i represents the i-th particle, x(i, j) represents the j-th position optimization parameter of the i-th particle, x(i, 1) represents latitude B, x(i, 2) represents longitude L, and x(i, 3) represents altitude H.

[0213] Furthermore, the fitness function is:

[0214]

[0215] Among them, s i is the position coordinate of the ith satellite, r is the position coordinate of the ground station, g(y) is the position coordinate of the distress beacon, s j is the position coordinate of the jth satellite, c is the speed of light, TDOA ij is the time difference between the distress signal passing through the ith satellite and the ith satellite and arriving at the ground station.

[0216] Specifically, the above-mentioned position coordinates are all coordinates in the Cartesian coordinate system, and f=g(y) is a conversion function that converts spatial geodetic coordinates into Cartesian coordinates. The formula is as follows:

[0217] The function f = g(y) is defined by the following coordinate transformation formula:

[0218]

[0219] in, is the radius of curvature of the geocentric circle; f = (x, y, z) is the geocentric coordinate in the Cartesian coordinate system; y = (B, L, H) is the three-dimensional coordinate in the geodetic ellipsoid coordinate system; e is the ellipsoid ellipsoidality.

[0220] Specifically, the designed fitness function covers 4 satellites and uses TDOA to minimize the error between every two satellites.

[0221] Furthermore, the particle position is updated in the constraint space, and the boundary conditions are used to limit the particle update position to satisfy the solution space. The specific content of the iteration is:

[0222] The particle velocity is updated as follows:

[0223] v i (t+1)=wv i (t)+c1r1(pBest i -x i (t))+c2r2(gBest-x i (t))

[0224] Among them, v i (t) is the velocity vector of particle i at the tth iteration, w is the inertia weight, which controls the inertia of the particle, that is, the influence of the velocity of the previous step on the current velocity, c1 is the cognitive coefficient, that is, the individual learning factor, which controls the degree to which the particle approaches its historical optimal position, c2 is the social coefficient, that is, the group learning factor, which controls the degree to which the particle approaches the global optimal position, r1 and r2 are random numbers between [0, 1], which increase the randomness and diversity of the search, and pBest iis the best historical position of particle i in the tth iteration, gBest is the global best historical position, which means the best position of all particles in the group in the tth iteration;

[0225] The particle position update formula is:

[0226] x i (t+1)=x i (t)+v i (t+1)

[0227] Among them, x i (t) is the position vector of particle i in the tth iteration, v i (t+1) is the velocity vector of particle i in the t+1th iteration.

[0228] Specifically, after each update is completed, it is checked whether the three update parameters of the particle are within the constraint range, and the longitude and latitude are checked to be limited within the constraint range by means of modulo operation, so as to realize cyclic mapping without complex logical judgment and ensure the continuity of parameters. No matter how large or small the input value is, it will be smoothly mapped to the corresponding restriction range without jumping or breakpoints, which is helpful to optimize the smooth operation of the algorithm.

[0229] Moreover, longitude and latitude are essentially periodic (characteristics of the earth's spherical coordinates). This method can make good use of this characteristic and cycle out-of-range values ​​back to the legal interval instead of simply truncating them, which is more in line with the actual meaning of geographic coordinates.

[0230] if x(i, 1)>90||x(i, 1)<-90

[0231] x(i,1)=mod(x(i,1)+90,180)-90;

[0232] end

[0233] if x(i, 2)>180||x(i, 2)<-180

[0234] x(i,2)=mod(x(i,2)+180,360)-180;

[0235] end

[0236] The height is adjusted directly using the initialization method, which is simple, intuitive and effective:

[0237] if x(i,3)>0.1||x(i,3)<-0.1

[0238] x(i,3)=0.2*rand()-0.1;

[0239] end

[0240] Specifically, during simulation, the fixed error is set to 3 microseconds for solving the ephemeris data, and the TDOA measurement error is set to 0.3 microseconds to meet the actual measurement conditions. The total number of TDOA algorithm simulations is 100 times, and the number of particle swarm algorithm iterations is set to 300 times and the number of particles is 150 for each TDOA run. Different signal-to-noise ratios will affect the TDOA measurement value and thus affect the positioning error. Therefore, see Figure 5 As shown in the figure, the positioning error under different signal-to-noise ratio conditions is simulated. However, the rules reflected by different signal-to-noise ratios are the same, see Figure 3 As shown in Figure 3a and 3b, the simulation results under the condition of 10dB signal-to-noise ratio are more in line with the actual signal-to-noise ratio conditions. The traditional algorithm does not set longitude, latitude, and altitude restrictions. Figure 4 As shown, 4a and 4b show the convergence curves of the two algorithms under 10dB signal-to-noise ratio conditions.

[0241] The various embodiments in this specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be referenced to each other.

[0242] The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for improving TDOA positioning accuracy of a MEOSAR system, characterized in that: The following steps are involved: Analyze precise ephemeris, obtain satellite position and satellite speed, and construct satellite position matrix and arrival time difference matrix of signals forwarded via multiple satellites; Initialize the particle swarm in the spatial geodetic coordinate system, set the particle position optimization parameters and the initialization range corresponding to the parameters; Evaluate the quality of particle positions according to the fitness function to obtain the best individual particle position and the global best position; Update the particle position in the constraint space, use boundary conditions to limit the particle update position to satisfy the solution space, and iterate; Update the parameters and determine whether the maximum number of iterations has been reached. If not, continue iterating. When the algorithm ends, output the global optimal solution.

2. A method for improving TDOA positioning accuracy of a MEOSAR system according to claim 1, characterized in that: The specific contents of obtaining satellite positions are: The steps to obtain satellite positions using orbital parameters are as follows: Calculate the average angular velocity of the satellite: Where n0 is the reference time t oe The average angular velocity of the earth is GM, which is the product of the gravitational constant G and the total mass of the earth M. GM = 3.986005 × 10 14 m 3 / s 2 ; n=n0+Δn Where Δn is the perturbation correction given by the broadcast ephemeris, and n is the average angular velocity of the satellite at the time of observation; Calculate the satellite's mean anomaly when the signal is transmitted: Δt=a0+a1(t 1 -t oc )+a2(t 1 -t oc ) 2 t=t 1 -Δt Among them, t 1 is the observation time, t oc is the reference time, t is the time after satellite clock correction, Δt is the satellite clock error, a0 is the satellite clock bias, a1 is the satellite clock drift, and a2 is the satellite clock drift rate; M k =M0+n(t-t oc ) Among them, M0 is the mean anomaly angle at the reference time, M k is the mean anomaly of the satellite when the signal is transmitted; Calculate the eccentric anastomosis angle: Yes k =M k +EsinE k Calculate the true anomaly angle: Where, e is the eccentricity of the satellite orbit; Calculate the ascending angle: Where ω is the periapsis angular distance; Compute the perturbation correction: Among them, C uc , C us , C rc , C rs , C ic , C is are 6 perturbation correction parameters, δ uk is the perturbation correction term of the ascending distance angle, δ rk is the satellite radial vector perturbation correction term, δ ik is the satellite orbit inclination perturbation correction term; Ascending intersection angle, satellite radius vector, and orbital inclination after perturbation correction: Where i0 is t oe The orbital inclination at the moment is given by the Kepler six parameters of the broadcast ephemeris; is the rate of change of orbital inclination i, given by the nine perturbation parameters in the broadcast ephemeris; A is the major radius of the satellite orbit; is the satellite radial vector perturbation correction term; Calculate the coordinates of the satellite in the orbital plane coordinate system: x' k =r k something k y' k =r k sin k Among them, x' k is the projection coordinate of the satellite on the orbital plane along the horizontal axis of the orbital plane, y' k is the projection coordinate of the satellite on the orbital plane along the longitudinal axis of the orbital plane; Calculate the right ascension of the ascending node Ω: in, is the rate of change of the ascending node with respect to time, Ω0 is the right ascension of the ascending node at the reference time, Ω0 and Provided by satellite ephemeris, Ω k To calculate the time t k The right ascension of the ascending node, is the angular velocity of the Earth's rotation, t k is the calculation time of the satellite position, t oe is the reference time of the satellite ephemeris; Calculate the coordinates of the satellite in the Earth-fixed coordinate system:

3. A method for improving TDOA positioning accuracy of a MEOSAR system according to claim 1, characterized in that: The specific content of obtaining satellite speed is: The steps to obtain satellite speed using orbital parameters are as follows: Calculate the first-order derivatives of the mean anomaly angle and the eccentric anomaly angle with respect to time respectively: in, is the first-order derivative of the mean anomaly with respect to time, is the first-order derivative of the true anomaly with respect to time; The first-order derivative of the true anomaly is: The first derivative of the angular distance of the ascending node: in, It is the first-order derivative of the true anomaly with respect to time, indicating the rate of change of the true position angle of the satellite in orbit with time. is the first derivative of the angular distance of the ascending node with respect to time; The first derivative of the perturbation correction term: in, is the first-order derivative of the perturbation correction term of the ascending distance angle with respect to time, is the first-order derivative of the satellite radial vector perturbation correction term with respect to time, is the first-order derivative of the satellite orbit inclination correction term with respect to time; The first-order derivatives of the ascending node angular distance, satellite radius vector, and orbital inclination after perturbation correction are: in, is the first-order derivative of the angular distance of the ascending node with respect to time, indicating the rate of change of the angular distance of the ascending node with time. is the first-order derivative of the satellite radius vector, which indicates the rate of change of the distance from the satellite to the center of the earth over time. is the first-order derivative of the orbital inclination with respect to time, indicating the rate of change of the orbital inclination with time, is the orbital inclination, that is, the angle between the orbital plane and the Earth's equatorial plane, is the first-order derivative with respect to time, indicating the rate at which the orbital inclination changes with time; Calculate the satellite's velocity in the orbital plane: in, is the first-order derivative of the satellite's projection coordinates on the orbital plane along the horizontal axis of the orbital plane with respect to time, that is, the velocity component of the satellite along the horizontal axis of the orbital plane. It is the first-order derivative of the projection coordinates of the satellite on the orbital plane along the longitudinal axis of the orbital plane with respect to time, that is, the velocity component of the satellite along the longitudinal axis of the orbital plane; The first derivative of the right ascension of the ascending node: in, is the first derivative of the right ascension of the ascending node, is the angular velocity of the ascending node right ascension, which refers to the precession velocity of the satellite orbit. is the rate of change of right ascension caused by the rotation of the Earth; Calculate the satellite's velocity in the ECEF coordinate system:

4. A method for improving TDOA positioning accuracy of a MEOSAR system according to claim 1, characterized in that: Initialize the particle swarm in the space geodetic coordinate system, set the particle position optimization parameters and the initialization range corresponding to the parameters as follows: Set the particle position optimization parameters to B, L, H, which are latitude, longitude, and altitude respectively; x(i,1)=180*rand()-90; % Initialize the latitude in the range of [-90,90]; x(i,2)=360*rand()-180; % Initialize longitude in the range of [-180,180]; x(i,3)=0.2*rand()-0.1;%Initialize the height within the range of [-100,100] meters; Among them, i represents the i-th particle, x(i,j) represents the j-th position optimization parameter of the i-th particle, x(i,1) represents latitude B, x(i,2) represents longitude L, and x(i,3) represents altitude H.

5. A method for improving TDOA positioning accuracy of a MEOSAR system according to claim 1, characterized in that: The fitness function is: Among them, s i is the position coordinate of the ith satellite, r is the position coordinate of the ground station, g(y) is the position coordinate of the distress beacon, s j is the position coordinate of the jth satellite, c is the speed of light, TDOA ij is the time difference between the distress signal passing through the i-th satellite and the j-th satellite and arriving at the ground station.

6. A method for improving TDOA positioning accuracy of a MEOSAR system according to claim 1, characterized in that: Update the particle position in the constraint space, use boundary conditions to limit the particle update position to satisfy the solution space, and the specific content of the iteration is: The particle velocity is updated as follows: v i (t+1)=wv i (t)+c1r1(pBest i -x i (t))+c2r2(gBest-x i (t)) Among them, v i (t) is the velocity vector of particle i at the tth iteration, w is the inertia weight, which controls the inertia of the particle, that is, the influence of the velocity of the previous step on the current velocity, c1 is the cognitive coefficient, that is, the individual learning factor, which controls the degree to which the particle approaches its historical optimal position, c2 is the social coefficient, that is, the group learning factor, which controls the degree to which the particle approaches the global optimal position, r1 and r2 are random numbers between [0,1], and pBest i is the best historical position of particle i in the tth iteration, gBest is the global best historical position, which means the best position of all particles in the group in the tth iteration; The particle position update formula is: x i (t+1)=x i (t)+v i (t+1) Among them, x i (t) is the position vector of particle i in the tth iteration, v i (t+1) is the velocity vector of particle i in the t+1th iteration.

Images