Method for separating SuperDARN radar multiple echoes based on high-resolution spectral analysis

Abstract

The invention discloses a method for separating multiple echoes of a SuperDARN radar based on high-resolution spectral analysis, and the method comprises the following steps: 1, obtaining the rawacf file data of the original data of the SuperDARN radar, and carrying out the preprocessing of the rawacf file data; 2, calculating the lag zero power of each range gate in the ACF data, comparing each integration time, and removing the range gates with the signal-to-noise ratio smaller than zero; 3, arranging the fitf file data obtained by fitting into a complete data set, and comparing the complete data set with echo information; step 4, putting the residual autocorrelation function matrix into the constructed root-music and esprit algorithm model, extracting Doppler frequency corresponding to each range gate, and obtaining Doppler velocity through the Doppler frequency; step 5, constructing power spectrum estimation of radar echoes through an ARMA model to obtain a power spectrum and a spectrum width; and writing the obtained spectrum width into a corresponding fitf file. According to the invention, the high-resolution spectrum analysis method is adopted to improve the resolution capability of the radar.

Description

Technical Field

[0001] The 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 more than 30 high-frequency radars that are primarily used to monitor plasma circulation in high-latitude regions in the northern and southern hemispheres. These radars use the ionosphere's refraction properties of high-frequency radio waves, causing the waves to intersect the Earth's magnetic field perpendicularly at high latitudes. This property allows the radar to measure irregular areas of electron density in the ionosphere that are aligned along the magnetic field and record changes in their backscattered Doppler shifts, typically in the 8 to 20 MHz frequency range. These irregular areas are considered tracers of plasma circulation.

[0003] Currently, the data of the Super Auroral Radar Network are distributed to users in the range-time format of the complex "autocorrelation" function. Users can process these ACF data through the Radar Software Toolkit (RST) or custom software. RST is the standard data analysis tool of the Super Auroral Radar Network, which contains all the routines required to convert the ACF data sets of multiple radars into high-latitude ionospheric potential maps. With the rapid development of spatial spectrum estimation, high-resolution spectrum estimation represented by the music algorithm and the esprits algorithm has developed very rapidly. As a major branch of spatial spectrum estimation, it has a wide range of applications and has shown very superior performance in many fields. For radar-based space physics research, such as the Super Auroral Radar Network (SuperDARN) data, high-resolution spectral analysis methods can greatly improve the analysis accuracy of data. 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, it is impossible to extract ionospheric echoes that are resolved as ground scattered echoes due to their too small spectrum width, causing 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 deficiencies of the existing method for extracting echo parameters from raw data, the present invention provides the following technical solutions: A method for separating multiple echoes of SuperDARN radar based on high-resolution spectrum analysis, comprising:

[0005] Step 1: Obtain rawacf file data of the raw data of the SuperDARN radar, and preprocess its autocorrelation function, wherein the preprocessing includes constructing a cross product for the ACF data of the rawacf file data, and performing phase unwrapping for 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 represented as the phase, and for the phase φ(τ) that may cause phase wrapping, phase unwrapping is performed on the phase whose lag period exceeds 2π.

[0007] Step 2: Calculate the hysteresis zero power of each range gate in the ACF data, and remove 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 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 according to the delayed zero power R(0) based on 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 increase the noise level by adding a correction factor μ, that is, when receiving the signal, by calculating 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 comparing the power of the two The ratio μ is calculated, and when the ratio μ is greater than or equal to 0.3Navg, the corresponding interference hysteresis is marked as "bad" and removed from the rawacf file data.

[0016] Step 3: The fitacf file data obtained by least squares fitting of the rawacf data using the rst software package is organized into a complete data set and compared with the echo information obtained by 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: The dimension p of the autocorrelation matrix is ​​selected. p = 8 is selected based on a 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 , a polynomial f(z) is constructed based on the eigenvector of the noise subspace. 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: U S =U S2 Φ, 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 the least squares method to minimize the objective function is as follows:

[0034]

[0035] Among them ∥·∥ 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 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 Once the Doppler frequency is obtained, it is converted into Doppler velocity.

[0038] Step 5: Use the ARMA model to construct the power spectrum estimation of the radar echo, obtain the power spectrum of 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 are the polynomial coefficients of the AR model, i is the order of the AR model, and 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(i) is the value of the autocorrelation function at time lag i;

[0046] S3: When the coefficients ai of the AR model and the coefficients b of the MA model are obtained i After that, 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, where 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 by the power spectrum using the full width at half maximum (FWHM), 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.

[0050] Step 6: Use the fanplot 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 by 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 mean square error RMSE, the normalized cross correlation NCC, and the 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, K(i,j) represents the pixel value of image Y at the grid point, 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 the image b in the i-th row and j-th column, μ 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, respectively, 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 prior art, the beneficial effects achieved by the present invention are as follows: based on the methods of spatial spectrum estimation and modern spectrum estimation, the present invention analyzes the frequency characteristics of the signal by means of a high-resolution spectrum analysis method and by establishing a statistical model of the signal to obtain its signal characteristics; compared with the standard process of SuperDARN, the present invention can extract ionospheric echoes that are distinguished as ground scattering echoes due to their too small spectrum width from ground scattering echoes, and expand the spatial coverage of ionospheric sources to the area where the standard method can only identify ground echoes; the array signal processing technology is pioneeringly used to extract the SuperDARN radar to distinguish different target echoes; the array signal processing technology and the modern spectrum estimation technology are pioneeringly applied to identify signal sources, reduce the interference of ground scattering on ionospheric data acquisition, and improve the accuracy and reliability of radar measurement. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0057] In the attached picture:

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

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

[0060] Figure 3 A schematic diagram of data comparison results provided by an embodiment of the present invention. DETAILED DESCRIPTION

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

[0062] See also Figures 1 to 3 In order to make up for the shortcomings of the existing method of 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:

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

[0064] 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 will cause phase wrapping, phase unwrapping is performed on the phase whose lag period exceeds 2π.

[0065] Step 2, calculate the lag zero power of each range gate, compare each integration time (beam), and eliminate the range gate with a signal-to-noise ratio Q less than zero;

[0066] In this embodiment, the real and imaginary parts of the ACF at each lag period are obtained by cross-producting the responses 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 calculation formula is as follows:

[0067]

[0068] 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,

[0069] 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 integration.

[0070]

[0071] Specifically, only the ten lowest hysteresis power values ​​are used as the upper limit of the integration and cannot reflect the entire noise level. The reduction of the upper limit of the integration will cause the lower value in the average estimate to be biased, 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.

[0072] Furthermore, when receiving signals, the measurement results are inaccurate due to the overlap or interference of signals in multiple range cells. 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. The ratio μ is calculated, and when the ratio μ is greater than or equal to 0.3Navg, the corresponding interference lag is marked as "bad" and removed from the rawacf file.

[0073] Step 3: The fitacf file data obtained by least squares fitting of the rawacf data using the rst software package is organized into a complete data set, and further compared with the echo information obtained by high-resolution spectral analysis;

[0074] Step 4: Put the remaining autocorrelation function matrix into the built 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.

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

[0076] By first analyzing a large number of simulated ACFs, the speed and spectrum width for different numbers of delays are systematically calculated; 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 ​​constructed into the Hermite covariance matrix R by the formula p , the specific formula is as follows:

[0077]

[0078] 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,

[0079] Furthermore, R p Perform eigenvalue decomposition R p =EΛE H , by observing 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 a polynomial f(z) based on the eigenvector of the noise subspace, and the formula of the polynomial f(z) is as follows:

[0080]

[0081] 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.

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

[0083] Solve the root 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:

[0084]

[0085] 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.

[0086] Furthermore, by using the preprocessed R p Put it into the esprit algorithm model. Unlike the music algorithm, esprit uses the construction of signal subspace. Its characteristics are small calculation amount 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: U S =U S2 Φ, 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:

[0087]

[0088] Among them ∥·∥ 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 to find Φ = VΛV -1 , where V is a diagonal matrix, V -1 is the inverse matrix of V, Λ is a diagonal matrix, and 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.

[0089] Step 5: Use the ARMA model to construct the power spectrum estimation of the radar echo, obtain the power spectrum of 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;

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

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

[0092]

[0093] Among them, a i are the polynomial coefficients of the AR model, i is the order of the AR model, and b j is the polynomial coefficient b of the MA model i , j is the order of the MA model;

[0094] Specifically, by determining the order of the AR model and the MA model, we can determine the order corresponding to the modulus value of acf 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 also need to find the polynomial coefficients b of the MA model. i , the formula is as follows:

[0095]

[0096] Where R(0) is the value of the autocorrelation function at zero lag, and R(i) is the value of the autocorrelation function at time lag i;

[0097] Furthermore, when the coefficients ai of the AR model and the coefficients b of the MA model are obtained i After that, the modern spectrum P(f) of the ARMA model is constructed, and the formula is as follows:

[0098]

[0099] 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:

[0100] Exemplarily, the spectrum width is defined by the power spectrum using the full width at half maximum (FWHM), 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 the 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 forms a regular sequence, 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.

[0101] Step 6: Use the fanplot 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 by the criteria of v < 30 m / s and w < 35 m / s.

[0102] 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 of RMSE, NCC, and PPMCC are as follows:

[0103]

[0104] 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, K(i,j) represents the pixel value of image Y at the grid point, 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 the image b in the i-th row and j-th column, μ 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, respectively, 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.

[0105] The present invention is based on the development of spatial spectrum estimation and array signal technology, adopts a high-resolution spectral analysis method, utilizes the massive data with rich information obtained by the SuperDARN radar, expresses the relationship between the data with a high-resolution spectral analysis with powerful fitting ability, extracts the SuperDARN radar to distinguish different target echoes, and improves the radar's resolution ability.

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

[0107] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is 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 can still modify the technical solutions described in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A method for separating multiple echoes of SuperDARN radar based on high-resolution spectrum analysis, characterized in that: The following steps are involved: Step 1: Obtain rawacf file data of the raw data of the SuperDARN radar, and preprocess its autocorrelation function, wherein the preprocessing includes constructing a cross product for the ACF data of the rawacf file data, and performing phase unwrapping for the phase with a lag period exceeding 2π; Step 2: Calculate the hysteresis zero power of each range gate in the ACF data, and remove the range gates with a signal-to-noise ratio less than zero by comparing each integration time; Step 3: The fitacf file data obtained by least squares fitting of the rawacf data using the rst software package is organized into a complete data set and compared with the echo information obtained by 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 estimation of the radar echo, obtain the power spectrum of 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 fanplot 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 by 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 multiple echoes of SuperDARN radar based on high-resolution spectrum analysis according to claim 1, characterized in that: The preprocessing of the autocorrelation function also 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 represented as the phase, and for the phase φ(τ) that may cause phase wrapping, phase unwrapping is performed on the phase whose lag period exceeds 2π.

3. The method for separating multiple echoes of SuperDARN radar based on high-resolution spectrum analysis according to claim 2, characterized in that: The process of removing the range gates with a signal-to-noise ratio less than zero further comprises the following steps: Step 201: Based on the constructed cross product R(τ)=|R(τ)|e jφ(τ) The specific 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 according to the delayed zero power R(0) based on 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: Further increase the noise level by adding a correction factor μ, that is, when receiving the signal, by calculating 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 comparing the power of the two Calculate the ratio μ, and when the ratio μ is greater than or equal to 0.3Navg, mark the corresponding interference lag as "bad" and remove it from the rawacf file data.

4. The method for separating multiple echoes of SuperDARN radar 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: The dimension p of the autocorrelation matrix is ​​selected. p = 8 is selected based on a 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 , a polynomial f(z) is constructed based on the eigenvector of the noise subspace. 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 multiple echoes of SuperDARN radar 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 multiple echoes of SuperDARN radar 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: U S =U S2 Φ, 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 the least squares method to minimize the objective function 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 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 Once the Doppler frequency is obtained, it is converted into Doppler velocity.

7. The method for separating multiple echoes of SuperDARN radar 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 are the polynomial coefficients of the AR model, i is the order of the AR model, and b j is the polynomial coefficient b of the MA model i , j 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 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: Where R(0) is the value of the autocorrelation function at zero lag, and R(i) is the value of the autocorrelation function at time lag i; S3: When the coefficients ai of the AR model and the coefficients b of the MA model are obtained i After that, 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 power spectrum using the full width at half maximum (FWHM), 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 multiple echoes of SuperDARN radar based on high-resolution spectrum analysis according to claim 1, characterized in that: The calculation formulas of the root mean 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, K(i,j) represents the pixel value of image Y at the grid point, 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 the image b in the i-th row and j-th column, μ 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, respectively, 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