A method for improving positioning accuracy of a MEOSAR system TDOA

Abstract

This invention discloses a method for improving the TDOA positioning accuracy of the MEOSAR system, relating to the field of satellite positioning technology. The method includes: analyzing precise ephemeris data to obtain satellite positions and velocities, constructing a satellite position matrix and a time difference of arrival (TDOA) matrix for signals relayed by multiple satellites; initializing a particle swarm in a geodetic coordinate system, setting particle position optimization parameters and their corresponding initialization ranges; evaluating the quality of particle positions based on a fitness function to obtain the optimal individual particle position and the global optimal position; updating particle positions within a constrained space, using boundary conditions to restrict particle position updates to satisfy the solution space, and iterating; updating parameters and determining if the maximum number of iterations has been reached; if not, continuing iteration; and finally, outputting the global optimal solution. Compared to the traditional PSO algorithm, this invention improves computational accuracy, significantly reduces positioning errors, and achieves better positioning results.

Description

Technical Field

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

[0002] The Medium Earth Orbit Search and Rescue Satellite System (MEOSAR), as a major future development direction of the COSPAS-SARSAT global satellite search and rescue system, has become a key technology for global maritime, air, and land rescue missions due to its efficient, real-time, and high-precision positioning capabilities. Within the MEOSAR system, TDOOA (Trace-to-Operation) positioning technology is widely used for locating distress signals, providing highly efficient and accurate positioning services. However, the TDOA positioning algorithm itself is a typical nonlinear optimization problem, often facing significant computational challenges in complex geographical environments and satellite coverage conditions. The traditional Particle Swarm Optimization (PSO) algorithm, with its global search capability, strong adaptability, and gradient independence, is widely used to solve nonlinear optimization problems such as TDOA positioning. In the TDOA positioning problem of the MEOSAR system, the PSO algorithm can find the positioning solution through an iterative optimization process, achieving a good balance between computational efficiency and accuracy.

[0003] However, the PSO algorithm cannot be directly used for passive satellite positioning because it is a nonlinear optimization problem constrained by the Earth's surface. PSO's update mechanism relies on the particle's position, velocity, individual optimal solution, and global optimal solution. Without constraints, the particle's movement depends entirely on these factors, generating solutions that do not meet the constraints, leading to 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 the MEOSAR system, improving 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 attempted to further improve the accuracy and reliability of TDOA positioning through methods such as algorithm optimization and the introduction of constraints. In summary, improving the existing PSO algorithm to adapt to the special constraints on the Earth's surface is of great significance and application value for enhancing the performance of the MEOSAR system.

[0005] Therefore, proposing a method to improve the positioning accuracy of the MEOSAR system TDOA in order to solve the difficulties existing in the prior art is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

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

[0007] To achieve the above objectives, the present invention provides the following technical solution:

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

[0009] Analyze precise ephemeris data to obtain satellite positions and velocities, and construct a satellite position matrix and a time difference matrix for signals relayed by multiple satellites.

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

[0011] The fitness function is used to evaluate the quality of particle positions, and the optimal position of each individual particle and the global optimal position are obtained.

[0012] Update the particle position within the constraint space, use boundary conditions to restrict the particle position update to satisfy the solution space, and iterate.

[0013] Update the parameters and check if 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 details for obtaining the satellite position are as follows:

[0015] The steps to obtain a satellite's position 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, GM is the product of the gravitational constant G and the total mass M of the Earth, 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 observation time;

[0021] Calculate the satellite's mean perihelion angle at the time of signal transmission:

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

[0023] t = t 1 -△t

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

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

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

[0027] Where M0 is the mean anterior angle at the reference time, M k The mean anomaly of the satellite at the time of signal transmission;

[0028] Calculate the near-point angle:

[0029] E k =M k +esinE k

[0030] Calculate the true anterior angle:

[0031]

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

[0033] Calculate the ascending distance angle:

[0034]

[0035] Where ω is the angular distance from the nearest point;

[0036] Calculate the perturbation correction term:

[0037]

[0038]

[0039] Among them, C uc C us C rc C rs C ic C is There are 6 perturbation correction parameters, δ uk δ is the perturbation correction term for the ascending intersection angle. rk For satellite radius perturbation correction term, δ ikThis is a correction term for satellite orbital inclination perturbation.

[0040] Perturbation-corrected ascending distance angle, satellite radius vector, and orbital inclination:

[0041]

[0042] Where i0 is t oe The orbital inclination at any given time is given by the Kepler-six parameters of the broadcast ephemeris; The rate of change of the orbital inclination i is given by the perturbation nine parameters in the broadcast ephemeris; A is the major radius of the satellite orbit; This is a satellite radius perturbation correction term;

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

[0044] x′ k =r k cosu k

[0045] y′ k =r k sinu k

[0046] Where, x′ k Let y′ be the projected coordinates of the satellite on the orbital plane along the horizontal axis of the orbital plane. k These are the projected coordinates 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, Let Ω0 be the rate of change of the ascending node with respect to time, and let Ω0 be the right ascension of the ascending node at the reference time. Provided by satellite ephemeris, Ω k To calculate time t k The right ascension of the ascending node, Let t be the Earth's rotational angular velocity. k t is the time for calculating the satellite's position. oe Reference time for satellite ephemeris;

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

[0051]

[0052] Optionally, the specific details for obtaining satellite velocity are as follows:

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

[0054] Calculate the first derivatives of the mean anterior angle and the deviated anterior angle with respect to time:

[0055]

[0056] in, The first derivative of the angle at the point of approach with respect to time is... is the first derivative of the true anterior angle with respect to time;

[0057] The first derivative of the true anterior angle is:

[0058]

[0059] First derivative of the angular distance from the ascending nodes:

[0060]

[0061] in, Let be the first derivative of the true anomaly angle with respect to time, representing the rate of change of the satellite's true position angle in its orbit over time. The first derivative of the angular distance from the ascending node with respect to time;

[0062] First derivative of the perturbation correction term:

[0063]

[0064] in, Let be the first derivative of the perturbation correction term of the ascending intersection angle with respect to time. The first derivative of the satellite radius perturbation correction term with respect to time. This is the first derivative of the satellite orbital inclination correction term with respect to time.

[0065] First derivatives of ascending node distance, satellite radius vector, and orbital inclination after perturbation correction:

[0066]

[0067] in, Let be the first derivative of the angular distance of the ascending node with respect to time, representing the rate of change of the angular distance of the ascending node with time. Let be the first derivative of the satellite's radius vector, representing the rate of change of the satellite's distance from the Earth's center with time. Let be the first derivative of the orbital inclination with respect to time, representing the rate of change of the orbital inclination over time. This refers to the orbital inclination, which is the angle between the orbital plane and the Earth's equatorial plane. It is the first derivative with respect to time, representing the rate at which the orbital inclination changes with time;

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

[0069]

[0070] in, Let be the first derivative of the satellite's projected coordinates along the horizontal axis of the orbital plane with respect to time, i.e., the velocity component of the satellite along the horizontal axis of the orbital plane. It is the first derivative of the satellite's projected coordinates along the longitudinal axis of the orbital plane with respect to time, i.e., 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, The first derivative of the right ascension of the ascending node. The angular velocity at the right ascension of the ascending node refers to the precession velocity of the satellite's orbit. The rate of change of right ascension caused by the Earth's rotation;

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

[0075]

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

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

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

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

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

[0081] Where 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 iLet be the position coordinates of the i-th satellite, r be the position coordinates of the ground station, g(y) be the position coordinates of the distress beacon, and s be the position coordinates of the satellite. j Let be the position coordinates of the j-th satellite, and c be the speed of light. TDOA ij Let be the time difference between the arrival time of the distress signal after passing through the i-th satellite and the -th satellite, respectively, at the ground station.

[0085] Optionally, the particle position is updated within the constraint space. Boundary conditions are used to restrict the particle position update to satisfy the solution space. The specific details of the iteration are 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) represents the velocity vector of particle i in the t-th iteration, w is the inertia weight, controlling the particle's inertia, i.e., the influence of the previous velocity on the current velocity, c1 is the cognitive coefficient, i.e., the individual learning factor, controlling the degree to which the particle approaches its own historical best position, c2 is the social coefficient, i.e., the group learning factor, controlling the degree to which the particle approaches the global best position, r1 and r2 are random numbers between [0, 1], increasing the randomness and diversity of the search, pBest i Let gBest be the best historical position of particle i in the t-th iteration, and gBest be the best global historical position, representing the best position of all particles in the swarm in the t-th iteration.

[0089] The particle position update formula is:

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

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

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

[0093] 1) Compared with the traditional PSO algorithm, it improves the calculation accuracy, significantly reduces the positioning error, and achieves better positioning results;

[0094] 2) By constraining the search space of particles within the range of restrictions, the ineffective exploration of the particle swarm is greatly reduced, saving resources and enabling particles to find the optimal solution faster, which greatly speeds up the convergence speed of the algorithm and has higher efficiency in practical applications.

[0095] 3) It has stronger noise resistance and adaptability, can significantly improve positioning accuracy when the signal-to-noise ratio is low, and continuously optimizes the results as the signal quality gradually improves. It shows stronger robustness and stability than the standard PSO algorithm, especially in the environment with weak signal, where it shows obvious advantages. Attached Figure Description

[0096] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

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

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

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

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

[0101] Figure 5 This invention provides a comparison chart of positioning errors under different signal-to-noise ratio conditions. Detailed Implementation

[0102] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

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

[0104] Analyze precise ephemeris data to obtain satellite positions and velocities, and construct a satellite position matrix and a time difference matrix for signals relayed by multiple satellites.

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

[0106] The fitness function is used to evaluate the quality of particle positions, and the optimal position of each individual particle and the global optimal position are obtained.

[0107] Update the particle position within the constraint space, use boundary conditions to restrict the particle position update to satisfy the solution space, and iterate.

[0108] Update the parameters and check if 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 to this method is the objective function to be optimized, along with various parameters including the number of particles, learning factor, maximum number of iterations, optimization objective dimension, satellite position matrix, ground station position, and arrival time difference matrix. The output is the global optimal position of the particle swarm (B, L, H), i.e., the distress beacon solution position in the geodetic coordinate system.

[0110] In the MEOSAR system, signals emitted by distress beacons travel through different medium-Earth orbit satellites to reach ground receiving stations. Ground station antennas receive these signals, record and compare their arrival times to obtain the time difference, which reflects the distance differences between the signal from the beacon and the different satellites. A set of time differences defines a hyperboloid, and the location of the signal source lies on this hyperboloid. Multiple sets of time differences can define multiple hyperboloid equations; therefore, three or more time difference equations can uniquely determine the location of the signal source.

[0111] Considering the equations in a geocentric Earth-fixed coordinate system, the arrival time observation equation for signals relayed by satellite i can be described as follows:

[0112]

[0113] Among them, (x i y i , z i ) for satellite S i The position coordinates, (x, y, z) are the position coordinates of the beacon signal, (x... g y g , z gLet ) represent the position coordinates of the ground receiving station, T0 represent the transmission time of the beacon signal (unknown), Δt represent the lead of the beacon clock relative to the ground station clock (unknown), and δL represent the position coordinates of the ground receiving station. i These are various error parameters, including ionospheric error, tropospheric error, and multipath effect error.

[0114] Subtracting the arrival time observation equations of two different satellites yields a set of arrival time difference equations, which can eliminate the unknowns Δt and T0, and also offset a large portion of the time error and errors caused by multipath effects.

[0115]

[0116] TDOA ij =TOA i -TOA j

[0117] Solving different time difference equations and the Earth's elliptic equation simultaneously: This allows us to determine the location of the distress beacon.

[0118] In the MEOSAR problem, the goal is to find the beacon location. Therefore, the construction of the fitness function should revolve around minimizing the localization error, and can be set as follows:

[0119]

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

[0121] The Earth is an approximate ellipsoid, and limiting latitude, longitude, and altitude within a specific range is intuitive in a geodetic coordinate system. There is an elevation difference between the Earth elliptical reference model and the actual Earth surface due to topographical features such as mountains, oceans, and other terrain features. Ignoring extreme terrain such as high mountains and deep valleys, the average elevation difference between most of the global land and ocean areas and the ellipsoid is typically around ±100 meters. Therefore, the ranges for latitude, longitude, and altitude are set as follows:

[0122] -90 ° ≤B≤90 °

[0123] -180 ° ≤L≤180 °

[0124] -100≤H≤100

[0125] It can satisfy the constraints for particle swarms to fly on or near the Earth's surface.

[0126] When evaluating particle positions, it is necessary to transform the spatial geodetic coordinate system to the Cartesian coordinate system before substituting the values ​​into the fitness function:

[0127]

[0128] The individual's optimal position and the global optimal position at the current iteration number are determined based on the fitness value, and the process 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] Wherein, the velocity of particle i at time t is v i (t), at position x i (t), the individual's optimal position is pBest i The globally optimal position is gBest.

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

[0133] The settings are dynamically adjusted during iteration, and gradually decreased as the number of iterations increases:

[0134]

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

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

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

[0138] During the update process, the following strategy is adopted to restrict particles to update within the constraints, ensuring that the algorithm uses 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, it outputs the currently found global best position and its corresponding fitness value as the optimal solution to the problem.

[0143] For details, see Figure 2 The diagram shown is a schematic of the particle confinement space.

[0144] Furthermore, the specific details for obtaining the satellite position are as follows:

[0145] The steps to obtain a satellite's position 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, GM is the product of the gravitational constant G and the total mass M of the Earth, 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 observation time;

[0151] Calculate the satellite's mean perihelion angle at the time of signal transmission:

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

[0153] t = t 1 -△t

[0154] Among them, t 1 t represents the observation time. ocHere, t is the reference time, Δt is the time corrected by the satellite clock, Δt is the satellite clock error, which represents the total deviation of the satellite clock relative to standard time (such as GPS system time or other reference time), a0 is the satellite clock deviation, a1 is the satellite clock drift (representing the rate of change of the satellite per unit time), and a2 is the satellite clock drift rate (representing the acceleration of clock drift over time, usually a second-order correction term used to compensate for nonlinear drift).

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

[0156] Where M0 is the mean anterior angle at the reference time, M k The mean anomaly of the satellite at the time of signal transmission;

[0157] Calculate the near-point angle:

[0158] E k =M k +esinE k

[0159] Calculate the true anterior angle:

[0160]

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

[0162] Calculate the ascending distance angle:

[0163]

[0164] Where ω is the angular distance from the nearest point;

[0165] Calculate the perturbation correction term:

[0166]

[0167] Among them, C uc C us C rc C rs C ic C is There are 6 perturbation correction parameters, δ uk δ is the perturbation correction term for the ascending intersection angle. rk For satellite radius perturbation correction term, δ ik This is a correction term for satellite orbital inclination perturbation.

[0168] Perturbation-corrected ascending distance angle, satellite radius vector, and orbital inclination:

[0169]

[0170]

[0171]

[0172] Where i0 is t oe The orbital inclination at any given time is given by the Kepler-six parameters of the broadcast ephemeris; The rate of change of the orbital inclination i is given by the perturbation nine parameters in the broadcast ephemeris; A is the major radius of the satellite orbit; This is a satellite radius perturbation correction term;

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

[0174] x′ k =r k cosu k

[0175] y′ k =r k sinu k

[0176] Where, x′ k Let y′ be the projected coordinates of the satellite on the orbital plane along the horizontal axis of the orbital plane. k These are the projected coordinates 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 a satellite, which is uniquely determined by the satellite's orbital elements (six Kepler parameters). These two coordinates describe the two-dimensional projection position of the satellite within the orbital plane, and their physical meaning is the satellite's components along the orbital plane in two directions.

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

[0179]

[0180] in, Let Ω0 be the rate of change of the ascending node with respect to time, and let Ω0 be the right ascension of the ascending node at the reference time. Provided by satellite ephemeris, Ω k To calculate time t k The right ascension of the ascending node, The Earth's rotational angular velocity (approximately 7.2921151467 × 10⁻⁶) 5 radians / second), t k t is the time for calculating the satellite's position. oe Reference time for satellite ephemeris;

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

[0182]

[0183] Furthermore, the specific details for obtaining satellite velocity are as follows:

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

[0185] Calculate the first derivatives of the mean anterior angle and the deviated anterior angle with respect to time:

[0186]

[0187] in, The first derivative of the angle at the point of approach with respect to time is... is the first derivative of the true anterior angle with respect to time;

[0188] The first derivative of the true anterior angle is:

[0189]

[0190] First derivative of the angular distance from the ascending nodes:

[0191]

[0192] in, Let be the first derivative of the true anomaly angle with respect to time, representing the rate of change of the satellite's true position angle in its orbit over time. The first derivative of the angular distance from the ascending node with respect to time;

[0193] First derivative of the perturbation correction term:

[0194]

[0195] in, Let be the first derivative of the perturbation correction term of the ascending intersection angle with respect to time. The first derivative of the satellite radius perturbation correction term with respect to time. This is the first derivative of the satellite orbital inclination correction term with respect to time.

[0196] First derivatives of ascending node distance, satellite radius vector, and orbital inclination after perturbation correction:

[0197]

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

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

[0200]

[0201] in, Let be the first derivative of the satellite's projected coordinates along the horizontal axis of the orbital plane with respect to time, i.e., the velocity component of the satellite along the horizontal axis of the orbital plane. It is the first derivative of the satellite's projected coordinates along the longitudinal axis of the orbital plane with respect to time, i.e., 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, The first derivative of the right ascension of the ascending node. The angular velocity at the right ascension of the ascending node refers to the precession velocity of the satellite's orbit. The rate of change of right ascension caused by the Earth's rotation;

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

[0206]

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

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

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

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

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

[0212] Where 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 Let be the position coordinates of the i-th satellite, r be the position coordinates of the ground station, g(y) be the position coordinates of the distress beacon, and s be the position coordinates of the satellite. j Let be the position coordinates of the j-th satellite, and c be the speed of light. TDOA ij Let be the time difference between the arrival time of the distress signal after passing through the i-th satellite and the -th satellite, respectively, at the ground station.

[0216] Specifically, all the coordinates mentioned above are in the Cartesian coordinate system. f = g(y) is the transformation function for converting spatial geodetic coordinates to Cartesian coordinates, as shown in the following formula:

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

[0218]

[0219] in, denoted as the radius of curvature of the ramid; f = (x, y, z) represents the geocentric coordinates in the Cartesian coordinate system; y = (B, L, H) represents the three-dimensional coordinates in the geodetic ellipsoidal coordinate system; and e represents the ellipsoidal ellipticity.

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

[0221] Furthermore, updating the particle positions within the constraint space, using boundary conditions to restrict the particle position updates to satisfy the solution space, and performing iterations are as follows:

[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) represents the velocity vector of particle i in the t-th iteration, w is the inertia weight, controlling the particle's inertia, i.e., the influence of the previous velocity on the current velocity, c1 is the cognitive coefficient, i.e., the individual learning factor, controlling the degree to which the particle approaches its own historical best position, c2 is the social coefficient, i.e., the group learning factor, controlling the degree to which the particle approaches the global best position, r1 and r2 are random numbers between [0, 1], increasing the randomness and diversity of the search, pBest iLet gBest be the best historical position of particle i in the t-th iteration, and gBest be the best global historical position, representing the best position of all particles in the swarm in the t-th iteration.

[0225] The particle position update formula is:

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

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

[0228] Specifically, after each update, it checks whether the three updated parameters of the particle are within the constraint range. It also checks whether the latitude and longitude are limited to the constraint range using modulo operations, thus achieving cyclic mapping. This eliminates the need for complex logical judgments and ensures the continuity of parameters. No matter how large or small the input value is, it will be smoothly mapped to the corresponding constraint range without any jumps or breakpoints, which helps to optimize the smooth operation of the algorithm.

[0229] Moreover, latitude and longitude are inherently periodic (a characteristic of Earth's spherical coordinates). This method can make good use of this characteristic, cycling out-of-range values ​​back into the legal range, rather than 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 can be 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 the simulation, a fixed error of 3 microseconds was set for calculating ephemeris data, and the TDOA measurement error was set to 0.3 microseconds to match real-world measurement conditions. The TDOA algorithm was simulated a total of 100 times, with each run using the particle swarm optimization algorithm iterating 300 times and employing 150 particles. Different signal-to-noise ratios will affect the TDOA measurement values ​​and thus the positioning error; therefore, please refer to [further details]. Figure 5 The simulation shows the positioning error under different signal-to-noise ratios. However, the pattern observed at different signal-to-noise ratios is the same; see [link to simulation]. Figure 3 As shown, 3a and 3b present simulation results under a 10dB signal-to-noise ratio condition, which better reflects real-world signal-to-noise ratio conditions. Traditional algorithms do not impose longitude, latitude, or altitude constraints. See also... Figure 4 As shown, 4a and 4b represent the convergence curves of the two algorithms under a signal-to-noise ratio of 10 dB.

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

[0242] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those 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 invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for improving the positioning accuracy of a MEOSAR system TDOA, characterized in that, Includes the following steps: Analyze precise ephemeris data to obtain satellite positions and velocities, and construct a satellite position matrix and a time difference matrix for signals relayed by multiple satellites. Initialize the particle swarm in a spatial geodetic coordinate system, and set the particle position optimization parameters and the corresponding initialization range of the parameters; The fitness function is used to evaluate the quality of particle positions, and the optimal position of each individual particle and the global optimal position are obtained. Update the particle position within the constraint space, use boundary conditions to restrict the particle position update to satisfy the solution space, and iterate. Update the parameters and check if the maximum number of iterations has been reached. If not, continue iterating. When the algorithm ends, output the global optimal solution. The specific steps for initializing the particle swarm in a geodetic spatial coordinate system, including setting the particle position optimization parameters and their corresponding initialization ranges, are as follows: Set the particle position optimization parameters to B, L, and H, which represent latitude, longitude, and altitude, respectively. x(i,1)=180 rand () - 90; % initialize latitude in [-90, 90] range; x(i,2)=360 rand()-180; %Initialize longitude within the range [-180, 180]; x(i,3)=0.2 rand()-0.1; % Initializes the height to the range of [-100, 100] meters; Where 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; The fitness function is: in, For the first i The position coordinates of the satellite, These are the location coordinates of the ground station. The coordinates of the distress beacon. For the first j The position coordinates of the satellite, The distress signal passed through the first i satellite and the first j The time difference between the arrival of each satellite at the ground station; The process involves updating particle positions within a constrained space, using boundary conditions to ensure that the updated particle positions satisfy the solution space, and performing the following iterations: The particle velocity is updated as follows: in, For particles In the The velocity vector of the next iteration Inertial weights control the particle's inertia, i.e., the influence of the previous step's velocity on the current velocity. The cognitive coefficient, or individual learning factor, controls how closely a particle moves closer to its historical best position. The social coefficient, or group learning factor, controls how much a particle approaches the globally optimal position. 、 for Random numbers between For particles In the The best historical position in the next iteration The best position in the global history represents the position of all particles in the swarm at the [missing information]. The optimal position in the next iteration; The particle position update formula is: in, For particles In the The position vector in the next iteration. For particles In the The velocity vector in the next iteration.

2. The method of claim 1, wherein the method is used in a MEOSAR system. The specific details for obtaining satellite positions are as follows: The steps to obtain a satellite's position using orbital parameters are as follows: Calculate the average angular velocity of the satellite: wherein is the average angular velocity of the earth at the reference time , is the product of the gravitational constant G and the total mass of the earth M , ; wherein, the perturbation correction given for broadcast ephemeris, is the mean angular velocity of the satellite at the observation time; Calculate the satellite's mean perihelion angle at the time of signal transmission: wherein, tobs is the observation time, tref is the reference time, tss is the time corrected for satellite clock, δsc is the satellite clock error, δsc is the satellite clock bias, δsc is the satellite clock drift, δsc is the satellite clock drift rate; wherein is the mean anomaly at the reference time, is the mean anomaly at the reference time, Calculate the near-point angle: Calculate the true anterior angle: wherein, is the eccentricity of the satellite orbit; Calculate the ascending distance angle: wherein is the close point angle distance; Calculate the perturbation correction term: wherein, 、 、 、 、 、 are 6 perturbation correction parameters, is a perturbation correction term of the ascending node angle, is a perturbation correction term of the satellite's radial, is a perturbation correction term of the satellite's inclination. Perturbation-corrected ascending distance angle, satellite radius vector, and orbital inclination: wherein is the mean motion of the satellite; is the inclination of the orbit at the time instant, given by the Keplerian six parameters of the broadcast ephemeris; is the inclination of the orbit is the rate of change of the inclination of the orbit, given by the perturbation nine parameters in the broadcast ephemeris; is the semi-major axis of the satellite orbit; is the satellite radial perturbation correction term; Calculate the satellite's coordinates in the orbital plane coordinate system: wherein is the projection coordinate of the satellite on the orbital plane along the transversal axis of the orbital plane, is the projection coordinate of the satellite on the orbital plane along the longitudinal axis of the orbital plane; computing the right ascension of the ascending node : wherein is the rate of change of the ascending node with respect to time, is the right ascension of the ascending node at the reference time, and is provided by the satellite ephemeris, is the right ascension of the ascending node at the computation time is the right ascension of the ascending node at the reference time, is the angular velocity of the Earth rotation, is the computation time of the satellite position, is the reference time of the satellite ephemeris; Calculate the satellite's coordinates in the Earth-fixed coordinate system: 。 3. The method of claim 1, wherein the method is a method for improving the positioning accuracy of a MEOSAR system TDOA, characterized in that, The specific details for obtaining satellite velocity are as follows: The steps to obtain satellite velocity using orbital parameters are as follows: Calculate the first derivatives of the mean anterior angle and the deviated anterior angle with respect to time: wherein is the first derivative of the mean anomaly with respect to time, is the first derivative of the true anomaly with respect to time; The first derivative of the true anterior angle is: First derivative of the angular distance from the ascending node: wherein is the first derivative of the true anomaly with respect to time, representing the rate of change of the true position angle of the satellite on the orbit with respect to time, is the first derivative of the ascending node angle with respect to time; First derivative of the perturbation correction term: wherein, is the first order derivative of the perturbation correction term for the right ascension with respect to time, is the first order derivative of the perturbation correction term for the satellite's argument of latitude with respect to time, is the first order derivative of the perturbation correction term for the satellite's inclination with respect to time; First derivatives of ascending node distance, satellite radius vector, and orbital inclination after perturbation correction: in, Let be the first derivative of the angular distance of the ascending node with respect to time, representing the rate of change of the angular distance of the ascending node with time. Let be the first derivative of the satellite's radius vector, representing the rate of change of the satellite's distance from the Earth's center with time. Let be the first derivative of the orbital inclination with respect to time, representing the rate of change of the orbital inclination over time. This refers to the orbital inclination, which is the angle between the orbital plane and the Earth's equatorial plane. The first derivative with respect to time represents the rate at which the orbital inclination changes with time; Calculate the satellite's velocity in the orbital plane: wherein is the first derivative of the satellite's projected coordinate on the orbital plane along the transversal axis of the orbital plane with respect to time, i.e. the velocity component of the satellite along the transversal axis of the orbital plane, is the first derivative of the satellite's projected coordinate on the orbital plane along the longitudinal axis of the orbital plane with respect to time, i.e. 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: wherein, is the first derivative of the right ascension of the ascending node, is the angular velocity of the right ascension of the ascending node, indicating the precession speed of the satellite orbit, is the speed of the change in right ascension caused by the rotation of the Earth; Calculate the satellite's velocity in the ECEF coordinate system: 。

Images