Sky wave over-the-horizon radar high-speed parallel coordinate registration method based on reference source assistance

Abstract

The invention relates to a sky wave over-the-horizon radar high-speed parallel coordinate registration method based on reference source assistance, and the method comprises the steps: constructing and estimating the posterior distribution of ionosphere activity parameters under a Bayesian framework based on the prior information of the ionosphere activity parameters and reference source data; secondly, according to posterior distribution of the ionosphere activity parameters, based on a near-end gradient descent method and an average field variational method, solving estimated values of the ionosphere activity parameters of all dimensions in parallel, and then substituting the estimated values of the ionosphere activity parameters of all dimensions into an international reference ionosphere model to obtain a posterior ionosphere; and finally, calculating an electromagnetic wave propagation path in the posterior ionosphere by using a three-dimensional ray tracing method to obtain a coordinate registration result of the to-be-positioned target. According to the method, timeliness can be improved while the accuracy of sky wave radar target positioning is ensured.

Description

Technical Field

[0001] The present invention relates to the technical field of over-the-horizon radar data processing, and in particular to a reference source-assisted skywave over-the-horizon radar high-speed parallel coordinate registration method. Background Art

[0002] Over-the-Horizon Radar (OTHR) is a radar system that exploits the ionosphere's reflection of high-frequency radio waves to detect targets beyond the horizon. It is a crucial component of modern military long-range early warning systems and environmental awareness systems. Its advantages include long detection range, wide coverage, strong anti-stealth capabilities, and high cost-effectiveness, making it a crucial component of long-range early warning systems.

[0003] Due to the operating principle of OTHR, radar performance is decisively influenced by the Earth's atmospheric ionosphere. Observation, modeling, and analysis of the ionosphere are prerequisites for radar coordinate registration and target tracking. Radar coordinate registration is the process of deriving target geographic information from the target's position information received by the radar in its radar coordinate system. OTHR error analysis shows that errors introduced by ionospheric activity are the primary source of coordinate registration errors. By introducing reference source information to assist in ionospheric activity prediction and formulating the ionospheric parameter inference task as a Bayesian inverse problem, coordinate registration performance can be significantly improved in the absence of conventional ionospheric detection methods.

[0004] However, the existing methods of introducing reference source information to assist in ionospheric activity prediction usually adopt a pure sampling approach. The pure sampling approach has a large amount of computation when solving high-dimensional ionospheric parameter estimation problems, and it is difficult to use modern computing platforms to improve computational efficiency, resulting in poor prediction timeliness. For example, the Chinese invention patent application with application number CN202310729538.5 and publication number CN117008110A discloses a reference source-assisted skywave over-the-horizon radar coordinate registration method. This method constructs the posterior distribution of the ionospheric activity parameters in the ionospheric region to be estimated under a Bayesian framework based on the prior information of the ionospheric activity parameters in the ionospheric region to be estimated and the reference source data. Then, based on the Markov chain, Gibbs sampling and importance sampling, the moment estimate of the posterior distribution of the ionospheric activity parameters in the ionospheric region to be estimated is solved. Then, the moment estimate of the posterior distribution of the ionospheric activity parameters in the ionospheric region to be estimated is substituted into the international reference ionosphere model to obtain the posterior ionosphere. Finally, the three-dimensional ray tracing method is used to calculate the electromagnetic wave propagation path in the posterior ionosphere to obtain the coordinate registration result of the target to be located. This method can significantly improve the accuracy of coordinate registration, thereby improving the accuracy of skywave over-the-horizon radar target positioning. However, in order to obtain a sufficient number of samples to estimate the posterior distribution of ionospheric activity parameters, this method requires repeated execution of the Markov chain walk step. This process requires repeated calls to the forward model, which uses the International Reference Ionosphere Model to update the electron concentration profile and uses ray tracing technology to calculate the propagation paths of a large number of radar signals in the ionosphere, making it difficult to ensure timeliness.

[0005] In general, the existing methods that use pure sampling methods to introduce reference source information to assist in ionospheric activity prediction have high computational complexity and poor timeliness when solving high-dimensional ionospheric parameter estimation problems. Summary of the Invention

[0006] Based on this, it is necessary to provide a high-speed parallel coordinate registration method for skywave over-the-horizon radar based on reference source assistance to address the above technical problems, which can improve timeliness while ensuring the accuracy of skywave radar target positioning.

[0007] The present invention provides a high-speed parallel coordinate registration method for skywave over-the-horizon radar based on reference source assistance, comprising the following steps:

[0008] The posterior distribution of ionospheric activity parameters is constructed and estimated based on the prior information of ionospheric activity parameters and reference source data in a Bayesian framework.

[0009] According to the posterior distribution of ionospheric activity parameters, the estimated values ​​of ionospheric activity parameters in all dimensions are solved in parallel based on the proximal gradient descent method and the mean field variation method;

[0010] The estimated values of the ionospheric activity parameters in all dimensions are substituted into the international reference ionosphere model to obtain a posterior ionosphere;

[0011] The electromagnetic wave propagation path is calculated in the posterior ionosphere by using a three-dimensional ray tracing method to obtain the coordinate registration result of the target to be positioned.

[0012] In one embodiment, the estimation of the ionospheric activity parameters in all dimensions is solved based on the posterior distribution of the ionospheric activity parameters according to the proximal gradient descent method and the mean field variational method, including the following steps:

[0013] A variational distribution family of the ionospheric activity parameters is selected under the mean field assumption;

[0014] An evidence lower bound function is constructed based on the posterior distribution of the ionospheric activity parameters and the variational distribution family;

[0015] The estimated values of the ionospheric activity parameters in each dimension are solved based on the evidence lower bound function in parallel;

[0016] The estimated values of the ionospheric activity parameters in all dimensions are obtained by combining the estimated values of the ionospheric activity parameters in each dimension.

[0017] In one embodiment, the estimated value of the ionospheric activity parameter in any dimension is solved, including the following steps:

[0018] The evidence lower bound function is decomposed and reorganized into a function composed of a first data item, a second data item and a third data item;

[0019] The proximal gradient descent method is used to optimize the first data item and the second data item to obtain the optimal variational distribution of the ionospheric activity parameter in the current dimension;

[0020] The ionospheric activity parameter is estimated according to the optimal variational distribution to obtain the estimated value of the ionospheric activity parameter in the current dimension;

[0021] The first data item is a data item that only uses a numerical sampling method to calculate gradient information, the second data item is a data item that can use an analytical method to calculate gradient information, and the third data item is a data item that does not participate in the optimization process.

[0022] In one embodiment, the evidence lower bound function is decomposed and reorganized into a function composed of a first data item, a second data item and a third data item, including the following steps:

[0023] The evidence lower bound function is decomposed by mean field approximation;

[0024] All dimensions obtained by mean field approximation are reorganized under the framework of proximal gradient optimization to obtain a function composed of a first data item, a second data item and a third data item.

[0025] In one embodiment, optimizing the first data item and the second data item using a proximal gradient descent method includes the following steps:

[0026] Calculating the gradient information of the first data item using the Monte Carlo method;

[0027] Using the quasi-Monte Carlo method to control the variance of the gradient information estimation of the first data item;

[0028] Establish a variational distribution parameter iteration formula within the framework of proximal gradient optimization;

[0029] The variational distribution parameters are calculated iteratively until the convergence condition is met to obtain the optimal variational distribution of ionospheric activity parameters.

[0030] In one embodiment, the function composed of the first data item, the second data item, and the third data item is:

[0031]

[0032] Where, is the first data item, is the second data item, elseterm is the third data item, x represents the ionospheric activity parameter vector, z represents the radar measurement vector of the reference source and the ground distance vector, x i represents the ionospheric activity parameter of the i-th dimension, x -i represents the ionospheric activity parameter vector of the remaining dimensions except the i-th dimension, λ represents the total variation distribution parameter, λ i Represents the variational parameter of the variational distribution of the current i-th dimension, q λ (x) is the variational distribution of all ionospheric activity parameters, is the variational distribution of the ionospheric activity parameter in the i-th dimension, logp(z|x) is the log-likelihood function, p(x i ,z -i ) is the prior distribution of all ionospheric activity parameters.

[0033] In one embodiment, the variational distribution parameter iteration formula is:

[0034]

[0035] Where, It is a hyperparameter introduced in the proximal gradient descent method, which is used to characterize the strong constraint effect of the proximal term on the optimization process. is the damping parameter, represents the variational distribution of the ionospheric activity parameter of the i-th dimension at the k-th iteration, g(λ k ) represents the first data item of the k-th iteration;

[0036] Among them, for a certain choice of variational distribution family, the following form of variational parameter iterative expression with damping form can be obtained:

[0037]

[0038] In the formula, for a certain choice of variational distribution family, the following form of variational parameter iteration expression with damping form can be obtained:

[0039]

[0040] Where, represents the variational parameter of the variational distribution of the i-th dimension at the k+1-th iteration, represents the variational parameter of the variational distribution of the i-th dimension at the k-th iteration, and l represents the learning rate.

[0041] The beneficial effects of the present invention are as follows: the present invention combines the mean field variational method, which is a deterministic method, with the sampling method to reduce the instability in the sampling process and improve computational efficiency. The introduction of the proximal gradient optimization method establishes an iterative update formula for ionospheric parameters to solve the estimated values ​​of ionospheric activity parameters in all dimensions. This allows for the parallel solution of ionospheric activity parameters in multiple dimensions, improving the parallel adaptability of solving ionospheric parameters, thereby improving the computational efficiency of the coordinate alignment method, and ensuring the timeliness of skywave radar coordinate alignment. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 This is one of the flow charts of a method for high-speed parallel coordinate registration of a skywave over-the-horizon radar based on reference source assistance provided by an embodiment of the present invention;

[0043] Figure 2 Schematic diagram of the ionospheric range and reference source distribution provided by an embodiment of the present invention;

[0044] Figure 3 Schematic diagram of the estimation accuracy of ionospheric parameters in the target area provided by an embodiment of the present invention;

[0045] Figure 4 Schematic diagram of the coordinate registration results of 12 targets using two methods provided in an embodiment of the present invention;

[0046] Figure 5 yes Figure 4 A magnified diagram of the coordinate registration results of targets 1-6;

[0047] Figure 6 yes Figure 4 A magnified diagram of the coordinate registration results of targets 7-12. DETAILED DESCRIPTION

[0048] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0049] In one embodiment, in one embodiment, Figure 1 As shown, Figure 1 This is one of the flow charts of the high-speed parallel coordinate registration method for skywave over-the-horizon radar based on reference source assistance provided by an embodiment of the present invention. The high-speed parallel coordinate registration method for skywave over-the-horizon radar based on reference source assistance of this embodiment includes the following steps:

[0050] S1. Construct and estimate the posterior distribution of ionospheric activity parameters based on their prior information and reference source data in a Bayesian framework.

[0051] It should be noted that the ionospheric activity parameters are the controlling parameters of the ionospheric electron concentration profile, and the reference source is a reference target with known geographical location information within the ionospheric region to be estimated. The reference source data includes the radar measurement vector and ground distance vector of the reference source.

[0052] In this embodiment, the posterior distribution of the ionospheric activity parameters in the ionospheric region to be estimated is:

[0053]

[0054] In the formula, y=(y 1 ,…,y r ),z=(z 1 ,…,z r ) represent the radar measurement vector and ground distance vector of the reference source, x is the ionospheric activity parameter vector in the ionospheric region to be estimated, r is the number of reference sources, β=(β 1 ,…,β r ) is the elevation angle information of each ray, ψ(x,β r ) and φ(x,β r ) represents the forward model, and π(x, β) represents the joint prior distribution of (x, β). For simplicity, the present invention will abbreviate the posterior distribution of ionospheric activity parameters as p(x|z), where x represents all ionospheric activity parameter vectors and z represents the radar measurement vector representing the reference source and the ground distance vector. Accordingly, p(x, z) is used to represent the joint distribution, and p(x) is used to represent the prior distribution of the latent variable.

[0055] S2. Based on the posterior distribution of ionospheric activity parameters, the estimated values ​​of ionospheric activity parameters in all dimensions are solved in parallel using the proximal gradient descent method and the mean field variation method.

[0056] In this embodiment, the mean field variation method, which is a deterministic method, is combined with the sampling method to reduce the instability in the sampling process and improve the calculation efficiency.

[0057] In addition, the introduction of the proximal gradient optimization method establishes an iterative update formula for ionospheric parameters to solve the estimated values ​​of ionospheric activity parameters in all dimensions, which can solve the ionospheric activity parameters in multiple dimensions in parallel. Solving the estimated values ​​of ionospheric activity parameters in all dimensions in parallel can greatly optimize the timeliness of skywave radar coordinate alignment. Specifically, step S102 includes the following steps:

[0058] S21. Select the family of variational distributions of ionospheric activity parameters under the mean field assumption.

[0059] Specifically, according to the basic method of mean field approximation, the selected variational distribution family should have conditional independence. To facilitate subsequent reasoning and calculation, the selected variational distribution family is preferably the normal distribution family. The variational distribution family in this embodiment is as follows.

[0060]

[0061] Among them, q λ (·) represents the selected variational distribution family, λ is the parameter of the entire variational distribution, and m represents the number of locations where ionospheric parameters need to be estimated, that is, the total dimension of ionospheric activity parameters. is the variational distribution of the ionospheric parameters in each dimension after mean field decomposition.

[0062] S22. Construct an evidence lower bound function based on the posterior distribution and variational distribution family of ionospheric activity parameters.

[0063] According to the basic method of variational inference, the constructed evidence lower bound function ELBO can be written as:

[0064]

[0065] Where q λ (x) is the variational distribution of all ionospheric activity parameters, logp(x,z) is the joint distribution of ionospheric activity parameters and observation data (radar measurement vector and ground distance vector of the reference source), logp(z|x) is the log-likelihood function, and p(x) is the prior distribution of all ionospheric activity parameters.

[0066] S23. Based on the evidence lower bound function, the estimated values ​​of the ionospheric activity parameters in each dimension are solved in parallel.

[0067] The estimation process of one-dimensional ionospheric activity parameters is executed in parallel for all dimensions, and the estimated values ​​of ionospheric activity parameters in all dimensions are collected. This iterative process is performed as an iterative value of the estimation result of all ionospheric activity parameters until the convergence condition is met, and the estimated values ​​of ionospheric activity parameters in each dimension of the ionosphere are obtained.

[0068] In this embodiment, solving the estimated value of the ionospheric activity parameter for any dimension specifically includes the following steps:

[0069] S231. Decompose and reorganize the evidence lower bound function into a function consisting of a first data item, a second data item, and a third data item.

[0070] Specifically, we first bring the expression of the selected mean field distribution family into the evidence lower bound function to perform mean field approximate decomposition of the evidence lower bound function, and obtain the following formula:

[0071]

[0072] Then, under the framework of proximal gradient optimization, all dimensions obtained by the mean field approximation are reorganized to obtain a function consisting of the first data item, the second data item, and the third data item. In this embodiment, for the variational distribution of the i-th dimension and its parameter λ i , the function composed of the first data item, the second data item and the third data item is:

[0073]

[0074] Where, is the first data item, is the second data item, else term is the third data item, x represents the ionospheric activity parameter vector, z represents the radar measurement vector of the reference source and the ground distance vector, x i represents the ionospheric activity parameter of the i-th dimension, x -i represents the ionospheric activity parameter vector of the remaining dimensions except the i-th dimension, λ represents the total variation distribution parameter, λ i Represents the variational parameter of the variational distribution of the current i-th dimension, q λ (x) is the variational distribution of all ionospheric activity parameters, is the variational distribution of the ionospheric activity parameter in the i-th dimension, log p(z|x) is the log-likelihood function, p(x i ,x -i ) is the prior distribution of all ionospheric activity parameters.

[0075] Among them, the first data item is a data item for calculating gradient information using only a numerical sampling method, the second data item is a data item for calculating gradient information using an analytical method, and the third data item is a data item that does not participate in the optimization process.

[0076] In this embodiment, from the perspective of computational convenience, the ELBO is decomposed into two parts. It is the first data item, which includes the log-likelihood function part in the probability model. In the embodiment, the likelihood function represents the probabilistic relationship between the observed variable and the latent variable, and is modeled with the help of a complex forward model. Specifically, it includes two models: the International Reference Ionospheric Model and the High-Performance Ray Tracing Model. This makes it difficult to use analytical methods to obtain gradient information and requires the use of numerical sampling schemes. However, the numerical sampling scheme requires a large number of Monte Carlo samples to estimate the gradient information, which is inherently computationally intensive and is not conducive to improving the timeliness of reasoning. Therefore, it should be regarded as a last resort when solving the gradient information. The gradient information that can be calculated using analytical methods of ELBO is organized into the second data item, specifically including the prior distribution and variational distribution of ionospheric activity parameters. Since these two items are derived from the problem modeling of the embodiment, they are specified in a form that is easy to calculate. Therefore, it is possible to use analytical methods to solve its gradient information. The rest that is not related to optimization can be organized into the third data item.

[0077] In addition to computational difficulty, two other factors must be considered when reorganizing the ELBO. First, convexity determines the local optimum and the choice of optimization method. The ELBO function often consists of both convex and non-convex components. For convex function optimization problems, well-established solutions exist, and these solutions are generally simpler than non-convex optimization methods. Therefore, analytical methods are considered for the convex component, while sampling methods are used for the non-convex component. This is because the randomness introduced by sampling methods naturally causes the optimization results to escape the local optimum. Second, conjugation is considered. Using conjugate computation, the variational distribution inference problem can be transformed into a closed-form update involving only the variational parameters and the likelihood function, significantly reducing the computational cost of the optimization process. Numerical sampling is used to solve the non-conjugate component, and proper conjugation can avoid unnecessary computation.

[0078] S232. Optimize the first data item and the second data item using the proximal gradient descent method to obtain the optimal variational distribution of the ionospheric activity parameters in the current dimension.

[0079] Specifically, the first data item and the second data item are optimized using the proximal gradient descent method, including the following steps:

[0080] S2321. Calculate the gradient information of the first data item using the Monte Carlo method.

[0081] Specifically, for the first data item g(λ), the Monte Carlo method is used to estimate its gradient, as shown in the following formula:

[0082]

[0083] For the expected expression shown in the formula, the MCMC method is used to estimate it, as shown in the following formula.

[0084]

[0085] in, The specific update method is: based on the variational distribution λ of the previous iteration, except for the current dimension x i Other dimension variables x -i Then, from the current dimension x i Sampling a set of S particles x from the variational distribution of i,s Next, each particle x i,s with x -i Combine to form a complete latent variable sample vector x i,s .

[0086] Subsequently, the forward model is evaluated and the likelihood is calculated. Finally, the results of all particles are weighted averaged to estimate the gradient of the objective function g(λ).

[0087] S2322. Use the quasi-Monte Carlo method to control the estimated variance of the gradient information of the first data item.

[0088] Quasi-Monte Carlo (QMC) methods are used to control variance. Variance control improves the convergence speed and stability of the algorithm and reduces the number of sampling particles used in importance sampling. Quasi-Monte Carlo methods achieve variance control by replacing the random numbers in classical Monte Carlo methods with low-deviation sequences, requiring little modification to the rest of the algorithm. This makes it a simple and widely used variance control strategy.

[0089] The specific process of controlling the variance of the gradient information estimation of the first data item using the quasi-Monte Carlo method is as follows:

[0090] is the variational distribution of the current i-th dimension Select the support domain. Since it is usually modeled as a one-dimensional normal distribution, we choose its support domain as μ±σ. For other types of prior distribution selection, the support domain that can cover its high probability region is generally selected.

[0091] Generate a set of quasi-Monte Carlo points on the selected support domain. In most cases, the QMC point set on the one-dimensional support domain is essentially a one-dimensional arithmetic progression within the support domain, which can be directly obtained using the lattice rule method or the number theory grid method.

[0092] After obtaining the quasi-Monte Carlo point set, we also use the method of step S2321 to estimate the gradient information of the first data item.

[0093] S2323. Establish a variational distribution parameter iteration formula under the framework of proximal gradient optimization.

[0094] After obtaining the gradient information estimate of the first data item, the update formula of the variational distribution parameter with damping is derived under the proximal gradient optimization framework, that is, the variational distribution parameter iterative formula is established. The specific process is as follows:

[0095] The proximal gradient optimization algorithm is an iterative method for solving optimization problems with nonsmooth regularization terms. It decomposes the objective function into smooth and nonsmooth components, using gradient descent to handle the smooth term and combining it with proximal mapping to handle the nonsmooth term. It boasts high computational efficiency and a wide range of applications, making it widely used in fields such as machine learning and signal processing.

[0096] Specifically, after obtaining the gradient information estimate of the first data item, the following unconstrained optimization problem is solved.

[0097]

[0098] Where λ represents the variational distribution parameter, the superscript represents the number of iterations, the subscript represents the dimension number, and if the subscript is missing, it represents the entire variational parameter. To solve the optimization problem shown in the above formula, in this embodiment, A, B, and C are used to mark the three parts of the objective function.

[0099] When solving this optimization problem, λ k It should be considered as a parameter that is irrelevant to the optimization problem, because its estimation has been achieved by MCMC method in the process of dealing with the difficult-to-differentiate part. The rest are normally distributed.

[0100] Taking the derivative of part A,

[0101]

[0102] Taking the derivative of part B,

[0103]

[0104] Taking the derivative of C,

[0105]

[0106] Then, setting the derivative to zero yields the extreme point satisfying

[0107]

[0108] That is, the iterative form of the variational distribution parameter.

[0109] Specifically, hyperparameter selection should adhere to the following principles: In the early stages of an optimization algorithm, gradients are typically large. Introducing strong regularization at this stage helps prevent algorithm divergence. Therefore, choosing a smaller β value can achieve a greater damping effect, thereby stabilizing the training process. As the algorithm approaches a critical point, appropriately reducing the damping coefficient β allows the algorithm to more fully explore the solution space and avoid becoming trapped in a local optimum. This strategy primarily aims to prevent algorithm oscillation. In practical calculations, strategies such as a fixed step size or a step size that decreases over time can be used to select β.

[0110] Else term is the remaining terms that are irrelevant to the optimization process. Then we can get the parameter iteration formula shown in the following formula.

[0111]

[0112] Where, represents the variational parameter of the variational distribution of the i-th dimension at the k+1-th iteration, represents the variational parameter of the variational distribution of the i-th dimension at the k-th iteration, and l represents the learning rate.

[0113] For example, if the selected prior distribution is a distribution with the form of exponential distribution family p(x)=h(x)exp{η p T(x)-I(η p )}. Where h(x) is the normalization function, η p is the natural parameter of the prior distribution, T(x) is the sufficient statistics function, I(η p ) is the logarithmic normalization function. The variational distribution family also chooses the exponential distribution family with mean field properties in η q is the natural parameter of the variational distribution. Then we have the following formula

[0114]

[0115] in,[·] (i) The expression takes the i-th dimension element of the vector. Then the natural parameter The iterative expression is

[0116]

[0117] in, in[·] (ii) It means taking the i-th element on the diagonal of the matrix.

[0118] In an optional embodiment, as before, under normal distribution conditions, record

[0119] Then for the variance parameter, there is the following iterative formula

[0120]

[0121] For the mean parameter, we have

[0122]

[0123]

[0124] These iterative formulas will be used to estimate the posterior distribution of the latent variables. Furthermore, if the family of variational distributions is limited to normal distributions with the same variance parameter, we only focus on estimating the mean parameter and write the precision matrix Λ = Σ -1 , potential vector η=Λμ, then the iterative formula can be rewritten as

[0125]

[0126] After solving the inference of the ionospheric activity parameters for a particular node, the inference method for all nodes can be easily obtained in parallel using Equation (13). This is because the formula has a significant damping form, and damped update is a simple method to improve the convergence of the iterative process and is widely used in various parallel computing methods that are improved from serial methods.

[0127] S2324. Calculate the variational distribution parameters iteratively until the optimal variational distribution of the ionospheric activity parameters meets the convergence condition.

[0128] For example, the convergence condition may be that the difference between the ionospheric activity parameters obtained in two adjacent iterations is less than a preset threshold.

[0129] S233. Estimate the ionospheric activity parameters according to the optimal variational distribution to obtain estimated values ​​of the ionospheric activity parameters in the current dimension.

[0130] Under the condition of the optimal variational parameter, the expectation of the variational distribution is obtained as the estimated value of the ionospheric activity parameter. For example, the selected variational distribution family is the normal distribution, and the selected variational parameter is the potential vector of the distribution. After the reference source-assisted skywave over-the-horizon radar high-speed parallel coordinate registration method of this embodiment converges, the mean vector of the normal distribution can be solved according to formula (20), and this mean vector is used as the estimated value of the ionospheric activity parameter.

[0131] S24. Combine the estimated values ​​of the ionospheric activity parameters of each dimension to obtain the estimated values ​​of the ionospheric activity parameters of all dimensions.

[0132] S3. Substituting the moment estimate of the posterior distribution of the ionospheric activity parameters in the ionospheric region to be estimated into the international reference ionosphere model to obtain the posterior ionosphere.

[0133] The International Reference Ionosphere model (IRI) is the internationally recommended ionospheric model. The IRI is a data-driven empirical model that characterizes key ionospheric parameters. Based on a long-term record of ionospheric observations from ground and space, it provides global ionospheric characteristics and activity based on global observations and a series of empirical formulas. The core model describes the monthly averages of key parameters such as electron density, electron temperature, ion temperature, and ion composition over the global altitude range of 60 to 2000 km.

[0134] The International Reference Ionosphere Model (IRIM) divides the ionosphere into six regions based on altitude. Within each region, a set of empirical formulas based on observational data is used to model output quantities such as electron concentration as a function of ionospheric parameters and spatial and temporal position. Ionospheric parameters characterize the influence of external factors such as solar activity and Earth's magnetic field on ionospheric behavior. These include the solar index (R) (based on observations of sunspot numbers), the Ionosonde Global Index (IG) (a modified R index based on ionospheric sounder data), and the Kp (magnetic field index). The IRI has been extensively validated and is widely used in radio communications, satellite navigation systems (such as GPS), radar systems, and other fields.

[0135] S4. Calculate the electromagnetic wave propagation path in the posterior ionosphere using a three-dimensional ray tracing method to obtain the coordinate registration result of the target to be detected.

[0136] It should be noted that the use of a three-dimensional ray tracing method to perform coordinate alignment in the ionosphere is well known to those skilled in the art and will not be described in detail here.

[0137] In a specific embodiment, we focus on a specific radar coordinate registration task, which includes inferring ionospheric parameters and then using the inferred ionospheric parameters to achieve target coordinate registration.

[0138] Due to the uncertainty and complexity of ionospheric activity, ionospheric parameters have a significant impact on the coordinate registration accuracy of OTHR. This embodiment uses the International Reference Ionosphere (IRI) model to model the complex behavior of the ionosphere. The IRI model uses the solar index R based on sunspot number observations, the ionospheric altimeter global index IG, and the magnetic field index K. pIonospheric parameters are used to model each layer of the ionosphere to illustrate the influence of external factors such as solar activity and the Earth's magnetic field on the behavior of the ionosphere. This embodiment regards ionospheric activity parameters (IAPs) as hidden variables and uses OTHR measurement data and reference sources to infer these hidden variables. It should be noted that the method proposed in this embodiment can also be applied to estimate other ionospheric activity parameters. After obtaining the posterior estimates of the IAPs, the ionospheric electron density profile is calculated by substituting these IAPs into the IRI model. Based on the ionospheric electron density profile, the position of the target is calculated using ray tracing technology. The coordinate registration error between the true position of the target and the estimated position is then calculated.

[0139] The specific settings of the embodiment are as follows: In the radar monitoring area, a rectangular area from 35 degrees to 42 degrees north latitude and 130 degrees to 140 degrees east longitude is selected for estimating ionospheric activity parameters. The grid modeling step size is 0.5 degrees. The specific settings of the ionospheric range and reference source distribution are as follows: Figure 2 As shown in the figure, the radar transmitter and receiver are located at 30 degrees north latitude and 145 degrees east longitude, operating at a frequency of 15 MHz. For each dimension of the ionospheric activity parameter, the means of the true and prior distributions are set to 120 and 90, respectively. A Gaussian Markov random field prior is used to introduce a spatially correlated prior distribution for the ionospheric activity parameters. Specifically, the diagonal elements of the precision matrix are assigned a value of 0.33. If two nodes are adjacent, the corresponding element of the precision matrix is ​​assigned a value of -0.3, and all other locations are zero. The radius of the Earth is assumed to be 6371 kilometers.

[0140] To demonstrate the effectiveness of the method proposed in this paper, a high-speed parallel coordinate registration method for skywave over-the-horizon radar, assisted by a reference source, was compared with an estimation method using the classic Gibbs algorithm. The method was evaluated from three perspectives: the accuracy of ionospheric parameter estimation in the target area, the improvement in the target coordinate registration longitude, and the simulation calculation time. In this example, the proposed method, named Par-PGVI, serves as scenario one, and the classic Gibbs method used for comparison serves as scenario two.

[0141] like Figure 3 As shown, Figure 3 This figure shows the accuracy of ionospheric parameter estimation in the target region. The left axis shows the estimated IAP (ionospheric activity parameter) values ​​in all dimensions as the number of iterations increases. Note that the left axis has no units because it is an exponential. The right axis shows the change in the root mean square error (RMSE) as the number of iterations increases. The horizontal axis represents the number of iterations in logarithmic terms.

[0142] First, from Figure 3It can be seen that, first, the method of the present invention achieves an inference accuracy comparable to that of the control group and is slightly improved in the current experimental environment. That is, as the iterative process proceeds, the estimated value of IAPs rapidly increases from the prior value IAPs to prior =90 converges to the true value IAP real =120, and the root mean square error (RMSE) in all dimensions drops rapidly from 30 to 0.

[0143] Second, the inference process of Par-PGVI is more stable. Under the current experimental parameter settings, Figure 3 In the figure, the blue line is more tightly packed than the yellow line before step 100. After step 100, when the parameters converge near their true values, the blue line also exhibits significantly less oscillation than the yellow line. This allows us to obtain inference results with fewer iterations, without having to wait for the Markov chain to continue traversing to collect sufficient samples, as in the control experiment. This demonstrates the computational efficiency and real-time performance of our method.

[0144] Secondly, the improvement of the calibration target coordinate registration accuracy is analyzed. We placed some calibration targets within the current radar field of view. The results of the high-speed parallel coordinate registration method of the skywave over-the-horizon radar based on the reference source and the coordinate registration using the classic Gibbs algorithm are shown in the figure below. Figure 4 The details of each number are shown in Figure 5 and Figure 6 Table 1 lists the specific improvements in positioning error for each calibration target. Figure 5 、 Figure 6 As can be seen from Table 1, after using our algorithm to calculate the posterior estimate of the ionospheric parameters, the ionospheric electron concentration profile is calculated by substituting it into the international reference ionospheric model. Ray tracing is performed on the calibration target in this electron concentration profile to calculate the posterior position of the calibration target. Since our algorithm can calculate very accurate ionospheric activity parameters, accurate ionospheric electron concentration inversion can be obtained after substituting it into the ionospheric model. This makes the posterior position of all calibration targets very close to the actual position of the calibration target, and on average, a 73.60% registration error improvement rate is achieved. This experimental result shows that the algorithm proposed in the present invention has good registration accuracy in terms of skywave radar coordinate registration and target tracking, as well as excellent timeliness.

[0145] Table 1 Coordinate registration error and its improvement rate data

[0146]

[0147] Finally, simulation time. Table 2 presents the running time comparison of the two scenarios on the simulation platform. It can be seen that Par-PGVI (the method of the present application) has a significant advantage in computing speed. After only a few dozen iterations, it takes about half an hour to estimate IAPs in all dimensions to quickly approach the true value. Compared with the control group, when the same estimation accuracy is reached, Par-PGVI only needs 2% of the running time of the Gibbs sampling method. No matter from the number of iterations needed to reach the same estimation accuracy or the overall operation time, the new method has achieved a significant advantage over the control method. Analysis of the entire sampling process shows that under the current hyperparameter configuration, after only a few dozen iteration steps and about half a minute, the estimated values of ionospheric parameters in all dimensions can quickly approach the actual values. This phenomenon further highlights the effectiveness of the parallelization method proposed in the present application in solving the inference problems posed by high-dimensional complex computational models.

[0148] Table 2 RMSE and running time data of two algorithms

[0149]

[0150] Note: The running time of scenario one is very small when the number of iterations is 1, which can be ignored, so it is 0.

[0151] The above examples only express several embodiments of the present application, and the description is more specific and detailed, but it cannot be understood as a limitation on the scope of the patent of the present application. It should be noted that for ordinary skilled persons in the art, without departing from the concept of the present application, a number of modifications and improvements can be made, which are all within the scope of protection of the present application. Therefore, the scope of protection of the present application should be subject to the appended claims.

Claims

1. A high-speed parallel coordinate registration method for skywave over-the-horizon radar based on reference source assistance, characterized in that: The following steps are involved: The posterior distribution of ionospheric activity parameters is constructed and estimated based on the prior information of ionospheric activity parameters and reference source data in a Bayesian framework. According to the posterior distribution of ionospheric activity parameters, the estimated values ​​of ionospheric activity parameters in all dimensions are solved in parallel based on the proximal gradient descent method and the mean field variation method; Substituting the estimated values ​​of ionospheric activity parameters in all dimensions into the international reference ionospheric model, the a posteriori ionosphere is obtained. The three-dimensional ray tracing method is used to calculate the electromagnetic wave propagation path in the a posteriori ionosphere to obtain the coordinate registration result of the target to be located.

2. The method for high-speed parallel coordinate registration of skywave over-the-horizon radar based on reference source assistance according to claim 1 is characterized in that: According to the posterior distribution of the ionospheric activity parameters, the estimation of the ionospheric activity parameters of all dimensions is solved in parallel based on the proximal gradient descent method and the mean field variation method, including the following steps: Select the family of variational distributions of ionospheric activity parameters under the mean field assumption; Constructing the evidence lower bound function based on the posterior distribution and variational distribution family of ionospheric activity parameters; Solving the estimated values ​​of the ionospheric activity parameters in each dimension in parallel based on the evidence lower bound function; The estimated values ​​of the ionospheric activity parameters in each dimension are combined to obtain the estimated values ​​of the ionospheric activity parameters in all dimensions.

3. The method for high-speed parallel coordinate registration of skywave over-the-horizon radar based on reference source assistance according to claim 2, characterized in that: To obtain the estimated values ​​of ionospheric activity parameters for any dimension, the following steps are involved: Decomposing and recombining the evidence lower bound function into a function consisting of a first data item, a second data item, and a third data item; The first data item and the second data item are optimized using the proximal gradient descent method to obtain the optimal variational distribution of the ionospheric activity parameters in the current dimension; According to the optimal variational distribution, the ionospheric activity parameters are estimated to obtain the estimated values ​​of the ionospheric activity parameters in the current dimension; The first data item is a data item for calculating gradient information using only a numerical sampling method, the second data item is a data item for calculating gradient information using an analytical method, and the third data item is a data item that does not participate in the optimization process.

4. The method for high-speed parallel coordinate registration of skywave over-the-horizon radar based on reference source assistance according to claim 3 is characterized in that: Decomposing and recombining the evidence lower bound function into a function consisting of a first data item, a second data item, and a third data item, comprising the following steps: Performing mean field approximation decomposition on the evidence lower bound function; In the framework of proximal gradient optimization, all dimensions obtained by mean field approximation are reorganized to obtain a function consisting of a first data term, a second data term and a third data term.

5. The method for high-speed parallel coordinate registration of skywave over-the-horizon radar based on reference source assistance according to claim 3 is characterized in that: Optimizing the first data item and the second data item using a proximal gradient descent method includes the following steps: Calculating the gradient information of the first data item using the Monte Carlo method; Using the quasi-Monte Carlo method to control the variance of the gradient information estimation of the first data item; Establish a variational distribution parameter iteration formula within the framework of proximal gradient optimization; The variational distribution parameters are iteratively calculated until the convergence condition is met to obtain the optimal variational distribution of ionospheric activity parameters.

6. The method for high-speed parallel coordinate registration of skywave over-the-horizon radar based on reference source assistance according to claim 4, characterized in that: The function composed of the first data item, the second data item, and the third data item is: Where, is the first data item, is the second data item, else term is the third data item, x represents the ionospheric activity parameter vector, z represents the radar measurement vector of the reference source and the ground distance vector, x i represents the ionospheric activity parameter of the i-th dimension, x -i represents the ionospheric activity parameter vector of the remaining dimensions except the i-th dimension, λ represents the total variation distribution parameter, λ i Represents the variational parameter of the variational distribution of the current i-th dimension, q λ (x) is the variational distribution of all ionospheric activity parameters, is the variational distribution of the ionospheric activity parameter in the i-th dimension, log p(z|x) is the log-likelihood function, p(x i ,x -i ) is the prior distribution of all ionospheric activity parameters.

7. The method for high-speed parallel coordinate registration of skywave over-the-horizon radar based on reference source assistance according to claim 5, characterized in that: The variational distribution parameter iteration formula is: Where, It is a hyperparameter introduced in the proximal gradient descent method, which is used to characterize the strong constraint effect of the proximal term on the optimization process. is the damping parameter, represents the variational distribution of the ionospheric activity parameter of the i-th dimension at the k-th iteration, g(λ k ) represents the first data item of the k-th iteration; Among them, for a certain choice of variational distribution family, the following form of variational parameter iterative expression with damping form can be obtained: Where, represents the variational parameter of the variational distribution of the i-th dimension at the k+1-th iteration, represents the variational parameter of the variational distribution of the i-th dimension at the k-th iteration, and l represents the learning rate.

Images