A method for separating multiple echoes of SuperDARN radars based on high-resolution spectral analysis

Abstract

The application discloses a method for separating multiple echoes of SuperDARN radar based on high-resolution spectral analysis, which comprises the following steps: step 1, obtaining rawacf file data of SuperDARN radar original data, and pre-processing the rawacf file data; step 2, calculating the lag zero power of each distance gate in ACF data, comparing each integration time, and eliminating the distance gate with a signal-to-noise ratio less than zero; step 3, arranging fitacf file data obtained by fitting into a complete data set, and comparing with echo information; step 4, placing the remaining autocorrelation function matrix into a root-music and esprit algorithm model built, extracting the Doppler frequency corresponding to each distance gate, and obtaining Doppler velocity through the Doppler frequency; and step 5, constructing power spectrum estimation of radar echo through an ARMA model, obtaining power spectrum and spectral width, and writing the obtained spectral width into a corresponding fitacf file. The application improves the resolution capability of the radar by using a high-resolution spectral analysis method.

Description

Technical Field

[0001] The present invention relates to the technical fields of ionospheric physics and atmospheric science, and in particular to a method for separating multiple echoes of a SuperDARN radar based on high-resolution spectrum analysis. Background Art

[0002] The Super Dual Auroral Radar Network (SuperDARN) is a global network of over 30 high-frequency radars primarily used to monitor plasma circulation in high latitudes in the Northern and Southern Hemispheres. These radars exploit the ionosphere's refraction of high-frequency radio waves, causing them to intersect perpendicularly with Earth's magnetic field at high latitudes. This property enables the radars to measure irregularities in the ionosphere's electron density aligned with the magnetic field, recording variations in the Doppler shift of their backscattered radiation, typically in the 8 to 20 MHz frequency range. These irregularities are considered tracers of plasma circulation.

[0003] Currently, Super DARN data are distributed to users in a range-time format representing complex autocorrelation functions. Users can process these ACF data using the Radar Software Toolkit (RST) or custom software. RST is the standard data analysis tool for the Super DARN, containing all the routines necessary to convert ACF datasets from multiple radars into high-latitude ionospheric potential maps. With the rapid development of spatial spectrum estimation, high-resolution spectrum estimation, represented by the MUSIC and ESPRITS algorithms, has developed rapidly. As a major branch of spatial spectrum estimation, it has a wide range of applications and demonstrates superior performance in many fields. For radar-based space physics research, such as Super DARN data, high-resolution spectral analysis methods can significantly improve data analysis accuracy. These algorithms can be applied to process radar echo data to improve the accuracy of Doppler velocity, spectrum width, and other related parameters, thereby obtaining a more detailed ionospheric potential distribution map. However, in the standard process of SuperDARN analyzing its frequency characteristics to obtain its signal characteristics, there is a problem of being unable to extract ionospheric echoes that are distinguished as ground scattered echoes due to their too small spectrum width. This causes ground scattering to interfere with ionospheric data collection, and the accuracy and reliability of radar measurements are not high. Summary of the Invention

[0004] The object of the present invention is to provide a method for separating multiple echoes of SuperDARN radar based on high-resolution spectrum analysis. In order to make up for the shortcomings of the existing method of extracting echo parameters from raw data, the present invention provides the following technical solution: A method for separating multiple echoes of SuperDARN radar based on high-resolution spectrum analysis, comprising:

[0005] Step 1: Obtain the rawacf file data of the SuperDARN radar raw data and preprocess its autocorrelation function, wherein the preprocessing includes constructing a cross product of the ACF data of the rawacf file data and performing phase unwrapping on the phase with a lag period exceeding 2π;

[0006] Preferably, the cross product R(τ) is constructed for the ACF data in the rawacf file data, R(τ)=|R(τ)|e jφ(τ) , where e j It is the exponential representation of a complex number, φ(τ) is expressed as the phase, and for the phase φ(τ) that may cause phase wrapping, phase unwrapping is performed on the phase with a lag period exceeding 2π.

[0007] Step 2: Calculate the hysteresis zero power of each range gate in the ACF data, and eliminate the range gates with a signal-to-noise ratio less than zero by comparing each integration time;

[0008] Preferably, the process of removing range gates with a signal-to-noise ratio less than zero includes:

[0009] Step 201: Based on the constructed cross product R(τ)=|R(τ)|e jφ(τ) The specific expansion calculation formula is as follows:

[0010]

[0011] in, is the complex conjugate of the signal s(t), t is time, τ is a fixed time lag, E is the mean, e -jωτ It is the exponential form of a complex number. Since the real and imaginary parts of the ACF data at each lag period are the cross product R(τ) of the response of a sequence on multiple pulse sequences, one of the sequences is a fixed time lag τ.

[0012] Step 202: Calculate the signal-to-noise ratio Q in decibels based on the delayed zero power R(0) according to the following formula. The calculation formula of the signal-to-noise ratio Q is as follows:

[0013]

[0014] Where R(0) is the lagging zero power, R n is the noise power estimate calculated from the ten lowest lag power values ​​for each integration time, where the integration time is the beam, p is the modulus of the ACF data, w is the probability density function, and d is the sign of the integration;

[0015] Step 203: Further improve the noise level by adding a correction factor μ, that is, when receiving the signal, calculate the zero-lag power Rcur(0) of the current ACF data and the zero-lag power Rint(0) of the ACF data from the potential interference range, and further compare the power of the two The ratio μ is calculated. When the ratio μ is greater than or equal to 0.3 Navg, the corresponding interference hysteresis is marked as "bad" and removed from the rawacf file data.

[0016] Step 3: The rawacf data are fitted with the least squares method of the rst software package to obtain the fitacf file data, which is then organized into a complete data set and compared with the echo information obtained from high-resolution spectral analysis.

[0017] Step 4: Put the remaining autocorrelation function matrix into the built root-music algorithm and esprit algorithm model, extract the Doppler frequency corresponding to each range gate, obtain the Doppler velocity through the Doppler frequency, and write the obtained Doppler velocity into the fitacf file;

[0018] Preferably, the root-music algorithm process includes:

[0019] a1: By analyzing a large amount of simulated ACF data, the speed and spectrum width of different delay numbers N are systematically calculated;

[0020] a2: Select the dimension p of the autocorrelation matrix. p = 8 is chosen based on the good compromise between the maximum number of signal sources expected in the signal and the number of points in the ACF data. The preprocessed ACF matrix is ​​used to construct the Hermite covariance matrix R p , the specific formula is as follows:

[0021]

[0022] Among them, Y n is the ACF matrix, is the complex conjugate of the ACF matrix, N is the number of time delays, and p is the dimension of the autocorrelation matrix;

[0023] a3: R p Perform eigenvalue decomposition R p =EΛE H , use the distribution of eigenvalues ​​to get the number of information sources, find the eigenvector corresponding to the noise eigenvalue through the number of information sources, and construct the noise subspace U N , construct the polynomial f(z) based on the eigenvector of the noise subspace, and the calculation formula of the polynomial f(z) is as follows:

[0024]

[0025] Among them, p(z) is a vector related to frequency f, z=e j2πf , M is the order of the signal, and H is the conjugate transpose.

[0026] Preferably, the process of extracting the Doppler frequency corresponding to each range gate and obtaining the Doppler velocity through the Doppler frequency includes:

[0027] Solve the root zk of the polynomial f(z) and get M roots zk. By focusing only on the roots on the unit circle or the MN roots zk closest to the unit circle, the angular frequency ω corresponding to the root zk of the polynomial is arg(zk), and the angular frequency is converted into the actual frequency Among them, f s is the sampling frequency of the superdarn radar. When the Doppler frequency is obtained, it is converted into the Doppler velocity v D , v D The calculation formula is as follows:

[0028]

[0029] Where c is the speed of light, f D is the Doppler frequency and f0 is the transmitting frequency of the radar.

[0030] Preferably, the esprit algorithm model is mainly to construct the signal subspace U S The specific process includes:

[0031] b1: The preprocessed Hermite covariance matrix R p Put it into the esprit algorithm model to extract the signal subspace U S ;

[0032] b2: Use the linear array structure to transform U S Divide matrix into subarrays and That is, using the shift invariance assumption, there is a rotation invariance relationship between the two subarrays: Where Φ is the rotation matrix, which contains the frequency or DOA information of the signal;

[0033] b3: Use the least squares method to minimize the following objective function to solve the rotation matrix Φ. The specific formula for minimizing the objective function using the least squares method is as follows:

[0034]

[0035] where ||·|| F is the Frobenius norm, which is the square root of the sum of the squares of the elements of the matrix;

[0036] b4: After obtaining the rotation matrix Φ, Φ=VΛV through eigenvalue decomposition -1 , where V is a diagonal matrix, V -1 It is the inverse matrix of V, Λ is a diagonal matrix, and the eigenvalue of Λ contains the frequency or direction information of the signal;

[0037] b5: Convert the angular frequency ω = arg(zk) corresponding to the eigenvalue zk of Λ into the actual frequency After obtaining the Doppler frequency, it is converted into Doppler velocity.

[0038] Step 5: Use the ARMA model to construct the power spectrum estimate of the radar echo, obtain the power spectrum for each range gate, use the full width at half maximum to obtain the spectrum width, and write the obtained spectrum width into the corresponding fitacf file;

[0039] Preferably, the calculation process of obtaining the spectrum width using the full width at half maximum includes:

[0040] s1: Construct ARMA model x(t). The specific calculation formula is as follows:

[0041]

[0042] Among them, a i is the polynomial coefficient of the AR model, i is the order of the AR model, b j is the polynomial coefficient b of the MA model i , j is the order of the MA model;

[0043] s2: By determining the order of the AR model and the MA model, we can determine the order corresponding to the modulus of the ACF data approaching 0. For the AR model, we can solve the Yule-Walker equation to obtain the polynomial coefficients ai of the AR model. After obtaining the polynomial of the AR model, we can further obtain the polynomial coefficients b of the MA model. i , polynomial coefficient b i The calculation formula is as follows:

[0044]

[0045] Where R(0) is the value of the autocorrelation function at zero lag, and R(j) is the value of the autocorrelation function at time lag j.

[0046] S3: When the coefficient a of the AR model is obtained i and the coefficient b of the MA model j Finally, the modern spectrum P(f) of the ARMA model is constructed. The calculation formula of the modern spectrum P(f) is as follows:

[0047]

[0048] in, Expressed as the variance of white noise, For the variance of white noise, a0 = 1. Therefore, given the autocorrelations r(0), r(1), ..., r(M), determine the variance of white noise

[0049] S4: Therefore, the spectrum width is defined using the full width at half maximum (FWHM) of the power spectrum, and the frequency range where the power spectrum reaches half of its maximum value is calculated. The corresponding delay of the calculated ACF data is then compared and verified with the spectrum width and Doppler velocity obtained by the standard process.

[0050] Step 6: Use the fanplot tool in the pydarn package to visualize the data. Given specific parameters, including Doppler velocity and spectrum width, project specific or all beams onto a geographic map in geographic coordinates. The fanplot distinguishes ionospheric scattered signals from ground scattered signals using the criteria of v < 30 m / s and w < 35 m / s.

[0051] Step 7: Use three evaluation indicators, namely root mean square error (RMSE), normalized cross correlation (NCC), and Pearson correlation coefficient (PPMCC), to evaluate the data effect and output the evaluation results.

[0052] Preferably, the calculation formulas of the root square error RMSE, normalized cross correlation NCC, and Pearson correlation coefficient PPMCC are as follows:

[0053]

[0054] Where m is the number of pixels in image a, n is the number of pixels in image b, I(i,j) represents the pixel value of image X at the grid point position, K(i,j) represents the pixel value of image Y at the grid point position, and x i is the pixel value or data value of image a in row i and column j, μ x is the pixel mean of image a, y i is the pixel value or data value of image b in row i and column j, μ y is the pixel mean of image b, ρ X,Y is the Pearson coefficient, σ X and σ Y are the standard deviations of X and Y, cov(X,Y) is the covariance of X and Y, μ X is the mean of X, μ Y is the mean of Y, X and Y are the datasets of two images respectively.

[0055] Compared with the existing technology, the beneficial effects achieved by the present invention are: based on spatial spectrum estimation and modern spectrum estimation methods, the present invention analyzes its frequency characteristics through high-resolution spectrum analysis methods and establishes a statistical model of the signal to obtain its signal characteristics. Compared with the standard process of SuperDARN, it can extract ionospheric echoes that are distinguished as ground scattered echoes due to their too small spectrum width from ground scattered echoes, and expand the spatial coverage range of ionospheric sources to the area where the standard method can only identify ground echoes. It pioneered the use of array signal processing technology to extract SuperDARN radar to distinguish different target echoes, and pioneered the application of array signal processing technology and modern spectrum estimation technology to identify signal sources, reducing the interference of ground scattering on ionospheric data acquisition, and improving the accuracy and reliability of radar measurements. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] The accompanying drawings are used to provide a further understanding of the present invention and constitute a part of the specification. Together with the embodiments of the present invention, they are used to explain the present invention and do not constitute a limitation of the present invention. In the accompanying drawings:

[0057] Figure 1 A schematic flow chart of the method steps for separating SuperDARN radar multiple echoes based on high-resolution spectrum analysis provided in an embodiment of the present invention;

[0058] Figure 2 A schematic diagram of an operating model provided by an embodiment of the present invention;

[0059] Figure 3 This is a schematic diagram of the data comparison results provided by an embodiment of the present invention. DETAILED DESCRIPTION

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

[0061] See also Figures 1 to 3 In order to make up for the shortcomings of the existing methods for extracting echo parameters from raw data, the present invention provides a technical solution: a method for separating multiple echoes of SuperDARN radar based on high-resolution spectrum analysis, comprising the following steps:

[0062] Step 1: Obtain the rawacf file data of the SuperDARN radar raw data and preprocess its autocorrelation function;

[0063] In this embodiment, the preprocessing process of preprocessing the autocorrelation function includes: for the ACF data of rawacf, constructing the cross product R(τ)=|R(τ)|e jφ(τ) , e j It is an exponential representation of a complex number, φ(τ) represents the phase, and based on the phase φ(τ) that causes phase wrapping, phase unwrapping is performed on the phase with a lag period exceeding 2π.

[0064] Step 2: Calculate the hysteresis zero power of each range gate and compare each integration time (beam), and eliminate the range gates with a signal-to-noise ratio Q less than zero;

[0065] In this embodiment, the real and imaginary parts of the ACF at each lag period are obtained by taking the cross product R(τ) of the response of a sequence on multiple pulse sequences, where one sequence is a fixed lag τ, i.e., R(τ) = |R(τ)|e constructed in step 1. jφ(τ) The specific expansion calculation formula is as follows:

[0066]

[0067] in, is the complex conjugate of the signal s(t), t is time, τ is a fixed time lag, E is the mean, e -jωτ is the exponential form of a complex number,

[0068] Exemplarily, the signal-to-noise ratio Q in decibels is calculated from the lagging zero power R(0) based on the following formula, where R n is the noise power estimate calculated from the ten lowest lag power values ​​for each integration time (beam), p is the modulus of the ACF, w is the probability density function, and d is the sign of the integration.

[0069]

[0070] Specifically, using only the ten lowest hysteresis power values ​​as the upper limit of integration cannot reflect the entire noise level. This reduction in the upper limit of integration will cause a deviation in the lower values ​​in the average estimate, and the estimated value of the noise power will be lower than the actual noise. Therefore, a correction factor μ needs to be added to improve the noise level.

[0071] Furthermore, when receiving signals, the measurement results may be inaccurate due to the overlap or interference of signals in multiple distance units. Therefore, it is necessary to process the span interference. First, it is necessary to calculate the zero-lag power Rcur(0) of the current ACF and the zero-lag power Rint(0) of the ACF from the potential interference range, and then compare the power of the two. Calculate the ratio μ. When the ratio μ is greater than or equal to 0.3 Navg, mark the corresponding interference hysteresis as "bad" and remove it from the rawacf file.

[0072] Step 3: The rawacf data were fitted with the least squares method of the rst software package to obtain the fitacf file data, which was then organized into a complete data set and further compared with the echo information obtained from high-resolution spectral analysis.

[0073] Step 4: Place the remaining autocorrelation function matrix into the established root-music and esprit algorithm models, extract the Doppler frequency corresponding to each range gate, obtain the Doppler velocity through the Doppler frequency, and write the obtained Doppler velocity into the fitacf file.

[0074] In this embodiment, the root-music algorithm process is as follows:

[0075] By first analyzing a large number of simulated ACFs, the speed and spectrum width are systematically calculated for different delay numbers; the total number of points in the ACF is 23, of which 15% are missing values ​​called "bad lag points", which also appear in real data due to reception interruptions or other disturbances during transmission. The dimension p of the autocorrelation matrix is ​​chosen to be p = 8, which represents a good compromise between the maximum number of signal sources expected in the signal and the number of points in the ACF. The preprocessed ACF matrix is ​​used to construct the Hermite covariance matrix R by the formula p , the specific formula is as follows:

[0076]

[0077] Among them, Y n is the ACF matrix, is the complex conjugate of the ACF, N is the number of time delays, p is the dimension of the autocorrelation matrix,

[0078] Furthermore, R p Perform eigenvalue decomposition R p =EΛE H , we use the distribution of eigenvalues ​​to obtain the number of information sources, find the eigenvector corresponding to the noise eigenvalue through the number of information sources, and construct the noise subspace U N , construct the polynomial f(z) based on the eigenvector of the noise subspace, and the formula of the polynomial f(z) is as follows:

[0079]

[0080] Among them, p(z) is a vector related to frequency f, z=e j2πf , M is the order of the signal, and H is the conjugate transpose.

[0081] In this embodiment, the process of obtaining the Doppler frequency is as follows:

[0082] Solve the roots of the polynomial f(z). Based on the M roots obtained, by only focusing on the roots on the unit circle or the MN roots zk closest to the unit circle, the angular frequency ω corresponding to the root zk of the polynomial is arg(zk), and then convert the angular frequency into the actual frequency Among them, f s is the sampling frequency of the superdarn radar. Once the Doppler frequency is obtained, it can be converted into the Doppler velocity v D , v D The calculation formula is as follows:

[0083]

[0084] Where c is the speed of light, f D is the Doppler frequency, f0 is the radar's transmitting frequency, based on which the Doppler velocity of the target echo is obtained.

[0085] Furthermore, by p Put it into the esprit algorithm model. Different from the music algorithm, esprit uses the construction of signal subspace. Its characteristics are small amount of calculation and no need for spectrum peak search. p After eigenvalue decomposition, the signal subspace U is extracted S , using the linear array structure to S Divide matrix into subarrays and Using the shift invariance assumption, there is a rotation invariance relationship between the two subarrays: Where Φ is the rotation matrix, which contains the frequency or DOA information of the signal. The following objective function can be minimized by the least squares method to solve the rotation matrix Φ. The specific formula for minimizing the objective function by the least squares method is as follows:

[0086]

[0087] where ||·|| F is the Frobenius norm, that is, the square root of the sum of the squares of the elements of the matrix. After we get the rotation matrix Φ, we can use the eigenvalue decomposition Φ=VΛV -1 , where V is a diagonal matrix, V -1 It is the inverse matrix of V, Λ is a diagonal matrix, the eigenvalue of Λ contains the frequency or direction information of the signal. Like music, the angular frequency ω=arg(zk) corresponding to the eigenvalue zk of Λ is converted into the actual frequency Once the Doppler frequency is obtained, it is converted into Doppler velocity.

[0088] Step 5: Use the ARMA model to construct the power spectrum estimate of the radar echo, obtain the power spectrum for each range gate, use the full width at half maximum to obtain the spectrum width, and write the obtained spectrum width into the corresponding fitacf file;

[0089] In this embodiment, the specific process of spectrum width calculation is as follows:

[0090] First, construct the ARMA model x(t):

[0091]

[0092] Among them, a i is the polynomial coefficient of the AR model, i is the order of the AR model, b j is the polynomial coefficient b of the MA model i , j is the order of the MA model;

[0093] Specifically, by determining the order of the AR model and the MA model, the order corresponding to the modulus of the acf approaching 0 is determined. For the AR model, the polynomial coefficient a of the AR model is obtained by solving the Yule-Walker equation. i After obtaining the polynomial of the AR model, we also need to find the polynomial coefficient b of the MA model. i , the formula is as follows:

[0094]

[0095] Where R(0) is the value of the autocorrelation function at zero lag, and R(j) is the value of the autocorrelation function at time lag j.

[0096] Furthermore, when the coefficient a of the AR model is obtained i and the coefficient b of the MA model j Finally, the modern spectrum P(f) of the ARMA model is constructed, and the formula is as follows:

[0097]

[0098] in, Expressed as the variance of white noise, For the variance of white noise, a0 = 1. Therefore, given the autocorrelations r(0), r(1), ..., r(M), the variance of white noise can be determined

[0099] Exemplarily, the spectrum width is defined by the full width at half maximum (FWHM) of the power spectrum, and the frequency range when the power spectrum reaches half of its maximum value is calculated; the raw ACF data used in the present invention is derived from SuperDARN radar detection data. The SuperDARN radar transmits a pulse sequence so that the time delay between any two pulses in the sequence is different, and a regular sequence is formed, thereby calculating the corresponding delay of the ACF, and performing data comparison and verification with the spectrum width and Doppler velocity obtained by the SuperDARN standard process.

[0100] Step 6: Use the fanplot tool in the pydarn package to visualize the data. Given specific parameters, including Doppler velocity and spectrum width, project specific or all beams onto a geographic map in geographic coordinates. The fanplot distinguishes ionospheric scattered signals from ground scattered signals using the criteria of v < 30 m / s and w < 35 m / s.

[0101] Step 7: Use the root mean square error (RMSE), normalized cross correlation (NCC), and Pearson correlation coefficient (PPMCC) to evaluate the data effect. The formulas for RMSE, NCC, and PPMCC are as follows:

[0102]

[0103] Where m is the number of pixels in image a, n is the number of pixels in image b, I(i,j) represents the pixel value of image X at the grid point position, K(i,j) represents the pixel value of image Y at the grid point position, and x i is the pixel value or data value of image a in row i and column j, μ x is the pixel mean of image a, y i is the pixel value or data value of image b in row i and column j, μ y is the pixel mean of image b, ρ X,Y is the Pearson coefficient, σ X and σ Y are the standard deviations of X and Y, cov(X,Y) is the covariance of X and Y, μ X is the mean of X, μ Y is the mean of Y, X and Y are the datasets of two images respectively.

[0104] The present invention is based on the development of spatial spectrum estimation and array signal technology, adopts high-resolution spectral analysis method, and uses the massive amount of information-rich data obtained by SuperDARN radar. The relationship between the data is expressed by high-resolution spectral analysis with powerful fitting ability, and the SuperDARN radar is used to distinguish the echoes of different targets, thereby improving the radar's resolution ability.

[0105] It should be noted that, in this document, relational terms such as first and second, etc., are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "comprises," "comprising," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that includes a list of elements includes not only those elements but also other elements not explicitly listed, or elements inherent to such process, method, article, or apparatus.

[0106] Finally, it should be noted that the above descriptions are merely preferred embodiments of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art will be able to modify the technical solutions described in the aforementioned embodiments or substitute equivalents for some of the technical features. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention shall be included within the scope of protection of the present invention.

Claims

1. A method for separating multiple echoes of SuperDARN radar based on high-resolution spectrum analysis, characterized by: The following steps are involved: Step 1: Obtain rawacf file data of SuperDARN radar raw data and preprocess its autocorrelation function, wherein the preprocessing includes constructing a cross product of the ACF data of the rawacf file data and performing phase unwrapping on the phase with a lag period exceeding 2π; Step 2: Calculate the hysteresis zero power of each range gate in the ACF data, and eliminate the range gates with a signal-to-noise ratio less than zero by comparing each integration time; Step 3: The rawacf data are fitted with the least squares method of the rst software package to obtain the fitacf file data, which is then organized into a complete data set and compared with the echo information obtained from high-resolution spectral analysis. Step 4: Put the remaining autocorrelation function matrix into the built root-music algorithm and esprit algorithm model, extract the Doppler frequency corresponding to each range gate, obtain the Doppler velocity through the Doppler frequency, and write the obtained Doppler velocity into the fitacf file; Step 5: Use the ARMA model to construct the power spectrum estimate of the radar echo, obtain the power spectrum for each range gate, use the full width at half maximum to obtain the spectrum width, and write the obtained spectrum width into the corresponding fitacf file; Step 6: Use the fan plot in the pydarn package to visualize the data. Given specific parameters, including Doppler velocity and spectrum width, project specific or all beams onto a geographic map in geographic coordinates. The fan plot distinguishes ionospheric scattered signals from ground scattered signals using the criteria of v < 30 m / s and w < 35 m / s. Step 7: Use three evaluation indicators, namely root mean square error (RMSE), normalized cross correlation (NCC), and Pearson correlation coefficient (PPMCC), to evaluate the data effect and output the evaluation results.

2. The method for separating SuperDARN radar multiple echoes based on high-resolution spectrum analysis according to claim 1, characterized in that: The pre-processing of the autocorrelation function further includes: constructing a cross product R(τ) for the ACF data in the rawacf file data, R(τ)=|R(τ)|e jφ(τ) , where e j It is the exponential representation of a complex number, φ(τ) is expressed as the phase, and for the phase φ(τ) that may cause phase wrapping, phase unwrapping is performed on the phase with a lag period exceeding 2π.

3. The method for separating SuperDARN radar multiple echoes based on high-resolution spectrum analysis according to claim 2, characterized in that: The process of eliminating range gates with a signal-to-noise ratio less than zero comprises the following steps: Step 201: Based on the constructed cross product R(τ)=|R(τ)|e jφ(τ) The specific expansion calculation formula is as follows: in, is the complex conjugate of the signal s(t), t is time, τ is a fixed time lag, E is the mean, e -jωτ It is the exponential form of a complex number. Since the real and imaginary parts of the ACF data at each lag period are the cross product R(τ) of the response of a sequence on multiple pulse sequences, one of the sequences is a fixed time lag τ. Step 202: Calculate the signal-to-noise ratio Q in decibels based on the delayed zero power R(0) according to the following formula. The calculation formula of the signal-to-noise ratio Q is as follows: Where R(0) is the lagging zero power, R n is the noise power estimate calculated from the ten lowest lag power values ​​for each integration time, where the integration time is the beam, p is the modulus of the ACF data, w is the probability density function, and d is the sign of the integration; Step 203: Improve the noise level by adding a correction factor μ, that is, when receiving the signal, calculate the zero-lag power Rcur(0) of the current ACF data and the zero-lag power Rint(0) of the ACF data from the potential interference range, and compare the two powers. Calculate the ratio μ. When the ratio μ is greater than or equal to 0.3 Navg, mark the corresponding interference hysteresis as "bad" and remove it from the rawacf file data.

4. The method for separating SuperDARN radar multiple echoes based on high-resolution spectrum analysis according to claim 1, characterized in that: The root-music algorithm process includes: a1: By analyzing a large amount of simulated ACF data, the speed and spectrum width of different delay numbers N are systematically calculated; a2: Select the dimension p of the autocorrelation matrix. p = 8 is chosen based on the good compromise between the maximum number of signal sources expected in the signal and the number of points in the ACF data. The preprocessed ACF matrix is ​​used to construct the Hermite covariance matrix R p , the specific formula is as follows: Among them, Y n is the ACF matrix, is the complex conjugate of the ACF matrix, N is the number of time delays, and p is the dimension of the autocorrelation matrix; a3: R p Perform eigenvalue decomposition R p =EΛE H , use the distribution of eigenvalues ​​to get the number of information sources, find the eigenvector corresponding to the noise eigenvalue through the number of information sources, and construct the noise subspace U N , construct the polynomial f(z) based on the eigenvector of the noise subspace, and the calculation formula of the polynomial f(z) is as follows: Among them, p(z) is a vector related to frequency f, z=e j2πf , M is the order of the signal, and H is the conjugate transpose.

5. The method for separating SuperDARN radar multiple echoes based on high-resolution spectrum analysis according to claim 4, characterized in that: The process of extracting the Doppler frequency corresponding to each range gate and obtaining the Doppler velocity through the Doppler frequency includes: Solve the root zk of the polynomial f(z) and get M roots zk. By focusing only on the roots on the unit circle or the MN roots zk closest to the unit circle, the angular frequency ω corresponding to the root zk of the polynomial is arg(zk), and the angular frequency is converted into the actual frequency Among them, f s is the sampling frequency of the superdarn radar. When the Doppler frequency is obtained, it is converted into the Doppler velocity v D , v D The calculation formula is as follows: Where c is the speed of light, f D is the Doppler frequency and f0 is the transmitting frequency of the radar.

6. The method for separating SuperDARN radar multiple echoes based on high-resolution spectrum analysis according to claim 5, characterized in that: The esprit algorithm model is mainly to construct the signal subspace U S The specific process includes: b1: The preprocessed Hermite covariance matrix R p Put it into the esprit algorithm model to extract the signal subspace U S ; b2: Use the linear array structure to transform U S Divide matrix into subarrays and That is, using the shift invariance assumption, there is a rotation invariance relationship between the two subarrays: Where Φ is the rotation matrix, which contains the frequency or DOA information of the signal; b3: Use the least squares method to minimize the following objective function to solve the rotation matrix Φ. The specific formula for minimizing the objective function using the least squares method is as follows: where ||·|| F is the Frobenius norm, which is the square root of the sum of the squares of the elements of the matrix; b4: After obtaining the rotation matrix Φ, Φ=VΛV through eigenvalue decomposition -1 , where V is a diagonal matrix, V -1 It is the inverse matrix of V, Λ is a diagonal matrix, and the eigenvalue of Λ contains the frequency or direction information of the signal; b5: Convert the angular frequency ω = arg(zk) corresponding to the eigenvalue zk of Λ into the actual frequency After obtaining the Doppler frequency, it is converted into Doppler velocity.

7. The method for separating SuperDARN radar multiple echoes based on high-resolution spectrum analysis according to claim 1, characterized in that: The calculation process of obtaining the spectrum width using the full width at half maximum includes: s1: Construct ARMA model x(t). The specific calculation formula is as follows: Among them, a i is the polynomial coefficient of the AR model, p is the order of the AR model, b j are the polynomial coefficients of the MA model, and q is the order of the MA model; s2: By determining the order of the AR model and the MA model, we can determine the order at which the modulus of the ACF data approaches 0. For the AR model, we can obtain the polynomial coefficient a of the AR model by solving the Yule-Walker equation. i , after obtaining the polynomial of the AR model, the polynomial coefficient b j The calculation formula is as follows: Where R(0) is the value of the autocorrelation function at zero lag, and R(j) is the value of the autocorrelation function at time lag j. S3: When the coefficient a of the AR model is obtained i and the coefficient b of the MA model j Finally, the modern spectrum P(f) of the ARMA model is constructed. The calculation formula of the modern spectrum P(f) is as follows: in, Expressed as the variance of white noise, For the variance of white noise, a0=1. Therefore, given the autocorrelations r(0), r(1),…, r(M), we can determine the variance of white noise. S4: Therefore, the spectrum width is defined by the full width at half maximum through the power spectrum, and the frequency range when the power spectrum reaches half of its maximum value is calculated. The corresponding delay of the calculated ACF data is then compared and verified with the spectrum width and Doppler velocity obtained by the standard process.

8. The method for separating SuperDARN radar multiple echoes based on high-resolution spectrum analysis according to claim 1, characterized in that: The calculation formulas of the root square error RMSE, normalized cross correlation NCC, and Pearson correlation coefficient PPMCC are as follows: Where m is the number of pixels in image a, n is the number of pixels in image b, I(i,j) represents the pixel value of image X at the grid point position, K(i,j) represents the pixel value of image Y at the grid point position, and x i is the pixel value or data value of image a in row i and column j, μ x is the pixel mean of image a, y i is the pixel value or data value of image b in row i and column j, μ y is the pixel mean of image b, ρ X,Y is the Pearson coefficient, σ X and σ Y are the standard deviations of X and Y, cov(X,Y) is the covariance of X and Y, μ X is the mean of X, μ Y is the mean of Y, X and Y are the datasets of two images respectively.

Images