Ionized layer non-stationary interference suppression method and device and computer readable storage medium

Abstract

The invention relates to an ionosphere non-stationary interference suppression method, which comprises the following steps of: 1, forming a multi-polarization time-frequency beam; and step 2, fractional domain filtering. According to the method disclosed by the invention, information such as a time domain, a frequency domain and a space domain is comprehensively utilized, and the polarization domain information is combined, so that the information utilization rate is improved. Firstly, for the condition that the difference between the signal and the interference is small in the frequency domain, the spatial domain and the Doppler domain, interference suppression can be achieved through the polarization information difference between the signal and the interference, and compared with an existing algorithm, the method has wider applicability.

Description

Technical Field

[0001] The present invention belongs to the field of ionospheric detection, interference suppression and array signal processing, and in particular relates to an ionospheric non-stationary interference suppression method, device and computer-readable storage medium in this field. Background Art

[0002] Ionospheric sounding aims to study and monitor the structure and changes of the ionosphere above Earth and is a key area of ​​scientific research and engineering applications. The ionosphere is a layer of Earth's atmosphere formed by ionized particles generated by the interaction between atmospheric molecules and solar radiation. The ionosphere is of great significance for many applications, including weather, communications, and navigation. By transmitting sounding signals and receiving echo signals, ionospheric sounding can reveal various parameters of the ionosphere, such as the ionospheric height, ionospheric density, and ionospheric layer structure. These parameters are crucial for studying electromagnetic wave propagation, space weather changes, communication signal transmission, and satellite navigation. As an important component of the Earth-space environment, the spatiotemporal changes of the ionosphere have a significant impact on the propagation of electromagnetic waves in space. Accurately monitoring, forecasting, and warning of the morphology, structure, and disturbances of the ionosphere are of great significance to shortwave communications, aerospace, navigation, meteorological monitoring, disaster warning, and national defense.

[0003] Modern ionospheric detection often uses techniques such as pulse compression and coherent integration. However, incoming signals are susceptible to interference from the environment, and these interferences are significantly non-stationary. To ensure detection performance, effective identification and suppression of these interferences are crucial. However, many practical challenges remain, such as the difficulty in automatically identifying non-stationary interference in incoming signals and poor suppression effectiveness. To address these challenges in suppressing ionospheric echo interference, electromagnetic vector sensors can be introduced to improve the suppression of ionospheric non-stationary interference.

[0004] Electromagnetic vector sensor arrays (EMVSEs) are a new type of array capable of acquiring spatial and polarizational information about electromagnetic signals. Leveraging their array position differences and polarization sensitivity, EMVSEs can perform functions such as polarization parameter estimation, Direction of Arrival (DOA) estimation, and polarization-spatial joint beamforming. These functions can suppress interfering signals with unique spatial or polarizational characteristics within incoming signals. With the rapid advancement of array signal processing technology, researchers have achieved valuable results in parameter estimation based on EMVSEs. Beamforming and parameter estimation algorithms based on various EMVSEs, such as orthogonal electric dipoles, orthogonal tripoles, orthogonal three-magnetic rings, and full EMVSEs, have been proposed. However, traditional beamforming methods for EMVSEs have limited effectiveness in suppressing ionospheric non-stationary interference and suffer from poor accuracy in estimating the spatial angle and polarization parameters of incoming signals. Therefore, it is of practical significance to propose an EMVSE array-based method for suppressing ionospheric non-stationary interference that offers improved performance and can automatically identify interference.

[0005] In the prior art, the methods for suppressing ionospheric echo interference are mainly divided into the following categories:

[0006] (1) Ionospheric echo interference suppression method based on matrix eigendecomposition method. This method generates a covariance matrix based on the received echo signal, then performs eigenvalue decomposition on it and decomposes it into noise subspace and signal subspace. The original echo signal is projected onto the noise subspace to separate the interference signal part in the ionospheric echo, thereby achieving the suppression of the interference signal.

[0007] The paper “Ionospheric Clutter Suppression in High-Frequency Ground Wave Radar Based on Time-Frequency Analysis” (Xiong Xinnong, 2008, Systems Engineering and Electronic Technology) proposes using the joint eigenvalue decomposition method in the array element-pulse domain and the pulse-range domain to suppress ionospheric interference. The paper “Multimode Clutter Suppression Method Based on Eigendecomposition” (Weina, 2016, Journal of Radio Science) proposes using a template filter to perform sliding filtering on the echo signal to determine the number of echo signal modes. This allows the target signal and noise subspace to be determined by eigendecomposition of the covariance matrix, thereby completing clutter suppression.

[0008] (2) Ionospheric clutter echo interference suppression method based on adaptive sidelobe cancellation. This method sends the interference signals entering from the sidelobe of the main antenna and the interference signals entering from the auxiliary antenna to the adaptive processor at the same time, and then calculates the optimal weights according to the corresponding algorithm. The obtained optimal weights make the weighted output of each auxiliary channel just enough to cancel the interference received by the main channel, thus achieving the effect of suppressing ionospheric echo interference. For example, the paper “Design of Adaptive Sidelobe Cancellation Algorithm under Mainlobe Interference” (Song Hu, 2016, Systems Engineering and Electronic Technology) proposed a sidelobe cancellation system design idea based on subarray data differentiation, that is, using the properties of linear arrays to change the correlation between the main and auxiliary antenna signals, thereby eliminating the influence of the interference signal incident from the main lobe; the paper “BSS-based ABF for heterogeneous ionospheric clutter mitigation in HFSWR” (Guo Yueyu, 2019, IET Radar, Sonar & Navigation) proposed sidelobe strip ionospheric echo interference suppression based on blind source separation. By pre-processing blind source separation, the information of independent sources is obtained, which can achieve better interference suppression effect.

[0009] (3) Ionospheric interference suppression based on space-time adaptation. The design concept of the space-time adaptive filter is similar to that of the minimum variance distortionless response (MVDR) filter, which ensures the constant output of the target signal while minimizing the total output energy. The paper "Ionospheric clutter suppression method based on STAP" (Zhang Xin, 2013, Systems Engineering and Electronic Technology) proposed a local joint processing optimization method to estimate the covariance matrix of the echo interference and obtain the optimal space-time weight vector, effectively improving the ionospheric clutter suppression problem and target detection problem when the target signal direction and the ionospheric clutter are close. The paper "Ionospheric clutter suppression algorithm for high-frequency ground wave radar based on spatial multi-beam" (Yao Di, 2017, Journal of Electronics and Information Technology) proposed a high-frequency ground wave radar suppression method based on space multi-beam, which uses the strong correlation between beams and the characteristic that the number of range units (or Doppler units) covered by the clutter is much larger than the number of units covered by the target to achieve the estimation and suppression of ionospheric echo interference information.

[0010] In summary, existing ionospheric non-stationary interference suppression technologies, both domestically and internationally, primarily exploit differences between the signal and the interference in the frequency, spatial, and Doppler domains. Polarization domain information is rarely utilized, resulting in a waste of polarization domain information. Furthermore, if the differences between the signal and the interference in the frequency, spatial, and Doppler domains are small, the interference suppression effect is weakened or even ineffective. For detection signals in a linear frequency modulation (LFM) system, if the interference is a linear interference signal in the same frequency band and the spatial and Doppler differences between the signal and the interference are small, existing suppression methods are unable to effectively suppress it. The desired signal and interference signals in ionospheric echoes are mostly non-stationary, and their covariance matrices vary over time, making it difficult to determine appropriate weight vectors. Therefore, traditional beamforming methods have poor suppression effectiveness and are unable to identify interference within the echo signal, resulting in suppression failure. Because they cannot determine whether the echo signal contains interference or noise, existing non-adaptive methods often lack flexibility and are less effective at suppressing time-varying interference types. For adaptive methods that require iterative optimization, such as the adaptive sidelobe cancellation method, when the algorithm is iterating, within the convergence range, the larger the single iteration step size, the faster the convergence. However, if the single iteration step size is too large, the recursive process will oscillate. The disadvantage of traditional fixed-step algorithms is that they have a slow convergence speed and cannot overcome the inherent contradiction between convergence speed and steady-state error. Under the premise of convergence, if the single iteration step size is larger, although the convergence speed can be accelerated, the steady-state error will increase accordingly. Conversely, the steady-state error will decrease, but the convergence speed will slow down. Summary of the Invention

[0011] The present invention addresses the technical problem of difficulty in suppressing non-stationary interference in ionospheric echo signals and provides a method, device and computer-readable storage medium for suppressing non-stationary interference in the ionospheric echo signal. The method utilizes the time domain, frequency domain, spatial domain and polarization domain information in the echo signal to perform time-frequency analysis and beamforming to achieve adaptive polarization-spatial domain filtering, thereby effectively suppressing interference and noise with different polarization-spatial and time-frequency characteristics from the desired signal.

[0012] The present invention adopts the following technical solutions:

[0013] A method for suppressing ionospheric non-stationary interference, the improvement of which is that it comprises the following steps:

[0014] Step 1: Multi-polarization time-frequency beamforming:

[0015] Assume that the array receives a total of M incoming signals, the first signal is the target signal, and the remaining M-1 signals are interference signals. After beamforming, the array output is:

[0016]

[0017] In the above formula, w is the beamforming weight vector, w H =[w1,w2,…],(·) H stands for conjugate transpose, s H+V,m (t)=[s H,m (t),s V,m (t)] T , a H,m is called the horizontal polarization manifold vector, a V,m is called the vertical polarization manifold vector, s H (t) is the horizontal polarization component, s V (t) is the vertical polarization component, n(t) is the additive Gaussian white noise;

[0018] (1) The incoming target signal is a fully polarized wave:

[0019] When the incoming target signal is a fully polarized wave, we have:

[0020]

[0021] In the above formula, Represents the spatial steering vector of the mth incoming signal; represents the polarization vector of the mth incoming signal; s m (t) is the mth incoming signal; the number of rows of x(t) is 6L, corresponding to the output of all channels of L electromagnetic vector sensors;

[0022] Perform short-time Fourier transform on each received signal in the above equation, and we have:

[0023]

[0024] In the above formula, (·) T stands for transpose;

[0025] X(t,f)=[X1(t,f),X2(t,f),…,X 6L (t,f)] T (7)

[0026]

[0027] N(t,f)=[N1(t,f),N2(t,f),…,N 6L (t,f)] T (10)

[0028]

[0029] In the above formula, x l(t) represents the additive noise of the lth output signal, n l (t) represents the additive noise of the l-th output, X l (t,f) represents the time-frequency distribution of the l-th output, S m (t,f) represents the time-frequency distribution of the mth signal, N l (t,f) represents the time-frequency distribution of the additive noise in the l-th output;

[0030] Then, ignoring the noise, calculate

[0031]

[0032] In the above formula, D XX (t,f) is called the space-time-frequency distribution matrix, D SS (t,f) is called the signal time-frequency distribution matrix, D XX (t,f) is calculated as follows:

[0033]

[0034] For the l-th output, according to formula (8), we will get a two-dimensional matrix Create a cell array D XX (t,f), each unit of which is a two-dimensional matrix, calculated according to the following formula:

[0035]

[0036] In the above formula, "Θ" represents the Hadamard product. For 2×2 matrices a and b, we have:

[0037]

[0038] Assume that the total length of the signal to be short-time Fourier transform is N x , the number of points of the discrete window function is N win , the number of overlapping points between windows is N over ,but The dimension of Where T=(N x -N win ) / (N win -N over )+1, Represents rounding down to an integer;

[0039] When the incoming signal is a non-stationary signal, it is necessary to select the spatial time-frequency distribution matrix D XX The spatial time-frequency covariance matrix is ​​obtained by averaging the frequency points with high energy at all moments in (t, f);

[0040] The space-time-frequency covariance matrix is ​​estimated by the following operation:

[0041] First calculate the frequency judgment matrix D judge :

[0042]

[0043] D judge It is still a two-dimensional matrix related to time and frequency, with a dimension of In the frequency judgment matrix D judge For time point t k , select f k A frequency point with high energy, this f k The set of frequency points is recorded as Ω k :

[0044] Ω k ={(t k ,f n )|n=1,2,…,f k} (17)

[0045] The total number of selected frequency points is:

[0046] Let the set of these frequency points be Ω, and we have:

[0047] Ω={Ω k |k=1,2,…,T} (18)

[0048] Then in D XX (t,f) In the two-dimensional matrix of each unit, the frequency point is selected according to the position of the frequency point in the set Ω, so each unit will select F frequency points. For the same frequency point (t k ,f n ), the frequency point selected by each unit corresponds to a value, and all units form a 6L×6L dimension matrix D XX (t k ,f n ), the final space-time-frequency covariance matrix It is obtained by the following formula:

[0049]

[0050] The polarization parameters of the target signal are estimated by the following formula:

[0051]

[0052] γ est =arctan{|p est (2) / p est (1)|} (23)

[0053] η est =arg{p est (2) / p est (1)} (24)

[0054] In the above formula, V is the The noise subspace obtained by eigenvalue decomposition; is the spatial steering vector of the incoming target signal direction; represents the smallest generalized eigenvalue of the matrix bundle; "arg" represents the principal value of the argument; p est is the polarization vector of the target obtained by estimation; γ est is the polarization auxiliary angle of the target obtained by estimation; η est is the polarization phase difference of the target obtained by estimation;

[0055] Then the target's space-polarization domain steering vector is obtained by the following formula:

[0056]

[0057] Assuming that the weight vector used in beamforming is w1, which is a 6L×1 dimensional complex vector, the power of the beamforming output signal is When the output power in the (θ1, φ1) direction is constant, w1 is designed according to the following criteria:

[0058]

[0059] According to the corresponding relationship between the real part and the imaginary part, the constraints Rewritten as:

[0060]

[0061] "Re" and "Im" represent the real part and the imaginary part respectively, then:

[0062]

[0063] Construct the following Lagrangian function:

[0064]

[0065] In the above formula, l1 and l2 are two real-valued Lagrange multipliers. Let l=l1+jl2, then l * =l1-jl2,l+l * =2l1, simplifying the above formula to:

[0066]

[0067] According to formula (31):

[0068]

[0069] The partial derivative operation of the Lagrangian function J(w1,l) with respect to the Lagrangian multiplier l is defined as follows:

[0070]

[0071] Obviously, if equation (35) or equation (36) holds, then equation (33) and equation (34) must be zero at the same time;

[0072] Similarly, the Lagrangian function J(w1,l) is The definition of partial derivative is:

[0073]

[0074] The following derivation and The result is, where ν is a 6L×1 dimensional column vector, and It is still a 6L×1 dimensional column vector, and its lth elements are:

[0075]

[0076]

[0077] According to the definition of formula (37), we have:

[0078]

[0079] Let equation (37) be zero, and we get:

[0080]

[0081] According to formula (35), we can get:

[0082]

[0083] Solved from this The solved l H are real numbers, l and l H Same, so Substitute it back In the equation, the weight vector w1 is obtained as:

[0084]

[0085] This is the weight vector of multi-polarization time-frequency minimum power distortion-free response beamforming when the incoming target signal is a fully polarized wave; by minimizing Suppress interference;

[0086] (2) The incoming target signal is a partially polarized wave:

[0087] When the incoming target signal is a partially polarized wave, the polarization parameters change over time, and its polarization information cannot be used. Therefore, beamforming must be performed on the horizontally polarized component and the vertically polarized component of the array receive signal. The array receive signal at this time is:

[0088]

[0089] Perform short-time Fourier transform on each received signal in the above equation, and we have:

[0090]

[0091] In the above formula, S H,m (t,f) represents the time-frequency distribution of the horizontal polarization component of the mth signal, S V,m (t,f) represents the time-frequency distribution of the vertical polarization component of the mth signal;

[0092] Ignoring the noise, we have:

[0093]

[0094] In the above formula, D SS,H (t,f) is called the time-frequency distribution matrix of the horizontal polarization component of the signal, D SS,V (t,f) is called the time-frequency distribution matrix of the vertical polarization component of the signal;

[0095] Assuming that the weight vector used for beamforming is w2, which is a 6L×1 dimensional complex vector, we still need to obtain the space-time-frequency covariance matrix first.

[0096]

[0097] When (t k ,f n )∈S H / V,1 (t,f), there is S H / V,2 (t,f)≈0, When (t k ,f n )∈S H / V,2 (t,f), there is S H / V,1 (t k ,f n )≈0, Now define:

[0098]

[0099] Then we have:

[0100]

[0101] Design w2 according to the following guidelines:

[0102]

[0103] According to the corresponding relationship between the real part and the imaginary part, the constraints Rewritten as:

[0104]

[0105] Then we have:

[0106]

[0107] So construct the Lagrangian function:

[0108]

[0109] In the above formula, l 1,1 、l 1,2 、l 2,1 and l 2,2 are all real-valued Lagrange multipliers, let l1=l 1,1 +jl 1,2 , l2=l 2,1 +jl 2,2 ,have The above formula is simplified to:

[0110]

[0111] According to formula (62), we have:

[0112]

[0113] Define the Lagrangian function J(w2,l1,l2) and find the partial derivatives of the Lagrangian multipliers l1 and l2 as follows:

[0114]

[0115] If equations (68) to (71) hold, then equations (64) to (67) must be zero at the same time;

[0116] definition:

[0117]

[0118] In the above formula, l is the Lagrange multiplier vector, and formula (63) can be rewritten as:

[0119]

[0120] According to the conclusions of formula (40) and formula (41), the Lagrangian function J(w2,l) is Find the partial derivative and set it to zero:

[0121]

[0122] According to formula (72), we can get:

[0123]

[0124] Solved from this Substitute it back In the equation, the weight vector w2 is:

[0125]

[0126] This is the weight vector for multi-polarization time-frequency minimum power distortion-free response beamforming when the incoming target signal is a partially polarized wave;

[0127] By minimizing Suppress interference;

[0128] Step 2, fractional domain filtering:

[0129] For the signal s(t), its fractional Fourier transform is:

[0130]

[0131] In the above formula, “F p " represents the fractional Fourier transform operator of order p, where p is any real number; K p (t,u) represents the transformation kernel, as shown below:

[0132]

[0133] In the above formula, κ=pπ / 2, which is the rotation angle of the time-frequency axis on the time-frequency distribution diagram of s(t).

[0134] S p The inverse transform of (u) is:

[0135]

[0136] When p=0,±1,±2, the fractional Fourier transform degenerates into special transforms, namely: F0[s(t)]=s(t), F1[s(t)]=FT[s(t)], F -1 [s(t)]=FT -1 [s(t)],F ±2 [s(t)]=s(-t), where "FT" stands for Fourier transform, -1” stands for inverse Fourier transform;

[0137] The fractional Fourier transform has the following properties:

[0138]

[0139] Therefore, F4[s(t)]=F2[F2[s(t)]]=s(t), and the period of p is 4;

[0140] For a linear frequency modulation signal s LFM (t), there is an optimal order p0, so that the spectrum after the p0-order Fourier transform becomes an impulse function. Therefore, in the fractional domain, the output y(t) = w after the multi-polarization time-frequency minimum power distortion-free response beamforming H x(t) is detected. If there is a spectrum peak only at the p0 order, it is considered that the interference is completely filtered out;

[0141] The optimal weight p0 is given by:

[0142]

[0143] In the above formula, K LFM For s LFM (t) frequency modulation slope, T LFM For s LFM (t) duration, f s is the sampling frequency;

[0144] Assume that the score domain is detected at p i There is still interference at the order, and fractional domain filtering is required. The process is as follows:

[0145] First, perform p on y(t) i The order Fourier transform is:

[0146]

[0147] Then build an autoregressive linear model:

[0148]

[0149] In the above formula, d i,n is the coefficient of the autoregressive model, q is the order of the autoregressive model, and n(u) is the noise;

[0150] The above d i,n Calculated using the Yule-Walker equation:

[0151]

[0152] In the above formula,

[0153] Assume Y pi (u) The total number of points is U, and q is determined by calculating the minimum value of the FPE function, where the FPE function is:

[0154]

[0155] Finally, Y pi The impulse of (u) is treated as zero, and then predicted using formula (84) to obtain Finally, according to formula (80), the inverse p i Fourier transform, we can filter out p i The result of interference at the order;

[0156] The above operations are performed on different interference orders to obtain the final filtering result.

[0157] An ionospheric non-stationary interference suppression device is improved in that it includes: a processor; a memory storing executable instructions of the processor; and the processor is configured to perform the steps of the above method by executing the executable instructions.

[0158] A computer-readable storage medium is used to store a program, wherein the improvement is that the steps of the above method are implemented when the program is executed.

[0159] The beneficial effects of the present invention are:

[0160] The method disclosed in the present invention comprehensively utilizes information such as the time domain, frequency domain, and spatial domain, and combines polarization domain information to improve information utilization. First, for situations where the difference between the signal and the interference in the frequency domain, spatial domain, and Doppler domain is small, the difference in polarization information between the signal and the interference can be used to achieve interference suppression, which has a wider applicability than existing algorithms. Secondly, compared with traditional beamforming methods, the method of the present invention can achieve effective suppression when facing non-stationary interference in ionospheric echoes. Thirdly, the method of the present invention can identify interference signals in ionospheric echo signals and adaptively adjust the size of weight vectors to achieve adaptive polarization-spatial domain filtering; finally, the method of the present invention does not require iterative methods such as gradient descent, and has lower time complexity and higher stability than existing solutions. BRIEF DESCRIPTION OF THE DRAWINGS

[0161] Figure 1 It is a schematic diagram of the overall process of the method of the present invention;

[0162] Figure 2 It is the overall processing flow chart of the multi-polarization time-frequency beamformer algorithm;

[0163] Figure 3 This is a schematic diagram of the frequency judgment matrix when the target signal is a fully polarized signal;

[0164] Figure 4 This is a schematic diagram of the frequency judgment matrix when the target signal is a partially polarized signal. DETAILED DESCRIPTION

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

[0166] Example 1. This example discloses a method for suppressing ionospheric non-stationary interference. The method includes time-frequency domain feature analysis of interference signals in received waves. This method addresses the difficulty in automatically identifying interference signals from ionospheric detection waves. The method analyzes the time-frequency domain features of interference signals received by an electromagnetic vector array based on simulation data and measured data. Multi-polarization space-time-frequency filtering and interference suppression addresses the problem of poor suppression of interference signals from ionospheric detection waves. This method proposes a suppression method for fully polarized and partially polarized target waves.

[0167] like Figure 1 As shown, the specific steps include:

[0168] Step 1: Multi-polarization time-frequency beamforming:

[0169] like Figure 2 As shown in Figure 1, assuming that the array receives a total of M incoming signals, the first signal is the target signal, and the remaining M-1 signals are interference signals. After beamforming, the array output is:

[0170]

[0171] In the above formula, w is the beamforming weight vector, w H =[w1,w2,…],(·) H stands for conjugate transpose, s H+V,m (t)=[s H,m (t),s V,m (t)] T , a H,m is called the horizontal polarization manifold vector, a V,m is called the vertical polarization manifold vector, s H (t) is the horizontal polarization component, s V (t) is the vertical polarization component, n(t) is the additive Gaussian white noise;

[0172] The goal of beamforming is to minimize interference and noise while retaining the desired signal. The received amplitude pattern and polarization matching pattern should have a high main lobe pointing in the direction of the desired signal and low side lobes. At the same time, a deep null or even a zero point can be formed in the direction of the incoming interference wave. For fully polarized interference, the following are true:

[0173]

[0174] Whether the array receives a fully polarized or partially polarized signal can be determined by the coherence of its horizontal and vertical polarization components. When the target signal received by the array is a fully polarized wave, its polarization information can be utilized. However, when the target signal is a partially polarized wave, its polarization information varies with time and cannot be utilized. Different distortionless response filtering methods must be designed for these two situations.

[0175] (1) The incoming target signal is a fully polarized wave:

[0176] When the incoming target signal is a fully polarized wave, the target's polarization parameters must be estimated and the interference can be suppressed in the multi-polarization time-frequency domain using the estimated polarization parameters.

[0177]

[0178] In the above formula, Represents the spatial steering vector of the mth incoming signal; represents the polarization vector of the mth incoming signal; s m (t) is the mth incoming signal; the number of rows of x(t) is 6L, corresponding to the output of all channels of L electromagnetic vector sensors;

[0179] For the first wave, is the spatial steering vector of the incoming target signal direction, which is known; a H,1 is called the horizontal polarization manifold vector of the first incoming signal, a V,1 It is called the vertical polarization manifold vector of the first incoming signal.

[0180] Perform short-time Fourier transform on each received signal in the above equation, and we have:

[0181]

[0182] In the above formula, (·) T stands for transpose;

[0183] X(t,f)=[X1(t,f),X2(t,f),…,X 6L (t,f)] T (7)

[0184]

[0185] N(t,f)=[N1(t,f),N2(t,f),…,N 6L (t,f)] T (10)

[0186]

[0187] In the above formula, x l (t) represents the additive noise of the lth output signal, n l (t) represents the additive noise of the l-th output, X l (t,f) represents the time-frequency distribution of the l-th output, S m (t,f) represents the time-frequency distribution of the mth signal, N l (t,f) represents the time-frequency distribution of the additive noise in the l-th output, both of which are two-dimensional functions of time and frequency;

[0188] Then, ignoring the noise, calculate

[0189]

[0190] In the above formula, D XX (t,f) is called the space-time-frequency distribution matrix, D SS (t,f) is called the signal time-frequency distribution matrix, D XX The calculation method of (t,f) is as follows. Similarly, D can be calculated SS (t,f).

[0191]

[0192] In actual processing, each output of the array can only obtain a limited number of snapshots. For the l-th output, according to formula (8), a two-dimensional matrix will be obtained Therefore, when actually processing, we can create a cell array D XX (t,f), each unit of which is a two-dimensional matrix, calculated according to the following formula:

[0193]

[0194] In the above formula, "Θ" represents the Hadamard product. For 2×2 matrices a and b, we have:

[0195]

[0196] Assume that the total length of the signal to be short-time Fourier transform is N x , the number of points of the discrete window function is N win , the number of overlapping points between windows is Nover ,but The dimension of Where T=(N x -N win ) / (N win -N over )+1, Represents rounding down to an integer;

[0197] When the incoming signal is a non-stationary signal, it is necessary to select the spatial time-frequency distribution matrix D XX The spatial time-frequency covariance matrix is ​​obtained by averaging the frequency points with high energy at all moments in (t, f);

[0198] In actual processing, according to the estimated spatial time-frequency distribution matrix D XX Each unit of (t,f) is a two-dimensional matrix, and the space-time-frequency covariance matrix needs to be estimated by the following operations:

[0199] First calculate the frequency judgment matrix D judge :

[0200]

[0201] D judge It is still a two-dimensional matrix related to time and frequency, with a dimension of In the frequency judgment matrix D judge For time point t k , select f k A frequency point with high energy, this f k The set of frequency points is recorded as Ω k :

[0202] Ω k ={(t k ,f n )|n=1,2,…,f k} (17)

[0203] The total number of selected frequency points is:

[0204] Let the set of these frequency points be Ω, and we have:

[0205] Ω={Ω k |k=1,2,…,T} (18)

[0206] Then in D XX (t,f) In the two-dimensional matrix of each unit, the frequency point is selected according to the position of the frequency point in the set Ω, so each unit will select F frequency points. For the same frequency point (t k ,f n), the frequency point selected by each unit corresponds to a value, and all units form a 6L×6L dimension matrix D XX (t k ,f n ), the final space-time-frequency covariance matrix It is obtained by the following formula:

[0207]

[0208] Take a fully polarized interference as an example to illustrate the following: The process and form.

[0209] According to formula (12):

[0210]

[0211] for Figure 3 D shown judge , it is necessary to select frequency points along two time-frequency lines, which belong to S1(t,f) and S2(t,f) respectively.

[0212] When (t k ,f n )∈S1(t,f), then S2(t k ,f n )≈0, When (t k ,f n )∈S2(t,f), then S1(t k ,f n )≈0, At this time there are:

[0213]

[0214] According to formula (19), the space-time-frequency covariance matrix is ​​obtained Each element is D XX The (t,f) corresponds to the average of the high-energy frequency points in the time-frequency distribution matrix of the unit, and contains the spatial and polarization domain information of the signal. According to the subspace theory of multi-polarization arrays, the polarization parameters of the target signal are estimated by the following formula:

[0215]

[0216] γ est =arctan{|p est (2) / p est (1)|} (23)

[0217] η est =arg{p est (2) / pest (1)} (24)

[0218] In the above formula, V is the The noise subspace obtained by eigenvalue decomposition; is the spatial steering vector of the incoming target signal direction, which is known; represents the smallest generalized eigenvalue of the matrix bundle; "arg" represents the principal value of the argument; p est is the polarization vector of the target obtained by estimation; γ est is the polarization auxiliary angle of the target obtained by estimation; η est is the polarization phase difference of the target obtained by estimation;

[0219] Then the target's space-polarization domain steering vector is obtained by the following formula:

[0220]

[0221] Assuming that the weight vector used in beamforming is w1, which is a 6L×1 dimensional complex vector, the power of the beamforming output signal is When the output power in the (θ1, φ1) direction is constant, w1 is designed according to the following criteria:

[0222]

[0223] According to formula (21), minimize Equivalent to minimizing That is, the interference is automatically suppressed. This problem is solved by the Lagrange multiplier method.

[0224] According to the corresponding relationship between the real part and the imaginary part, the constraints Rewritten as:

[0225]

[0226] "Re" and "Im" represent the real part and the imaginary part respectively, then:

[0227]

[0228] Construct the following Lagrangian function:

[0229]

[0230] In the above formula, l1 and l2 are two real-valued Lagrange multipliers. Let l=l1+jl2, then l * =l1-jl2,l+l * =2l1, simplifying the above formula to:

[0231]

[0232] According to formula (31):

[0233]

[0234] The partial derivative operation of the Lagrangian function J(w1,l) with respect to the Lagrangian multiplier l is defined as follows:

[0235]

[0236] Obviously, if equation (35) or equation (36) holds, then equation (33) and equation (34) must be zero at the same time;

[0237] Similarly, the Lagrangian function J(w1,l) is The definition of partial derivative is:

[0238]

[0239] The following derivation and The result is, where ν is a 6L×1 dimensional column vector, and It is still a 6L×1 dimensional column vector, and its lth elements are:

[0240]

[0241] According to the definition of formula (37), we have:

[0242]

[0243] Let equation (37) be zero, and we get:

[0244]

[0245] According to formula (35), we can get:

[0246]

[0247] Solved from this The solved l H are real numbers, l and l H Same, so Substitute it back In the equation, the weight vector w1 is obtained as:

[0248]

[0249] This is the weight vector of multi-polarization time-frequency minimum power distortion-free response beamforming when the incoming target signal is a fully polarized wave; by minimizing Suppress interference;

[0250] (2) The incoming target signal is a partially polarized wave:

[0251] When the incoming target signal is a partially polarized wave, the polarization parameters change over time, and its polarization information cannot be used. Therefore, beamforming must be performed on the horizontally polarized component and the vertically polarized component of the array receive signal. The array receive signal at this time is:

[0252]

[0253] Perform short-time Fourier transform on each received signal in the above equation, and we have:

[0254]

[0255]

[0256] In the above formula, the definitions of X(t,f) and N(t,f) are the same as those when the incoming signal is a fully polarized wave. H,m (t,f) represents the time-frequency distribution of the horizontal polarization component of the mth signal, S V,m (t,f) represents the time-frequency distribution of the vertical polarization component of the mth signal, both of which are two-dimensional functions of time and frequency;

[0257] Ignoring the noise, we have:

[0258]

[0259] In the above formula, D SS,H (t,f) is called the time-frequency distribution matrix of the horizontal polarization component of the signal, D SS,V (t,f) is called the time-frequency distribution matrix of the vertical polarization component of the signal, D XX The definition of (t,f) is the same as that when the incoming signal is a fully polarized wave;

[0260] Assume that the weight vector used for beamforming is w2, which is a 6L×1 dimensional complex vector. First, we still need to use the method introduced in the case of fully polarized wave to obtain the spatial time-frequency covariance matrix. Using a partial polarization interference description we get The process and In this case, Equation (20) becomes:

[0261]

[0262] for Figure 4 D shown judge , it is necessary to select frequency points along two time-frequency lines, which belong to S H / V,1 (t,f) and S H / V,2(t,f).

[0263] When (t k ,f n )∈S H / V,1 (t,f), there is S H / V,2 (t,f)≈0, When (t k ,f n )∈S H / V,2 (t,f), there is S H / V,1 (t k ,f n )≈0, Now define:

[0264]

[0265] Then we have:

[0266]

[0267] Design w2 according to the following guidelines:

[0268]

[0269] The problem is still solved by the Lagrange multiplier method. Note that there are two constraints at this time. According to the correspondence between the real part and the imaginary part, the constraints are Rewritten as:

[0270]

[0271] Then we have:

[0272]

[0273] So construct the Lagrangian function:

[0274]

[0275] In the above formula, l 1,1 、l 1,2 、l 2,1 and l 2,2 are all real-valued Lagrange multipliers, let l1=l 1,1 +jl 1,2 , l2=l 2,1 +jl 2,2 ,have The above formula is simplified to:

[0276]

[0277] According to formula (62), we have:

[0278]

[0279] Define the Lagrangian function J(w2,l1,l2) and find the partial derivatives of the Lagrangian multipliers l1 and l2 as follows:

[0280]

[0281] If equations (68) to (71) hold, then equations (64) to (67) must be zero at the same time;

[0282] definition:

[0283]

[0284] In the above formula, l is the Lagrange multiplier vector, and formula (63) can be rewritten as:

[0285]

[0286] According to the conclusions of formula (40) and formula (41), the Lagrangian function J(w2,l) is Find the partial derivative and set it to zero:

[0287]

[0288] According to formula (72), we can get:

[0289]

[0290] Solved from this Substitute it back In the equation, the weight vector w2 is:

[0291]

[0292] This is the weight vector for multi-polarization time-frequency minimum power distortion-free response beamforming when the incoming target signal is a partially polarized wave;

[0293] By minimizing Suppress interference;

[0294] Step 2, fractional domain filtering:

[0295] Considering that interference suppression solely in the multi-polarization time-frequency domain generally cannot completely eliminate interference, it is possible to further consider secondary interference suppression on the PTF-MPDR output based on the characteristics of the incoming signal in other domains. Since the incoming target signal is a linear frequency modulated signal, it can form an impulse pulse in the fractional domain of a specific order. This property allows fractional-domain detection of the PTF-MPDR (multi-polarization time-frequency minimum power distortionless response beamforming) output. The detection results determine whether further processing is required. If spectral peaks are present only at the order corresponding to the target signal, the interference is considered nearly completely suppressed and no further processing is required. If significant spectral peaks at other orders exist in addition to the impulse at the target order, further fractional-domain filtering of these peaks can be considered. The following describes the fractional-domain detection and filtering process.

[0296] First, let’s introduce the concept of fractional Fourier transform. For a signal s(t), its fractional Fourier transform is:

[0297]

[0298] In the above formula, “F p " represents the fractional Fourier transform operator of order p, where p is any real number; K p (t,u) represents the transformation kernel, as shown below:

[0299]

[0300] In the above formula, κ = pπ / 2, which is the rotation angle of the time-frequency axis on the time-frequency distribution diagram of s(t), which can also be called the rotation angle of the fractional Fourier transform;

[0301] S p The inverse transform of (u) is:

[0302]

[0303] It can be seen from the above formula that the signal s(t) is composed of a set of weight coefficients S p (u) orthogonal basis function K -p (t,u) representation, these basis functions are complex exponential functions in the form of linear frequency modulation.

[0304] When p=0,±1,±2, the fractional Fourier transform degenerates into special transforms, namely: F0[s(t)]=s(t), F1[s(t)]=FT[s(t)], F -1 [s(t)]=FT -1 [s(t)],F ±2 [s(t)]=s(-t), where "FT" stands for Fourier transform, -1 ” stands for inverse Fourier transform;

[0305] The fractional Fourier transform has the following properties:

[0306]

[0307] Therefore, F4[s(t)]=F2[F2[s(t)]]=s(t), and the period of p is 4;

[0308] For a linear frequency modulation signal s LFM (t), there is an optimal order p0, so that the spectrum after the p0-order Fourier transform becomes an impulse function. Therefore, in the fractional domain, the output y(t) = w after the multi-polarization time-frequency minimum power distortion-free response beamforming H If the spectrum peak is only at the p0 order, the interference is considered to be completely filtered out. It should be noted that for s LFM (t) for example, LFM (t) = s LFM (-t), so we only need to scan the part p∈[0,2].

[0309] The optimal weight p0 is given by:

[0310]

[0311] In the above formula, K LFM For s LFM (t) frequency modulation slope, T LFM For s LFM (t) duration, f s is the sampling frequency;

[0312] Assume that the score domain is detected at p i There is still interference at the order, and fractional domain filtering is required. The process is as follows:

[0313] First, perform p on y(t) i The order Fourier transform is:

[0314]

[0315] Then build an autoregressive linear model:

[0316]

[0317] In the above formula, d i,n is the coefficient of the autoregressive model, q is the order of the autoregressive model, and n(u) is the noise;

[0318] The above d i,n Calculated using the Yule-Walker equation:

[0319]

[0320] In the above formula,

[0321] Assumptions The total number of points is U, and q is determined by calculating the minimum value of the FPE (Final Prediction Error) function, where the FPE function is:

[0322]

[0323] Finally, The impulse of is treated as zero, and then predicted using formula (84) to obtain Finally, according to formula (80), the inverse p i Fourier transform, we can filter out p i The result of interference at the order;

[0324] The above operations are performed on different interference orders to obtain the final filtering result.

[0325] This embodiment also discloses an ionospheric non-stationary interference suppression device, comprising: a processor; a memory storing executable instructions of the processor; and the processor being configured to perform the steps of the above method by executing the executable instructions.

[0326] This embodiment also discloses a computer-readable storage medium for storing a program, which implements the steps of the above method when executed.

Claims

1. A method for suppressing ionospheric non-stationary interference, characterized in that: The steps include: Step 1: Multi-polarization time-frequency beamforming: Assume that the array receives a total of M incoming signals, the first signal is the target signal, and the remaining M-1 signals are interference signals. After beamforming, the output of the array is: In the above formula, w is the beamforming weight vector, w H =[w1,w2,…],(·) H stands for conjugate transpose, s H+V,m (t) = [s H,m (t),s V,m (t)] T , a H,m is called the horizontal polarization manifold vector, a V,m is called the vertical polarization manifold vector, s H (t) is the horizontal polarization component, s V (t) is the vertical polarization component, n(t) is the additive Gaussian white noise; (1) The incoming target signal is a fully polarized wave: When the incoming target signal is a fully polarized wave, we have: In the above formula, Represents the spatial steering vector of the mth incoming signal; represents the polarization vector of the mth incoming signal; s m (t) is the mth incoming signal; the number of rows of x(t) is 6L, corresponding to the output of all channels of L electromagnetic vector sensors; Each received signal in the above formula is subjected to short-time Fourier transform, and we have: In the above formula, (·) T stands for transpose; X(t,f)=[X1(t,f),X2(t,f),…,X 6L (t,f)] T (7) N(t,f)=[N1(t,f),N2(t,f),…,N 6L (t,f)] T (10) In the above formula, x l (t) represents the additive noise of the lth output signal, n l (t) represents the additive noise of the l-th output, X l (t,f) represents the time-frequency distribution of the lth output, S m (t,f) represents the time-frequency distribution of the mth signal, N l (t,f) represents the time-frequency distribution of the additive noise in the l-th output; Then, ignoring the noise, we calculate In the above formula, D XX (t,f) is called the space-time-frequency distribution matrix, D SS (t,f) is called the signal time-frequency distribution matrix, D XX (t,f) is calculated as follows: For the l-th output, according to formula (8), we will get a two-dimensional matrix Create a cell array D XX (t,f), each unit of which is a two-dimensional matrix, calculated according to the following formula: In the above formula, "Θ" represents the Hadamard product. For 2×2 dimensional matrices a and b, we have: Assume that the total length of the signal to be short-time Fourier transformed is N x , the number of points of the discrete window function is N win , the number of overlapping points between windows is N over ,but The dimension of Where T = (N x -N win ) / (N win -N over )+1, Represents rounding down; When the incoming signal is a non-stationary signal, it is necessary to select the spatial time-frequency distribution matrix D XX The spatial time-frequency covariance matrix is ​​obtained by averaging the frequency points with high energy at all moments in (t,f); The space-time-frequency covariance matrix is ​​estimated by the following operation: First, calculate the frequency judgment matrix D judge : D judge It is still a two-dimensional matrix related to time and frequency, with a dimension of In the frequency point judgment matrix D judge For time point t k , select f k A frequency point with high energy, this f k The set of frequency points is recorded as Ω k : Ω k ={(t k ,f n )|n=1,2,…,f k } (17) The total number of selected frequencies is: Let the set of these frequency points be Ω, and we have: Oh={Oh k |k=1,2,…,T} (18) Then in D XX (t,f) In the two-dimensional matrix of each unit, the frequency point is selected according to the position of the frequency point in the set Ω, so each unit will select F frequency points. For the same frequency point (t k ,f n ), each unit selects a frequency point corresponding to a value, and all units form a 6L×6L-dimensional matrix D XX (t k ,f n ), the final space-time-frequency covariance matrix It is obtained by the following formula: The polarization parameters of the target signal are estimated by the following formula: γ est =arctan{|p est (2) / p est (1)|} (23) or est =arg{p est (2) / p est (1)} (24) In the above formula, V is the The noise subspace obtained by eigenvalue decomposition; is the spatial steering vector of the incoming target signal direction; represents the smallest generalized eigenvalue of the matrix bundle; "arg" represents the argument principal value; p est is the polarization vector of the target obtained by estimation; est is the polarization auxiliary angle of the target obtained by estimation; η est is the polarization phase difference of the target obtained by estimation; Then the target's null-polarization domain steering vector is obtained by the following formula: Assuming that the weight vector used in beamforming is w1, which is a 6L×1 dimensional complex vector, the power of the beamforming output signal is When the output power in the (θ1, φ1) direction is constant, w1 is designed according to the following criteria: According to the corresponding relationship between the real part and the imaginary part, the constraints Rewritten as: "Re" and "Im" represent the real part and the imaginary part respectively, then: Construct the following Lagrangian function: In the above formula, l1 and l2 are two real-valued Lagrange multipliers. Let l = l1 + jl2, then l * =l1-jl2,l+l * =2l1, simplify the above formula to: According to formula (31), The partial derivative operation of the Lagrangian function J(w1,l) with respect to the Lagrangian multiplier l is defined as follows: Obviously, if equation (35) or equation (36) holds, then equation (33) and equation (34) must be zero at the same time; Similarly, the Lagrangian function J(w1,l) is The definition of partial derivative is: The following derivation and The result is, where ν is a 6L×1-dimensional column vector, and It is still a 6L×1 dimensional column vector, and its lth elements are: According to the definition of formula (37), we have: Let (37) be zero, and we get: According to formula (35), we can get: Solved from this The solved l H are real numbers, l and l H Same, so Substitute it back In the equation, the weight vector w1 is obtained as: This is the weight vector of multi-polarization time-frequency minimum power distortion-free response beamforming when the incoming target signal is a fully polarized wave; by minimizing Suppress interference; (2) The incoming target signal is a partially polarized wave: When the incoming target signal is a partially polarized wave, the polarization parameters will change with time, and its polarization information cannot be used. It is necessary to beamform the horizontal polarization component and the vertical polarization component of the array receiving signal respectively. At this time, the array receiving signal is: Each received signal in the above formula is subjected to short-time Fourier transform, and we have: In the above formula, S H,m (t,f) represents the time-frequency distribution of the horizontal polarization component of the mth signal, S V,m (t,f) represents the time-frequency distribution of the vertical polarization component of the mth signal; Without considering the noise, we have: In the above formula, D SS,H (t,f) is called the time-frequency distribution matrix of the horizontal polarization component of the signal, D SS,V (t,f) is called the time-frequency distribution matrix of the vertical polarization component of the signal; Assuming that the weight vector used for beamforming is w2, which is a 6L×1 dimensional complex vector, we still need to obtain the space-time-frequency covariance matrix When (t k ,f n )∈S H / V,1 (t,f), we have When (t k ,f n )∈S H / V,2 (t,f), we have Now define: Then we have: Design w2 according to the following guidelines: According to the corresponding relationship between the real part and the imaginary part, the constraints Rewritten as: Then we have: So construct the Lagrangian function: In the above formula, l 1,1 , l 1,2 , l 2,1 and l 2,2 are all real-valued Lagrange multipliers, let l1 = l 1,1 +jl 1,2 , l2=l 2,1 +jl 2,2 ,have The above formula is simplified to: According to formula (62), we have: Define the Lagrangian function J(w2,l1,l2) and find the partial derivatives of the Lagrangian multipliers l1 and l2 as follows: If equations (68) to (71) hold, then equations (64) to (67) must be zero at the same time; definition: In the above formula, l is the Lagrange multiplier vector, then formula (63) can be rewritten as: According to the conclusions of equations (40) and (41), the Lagrangian function J(w2,l) is Find the partial derivative and set it to zero: According to formula (72), we can get: Solved from this Substitute it back In the equation, the weight vector w2 is obtained as: This is the weight vector of multi-polarization time-frequency minimum power distortion-free response beamforming when the incoming target signal is a partially polarized wave; By minimizing Suppress interference; Step 2, fractional domain filtering: For the signal s(t), its fractional Fourier transform is: In the above formula, "F p " represents a fractional Fourier transform operator of order p, where p is an arbitrary real number; K p (t,u) represents the transformation kernel, as shown below: In the above formula, κ = pπ / 2, which is the rotation angle of the time-frequency axis on the time-frequency distribution diagram of s(t). S p The inverse transform of (u) is: When p = 0, ±1, ±2, the fractional Fourier transform degenerates into special transforms, namely: F0[s(t)] = s(t), F1[s(t)] = FT[s(t)], F -1 [s(t)]=FT -1 [s(t)],F ±2 [s(t)] = s(-t), where "FT" stands for Fourier transform, -1 ” stands for inverse Fourier transform; The fractional Fourier transform has the following properties: Therefore, F4[s(t)]=F2[F2[s(t)]]=s(t), and the period of p is 4; For a linear frequency modulated signal s LFM (t), there is an optimal order p0, so that its spectrum after p0-order Fourier transform becomes an impulse function. Therefore, in the fractional domain, the output y(t) = w after multi-polarization time-frequency minimum power distortion-free response beamforming H x(t) is detected. If there is a spectrum peak only at the p0 order, it is considered that the interference is completely filtered out; The optimal weight p0 is given by: In the above formula, K LFM For LFM (t) frequency modulation slope, T LFM For LFM (t) duration, f s is the sampling frequency; Assume that the score domain is detected at p i There is still interference at the order, and fractional domain filtering is required. The process is as follows: First, p y(t) i The order Fourier transform is: Then build an autoregressive linear model: In the above formula, d i,n is the coefficient of the autoregressive model, q is the order of the autoregressive model, and n(u) is the noise; The above i,n Calculated using the Yule-Walker equation: In the above formula, Assumptions The total number of points is U, and q is determined by calculating the minimum value of the FPE function, where the FPE function is: Finally, The impulse of is treated as zero, and then predicted using equation (84) to obtain Finally, according to formula (80), the inverse p i Fourier transform, we can filter out p i The result of interference at the order; The above operation is performed on different interference orders to obtain the final filtering result.

2. An ionospheric non-stationary interference suppression device, characterized in that: include: processor; A memory in which executable instructions of the processor are stored; the processor is configured to perform the steps of the method of claim 1 by executing the executable instructions.

3. A computer-readable storage medium for storing a program, characterized in that: When the program is executed, the steps of the method according to claim 1 are implemented.

Images