CN113238193A - Multi-component combined reconstruction SAR echo broadband interference suppression method - Google Patents

Multi-component combined reconstruction SAR echo broadband interference suppression method Download PDF

Info

Publication number
CN113238193A
CN113238193A CN202110443842.4A CN202110443842A CN113238193A CN 113238193 A CN113238193 A CN 113238193A CN 202110443842 A CN202110443842 A CN 202110443842A CN 113238193 A CN113238193 A CN 113238193A
Authority
CN
China
Prior art keywords
time
frequency
sar echo
interference
component
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110443842.4A
Other languages
Chinese (zh)
Other versions
CN113238193B (en
Inventor
周峰
梅诗国
韩文畅
樊伟伟
田甜
石晓然
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN202110443842.4A priority Critical patent/CN113238193B/en
Publication of CN113238193A publication Critical patent/CN113238193A/en
Application granted granted Critical
Publication of CN113238193B publication Critical patent/CN113238193B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/36—Means for anti-jamming, e.g. ECCM, i.e. electronic counter-counter measures
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

The invention belongs to the technical field of radar signal processing, and discloses a multi-component combined reconstruction SAR echo broadband interference suppression method; the method is used for constructing a two-dimensional time-frequency diagram based on short-time Fourier transform of self-adaptive window width and estimating instantaneous frequency; in addition, ridge path detection based on a Viterbi algorithm and the proposed re-detection method is designed, and instantaneous frequency estimation is carried out; finally, the multi-component combined reconstruction broadband interference suppression method is provided for broadband interference suppression of the synthetic aperture radar, can effectively suppress broadband interference, and lays a certain foundation for synthetic aperture radar image acquisition in an interference environment.

Description

Multi-component combined reconstruction SAR echo broadband interference suppression method
Technical Field
The invention relates to the technical field of Radar signal processing, in particular to a multi-component combined reconstruction SAR echo broadband interference suppression method, which is used for recovering Synthetic Aperture Radar (SAR) echoes and effectively suppressing broadband interference so as to obtain a high-quality SAR image.
Background
The synthetic aperture radar can work all day long and all day long, and is widely applied in the fields of national defense and military (battlefield reconnaissance, situation tracking and the like), national economy (such as transportation, weather forecast, resource detection and the like) and scientific research (such as aerospace, atmospheric physics, celestial body research and the like).
As electromagnetic environments become more complex, electromagnetic interference inevitably destroys the SAR echo. These interferences may be referred to as radio frequency interferences, which degrade the signal-to-interference ratio of the measurement data and seriously affect the imaging quality of the synthetic aperture radar system. In the presence of radio frequency interference, the estimation of the doppler parameters will be inaccurate, which may lead to image defocus. Furthermore, there are some artifacts in SAR images, making it difficult to detect and identify targets. If the SAR data for radio frequency interference suppression is discarded, a great waste of time and resources is caused. Therefore, it is necessary to develop an effective interference suppression method by making full use of precious synthetic aperture radar data.
These interferences can be classified into two categories according to their bandwidths: narrowband interference (NBI) and Wideband interference (WBI). NBI means that the bandwidth of the interference is only a small fraction of the SAR echo bandwidth. In contrast, WBI occupies a large portion of the SAR echo bandwidth.
Over the past few years, researchers have proposed a number of approaches to inhibit both NBI and WBI, some of which have been developed from methods of inhibiting NBI, and therefore, methods of NBI inhibition have been briefly introduced, followed by methods of WBI inhibition. Some researchers use mathematical models to characterize NBI. NBI was modeled by l.nguyen and m.soumekh as the sum of some complex sine waves. J.yi and x.wan et al derive a closed form of approximate maximum likelihood estimation of the estimated interference parameters. However, model errors can lead to a significant degradation in performance. To this end, some researchers have proposed methods based on filter or transform domain analysis. Reigber and l.ferro-Fami propose notch filters for suppressing interference. Wang and W.Y et al propose another method to obtain the radio frequency interference subspace by approximate spectral decomposition, effectively reconstruct the radio frequency interference signal, and then remove the interference from the originally received SAR pulse. Feng and h.zheng et al apply a new method of using subband-spectral cancellation. Zhou and r.wu propose suppression of NBI based on feature subspace filtering. M.spencer and c.chen propose slow time thresholding methods for detecting and suppressing ground radio frequency interference. Smith and r.d.hill propose Least Mean Square (LMS) filters and wiener (Weiner) filters. Similarly, c.t.c.le, s.hensley and e.chapin propose LMS adaptive filters for suppressing radio frequency interference in broadband synthetic aperture radar signals. In addition, y.huang and g.liao et al also propose three methods of jointly utilizing NBI sparsity and low rank property.
For NBI, the bandwidth of NBI is relatively small and signal loss after interference suppression is tolerable. However, this is not applicable to WBI, which is much more complex to inhibit than NBI. There are some non-parametric WBI suppression methods, such as temporal spectral notch filters and subspace projections, unlike the corresponding NBI suppression methods, m.tao and f.zhou indicate that synthetic aperture radar pulses are transformed to the two-dimensional time-frequency domain by short-time fourier transform (STFT), and that notch and feature subspace-based filtering can be used for the temporal spectrum. S.zhang and m.xing propose an interference suppression method based on time-frequency transformation and wavelet transformation. Yang and w.du et al propose a WBI suppression method using an iterative adaptive method and an orthogonal subspace projection method. The main difficulty with these non-parametric methods is that the characteristics of WBI are not fully exploited, which may lead to signal loss in the synthetic aperture radar echo. The tensor technique of suppressing both narrowband and wideband interference is proposed by y.huang and l.zhang, et al, which has a good effect on complex interference. However, this method employs a multi-view approximation of single channel data, which is based on the assumption that the interference has approximately stable frequencies or frequency bands within one synthetic aperture time. Therefore, this method is not suitable for variable-band interference.
Disclosure of Invention
Aiming at the problems in the prior art, the invention aims to provide a multi-component combined reconstruction SAR echo broadband interference suppression method, which is characterized in that a two-dimensional time-frequency graph is constructed based on STFT (self-adaptive window width), ridge path detection based on a Viterbi (Viterbi) algorithm and a redetection method is designed and used for precise estimation of instantaneous frequency; finally, the multi-component combined reconstruction broadband interference suppression method is provided for broadband interference suppression of the synthetic aperture radar, can effectively suppress broadband interference, and lays a certain foundation for image acquisition of the synthetic aperture radar.
The basic idea of the invention is as follows: first, synthetic aperture radar returns in the presence of broadband interference are converted to a two-dimensional time-frequency representation by STFT. Secondly, estimating the instantaneous frequency of the interference component through the obtained two-dimensional time-frequency representation; the instantaneous frequency estimation can be divided into ridge path detection and ridge path recombination. Afterwards, a multi-component joint reconstruction method is utilized to extract WBI.
In order to achieve the purpose, the invention is realized by adopting the following technical scheme.
A multi-component combined reconstruction SAR echo broadband interference suppression method comprises the following steps:
step 1, carrying out short-time Fourier transform of self-adaptive window width on an original SAR echo signal to obtain a corresponding two-dimensional time-frequency graph;
step 2, carrying out instantaneous frequency estimation on the interference component on the two-dimensional time-frequency diagram to obtain an instantaneous frequency estimation value of the interference component;
and 3, performing joint reconstruction on the instantaneous frequency estimation value of the interference component by using a multi-component joint reconstruction method to complete broadband interference suppression on the original SAR echo, and obtaining the SAR echo after the broadband interference suppression, namely the reconstructed SAR echo.
Compared with the prior art, the invention has the beneficial effects that:
(1) the time-frequency diagram is obtained by adopting the STFT with the adaptive window width, the time-frequency localization information of the signal can be better provided, and the high-resolution time-frequency diagram can be obtained by adopting the adaptive window width.
(2) The invention carries out instantaneous frequency estimation on a plurality of components in the WBI, and comprises ridge path detection and ridge path recombination by using a generalized Viterbi algorithm and a redetection method, so that the obtained instantaneous frequency estimation value of each WBI component well reveals the essential attribute of interference, and provides powerful support for subsequent WBI filtering.
(3) The invention designs a broadband interference suppression method of multi-component combined reconstruction, filters WBI with high precision and recovers SAR echo.
Drawings
The invention is described in further detail below with reference to the figures and specific embodiments.
FIG. 1 is a flow chart of a method for suppressing wideband interference based on multi-component joint reconstruction according to the present invention;
FIG. 2 is a time-frequency diagram and a frequency spectrum of three different types of WBI examples constructed according to an embodiment of the present invention, where (a) is a time-frequency diagram of a chirp WBI; (b) is the spectrum of the linear frequency modulation WBI; (c) is a time-frequency diagram of a generalized sinusoidal frequency modulation WBI; (d) is a frequency spectrum of a generalized sinusoidal frequency modulation WBI; (e) is a time-frequency diagram of an arbitrary curve frequency modulation WBI; (f) frequency-modulating the WBI spectrum for an arbitrary curve;
FIG. 3 is a schematic diagram of frequency spectrum and time-frequency diagrams of measured SAR echoes and SAR echoes with WBI constructed according to an embodiment of the present invention; wherein (a) is the spectrum of the SAR echo; (b) a time-frequency graph of the SAR echo is obtained; (c) is the spectrum of the SAR echo containing WBI; (d) a time-frequency diagram of SAR echoes containing WBI;
FIG. 4 is a schematic diagram of ridge path detection and ridge path reconstruction constructed in accordance with an embodiment of the present invention; wherein, (a) is an instantaneous frequency estimation graph obtained for a ridge path based on a Viterbi algorithm and a re-detection method; (b) a plot of instantaneous frequency crossover area for a disconnect; (c) reconstructing a result graph for a ridge path;
fig. 5 is a schematic diagram of four WBI components extracted from a SAR echo containing WBI constructed according to an embodiment of the present invention;
FIG. 6 is a schematic diagram of a recovered SAR echo time-frequency graph constructed according to an embodiment of the present invention;
FIG. 7 is a diagram illustrating SAR echoes and recovery results thereof with WBI at a signal-to-interference ratio of-3 dB, constructed in accordance with an embodiment of the present invention; wherein, (a) is a time-frequency diagram of the SAR echo in the presence of WBI; (b) is a graph of the instantaneous frequency estimation result in this case; (c) a time-frequency diagram of the recovered SAR echo;
FIG. 8 is a diagram of measured radar echo imaging results performed in accordance with an embodiment of the present invention; wherein, (a) the image is a non-interference SAR echo imaging result image, and (b) the image is a SAR echo imaging result image with-12 dB signal-to-interference ratio WBI;
FIG. 9 is a graph comparing the results of four methods performed in accordance with one embodiment of the present invention; the method comprises the steps of (a) obtaining an instantaneous frequency spectrum notch suppression result graph, (b) obtaining an instantaneous characteristic subspace filtering suppression result graph, (c) obtaining an interference suppression method suppression result graph based on time-frequency transformation and wavelet transformation, and (d) obtaining a multi-component combined reconstruction method suppression result graph;
Detailed Description
Embodiments of the present invention will be described in detail below with reference to examples, but it will be understood by those skilled in the art that the following examples are only illustrative of the present invention and should not be construed as limiting the scope of the present invention.
Referring to fig. 1, the method for suppressing the SAR echo broadband interference through multi-component joint reconstruction provided by the present invention includes the following steps:
step 1, performing short-time Fourier transform (STFT) of self-adaptive window width on an original SAR echo signal to obtain a corresponding two-dimensional time-frequency graph;
fourier transforms are widely used for signal analysis and processing, and can reveal the characteristics of a signal in the frequency domain. However, the original fourier transform is performed over the entire time of the signal, and time information is lost. As an extension of the fourier transform, STFT provides time-frequency localization information of a signal, which is a common method in time-frequency analysis. STFT is widely used for generating time-frequency maps, and unlike quadratic time-frequency distributions, which have no cross terms and are more suitable for multi-component signal analysis, the basic idea of STFT is to perform a fourier transform on a portion of the signal, usually using a smooth window for weighting.
(1.1) for a discrete-time signal, the STFT may be expressed as:
Figure BDA0003036025450000061
where j denotes an imaginary unit, N denotes a time sampling number, k denotes a frequency sampling number, N is a sample number, and s (N) and h (-) denote a discrete time signal and a window function, respectively.
In the present invention, STFT is used to generate a time-frequency diagram, as shown in fig. 2, which shows a time-frequency diagram and a frequency spectrum diagram of WBI of different types (chirp, generalized sinusoidal chirp, arbitrary curve chirp), the frequency spectrum is generated by normalizing Fast Fourier Transform (FFT) results with the number of samples, and as can be seen from the figure, the broadband interference has a wide bandwidth, the frequency spectrum cannot describe frequency variation in detail, and the two-dimensional time-frequency curve can more accurately represent the broadband interference.
The window width is crucial for obtaining a high resolution time-frequency representation, and in order to obtain a higher frequency resolution, a long window should be chosen, but this reduces the time resolution. Therefore, there must be a compromise between frequency resolution and time resolution.
(1.2) the Gaussian window is widely used in time-frequency analysis because it has the smallest time-bandwidth product, and therefore the STFT is performed using the Gaussian window as a window function. The Gaussian window is defined as
Figure BDA0003036025450000071
Where n is the time sample number and σ is the standard deviation. The gaussian window is defined on (— infinity, + ∞), but it must be truncated to achieve numerical computation. The present invention contemplates a Gaussian window over [ -3 σ,3 σ ]. For a gaussian window, the standard deviation can be considered as a measure of the window width.
(1.3) in order to obtain a time-frequency graph with high aggregation degree, solving the following optimization problem:
Figure BDA0003036025450000072
wherein CM represents a measure of concentration，STFTs，σIs a Gaussian window STFT of the signal s whose standard deviation is σ, and σoptIs the best standard deviation.
In order to compress the search range and improve the efficiency, the standard deviation is optimized under the logarithmic scale; the specific process is as follows:
(1.3a) setting the search interval [ sigma ] in advance1，σ2]Then the search interval on a logarithmic scale is converted to [ log ]uσ1，loguσ2]Wherein u > 1;
(1.3b) dividing the search space [ log ]uσ1，loguσ2]Dividing the cell into N-1 cells with the same length, and then calculating the aggregation degree measure of two end points of each cell; in the present invention, let σ 11 and σ2N/6, where N is the number of samples.
The concentration measure can be used for evaluating the quality of the time-frequency graph, and the concentration measure under STFT is expressed as:
Figure BDA0003036025450000073
wherein the content of the first and second substances,
Figure BDA0003036025450000081
is a time-frequency diagram WsNormalized representation of (n, k), i.e.
Figure BDA0003036025450000082
Wherein the content of the first and second substances,
Figure BDA0003036025450000083
using STFT to generate time-frequency diagram, therefore let Ws(n，k)＝STFTs(n，k)。
In this embodiment, α is 0.5 and β is 1.1.
(1.3c) constructing a new search interval by adopting two side points of the point corresponding to the maximum value of the concentration degree measure, wherein the new search interval is smaller than the search interval in the step (1.3 a);
repeating steps (1.3b) - (1.3c) until the length of the new search interval is less than a given threshold loguEta, using the finally obtained new search interval as the optimal standard deviation sigmaoptAnd taking the optimal standard deviation as the standard deviation of a Gaussian window, and performing STFT on the original SAR echo signal to obtain a corresponding high-aggregation time-frequency graph.
The above iterative process is explained below: the number of iterations of the above-described optimal standard deviation optimization process is independent of u, as demonstrated in detail below. The length of the initial interval can be written as:
Figure BDA0003036025450000084
for the next iteration, the interval length becomes
Figure BDA0003036025450000085
Wherein, a ', b', σ1' and σ2' is an updated variable. An interval shrinkage rate of
Figure BDA0003036025450000086
The iteration will stop in the following cases
Figure BDA0003036025450000087
The above formula is equivalent to (consider u > 1)
Figure BDA0003036025450000091
The interval shrinkage rate and the termination condition of the algorithm and the number of iterations are independent of u. In the invention, u is 10, and a time-frequency diagram with high aggregation degree is finally obtained.
Step 2, carrying out instantaneous frequency estimation of interference components on the two-dimensional time-frequency diagram to obtain instantaneous frequency estimation;
comprising the following substeps:
(2.1) spine path detection based on a Viterbi algorithm and a re-detection method;
the instantaneous frequency can be estimated by detecting the ridges from the time-frequency diagram, which is non-parametric and adaptive. A simple ridge path detection method is represented as:
Figure BDA0003036025450000092
wherein the content of the first and second substances,
Figure BDA0003036025450000093
representing the ridge path, is an estimate of the true instantaneous frequency, Ws(n, k) is a time-frequency diagram of the signal, IfIs an indexed set of frequency points.
However, the obtained ridge path may be distorted in a noisy environment. Furthermore, for multi-component signals, the method may extract points due to other components. Using the (generalized) Viterbi algorithm a smooth ridge path is obtained which passes through points in the time-frequency diagram with large amplitudes.
(2.1a) time-frequency diagram W of the original SAR echo signal using Viterbi algorithms(n, k) performing ridge path detection to obtain an initial ridge path and a corresponding instantaneous frequency estimate.
For a time interval n1，n2]The set P may be defined as n1To n2The set of all possible paths in between.
The instantaneous frequency estimate may be expressed as
Figure BDA0003036025450000101
Wherein L (p (n); n)1，n2) Represents a time interval [ n ]1，n2]Middle path p: (n), L (p (n); n is1，n2) Is the penalty function f (W (n, p (n))) and g (p (n), p (n +1)) in the time interval [ n1，n2]The sum of the inner paths. f (-) and g (-) are penalty functions, W, respectivelys(n, p (n)) represents a time-frequency diagram;
the points with large amplitude of the time-frequency diagram reveal the instantaneous frequency to a great extent, and the penalty function f (x) can be constructed according to the following method: the time-frequency graph values at the time index n are ordered in a non-increasing order:
W(n，k1)≥W(n，k2)≥…≥W(n，kj′)≥…≥W(n，kM)
kj′∈If，j′＝1，2，…，M
wherein M is IfThe number of elements in (b), j' is an index in the sorted sequence, and the function f (x) can be expressed as:
f(Ws(n，kj′))＝j′-1
using function f (x), points with larger magnitudes are more likely to be selected as candidate points for instantaneous frequency estimation. Another penalty function g (p (n), p (n +1)) is related to the smoothness of the path. If the function g (x, y) is constant, the instantaneous frequency estimation is reduced to a simple ridge detection method, but this may be severely affected by noise. To obtain a smoother path, g (x, y) may be expressed as
Figure BDA0003036025450000102
Where δ is a threshold value representing the maximum instantaneous frequency variation allowed between successive points, and ρ is a penalty factor. If the change is less than δ, it is considered reasonable with no penalty. Otherwise, the penalty will be proportional to the fraction that exceeds the threshold. The threshold δ may be set to a value corresponding to several points to obtain a desired result. For δ → ∞, g (x, y) ═ 0, at which time the instantaneous frequency estimation is simplified to a simple ridge path detection method. ρ controls the degree of penalty. In the present invention, δ is 2 and ρ is 10.
For having M frequenciesPoint and time-frequency diagram of N moments, the possible path number between two end points is MNThis makes it very difficult to directly search for an optimal path. The present invention uses a recursive method based on the generalized Viterbi algorithm, which greatly reduces the computational burden. The method comprises the following concrete steps:
first, for a time interval [ n ]1，ni]The partial best path is represented as:
Figure BDA0003036025450000111
wherein, Pij′Is comprised of a time n1And point (n)1，kj′) A set of all possible paths in between; the corresponding instantaneous frequency estimates are:
Figure BDA0003036025450000112
secondly, for the next instant ni+1The partial best path is represented as:
pi+1(n，kj′)＝[pi(n；kl)，(ni+1，kj′)]
wherein l can be obtained by the following formula
Figure BDA0003036025450000113
For time n > n1For M points at time n, M points need to be searched2A candidate point, which causes a large computational burden; to reduce computational complexity, the space [0, M-1 ] is searched]Modified as [ j '-delta, j' + delta]. Thus, a smooth ridge path can be obtained through the recursive implementation of the generalized Viterbi algorithm, and the ridge path is used as the estimation of the instantaneous frequency, and the calculation time is greatly reduced.
As shown in fig. 3(a) and (b), are the frequency spectrum and time-frequency plots of the measured SAR echo, while for multi-component signals, as shown in fig. (c) and (d), SAR echoes containing WBIs may have time-frequency cross-components. In this case, the paths are detected one by one, and the frequency points around the path just detected are set to zero to reduce the influence on the remaining component detection.
However, the instantaneous frequency estimate around the crossover point is distorted by the gap. The present invention proposes a re-detection method to reduce distortion.
(2.1b) estimate for initial instantaneous frequency
Figure BDA0003036025450000121
From the time-frequency diagram W of the corresponding raw SAR echo signali(n, k) constructing a new time-frequency representation W'i(n，k)：
Figure BDA0003036025450000122
(2.1c) representing W 'for each new time-frequency using Viterbi Algorithm'i(n, k) carrying out ridge path detection in the step (2.1a) to obtain a ridge path after the redetection, namely the instantaneous frequency estimation after the redetection.
After the above re-detection method is adopted, the distortion of instantaneous frequency estimation around the ridge path intersection point is well reduced, and finally the ridge path after the time-frequency image re-detection is obtained.
(2.2) cutting off and recombining the redetected ridge paths to obtain new ridge paths, wherein each new ridge path corresponds to a final instantaneous frequency estimation of a WBI component;
for a multi-component signal, as shown in fig. 4(a), each detected ridge path may contain a different component because the direction of change of the ridge path is not considered in the above method. A ridge-path-slope based reorganization approach is therefore used to solve this problem. The resulting ridge paths are recombined as shown in fig. 4(b) and (c).
The intersection points of the re-examined ridge paths are cut off, and then re-connected by using a re-combination method based on the ridge path slope to form new ridge paths, as shown in fig. 5, so as to obtain the final instantaneous frequency estimate of each WBI component.
And 3, performing joint reconstruction on the final instantaneous frequency estimation of each WBI component by using a multi-component joint reconstruction method to complete broadband interference suppression on the original SAR echo, and obtaining the SAR echo after the broadband interference suppression, namely the reconstructed SAR echo.
The instantaneous frequency reveals the essential nature of the interference and a filter based on the instantaneous frequency can concentrate the interference energy and filter it out.
In general, a discretely sampled SAR echo containing Mc interference components may be represented as
Figure BDA0003036025450000131
Wherein, IFm(u) an instantaneous frequency estimate corresponding to the mth interference component; n is the number of samples, z (t) is the sum of useful echo and noise,
Figure BDA0003036025450000132
is the complex envelope of the mth interference component, which, assuming it is a band-limited signal, can be expanded to
Figure BDA0003036025450000133
Wherein the content of the first and second substances,
Figure BDA0003036025450000134
complex fourier coefficients that are the complex envelope of the mth component; f. of0＝fs/(QN) is the fundamental frequency, Q is a positive integer, fsIs the sampling frequency. In the present invention, Q is 4.
Substituting the above formula into the expression of s (t) to obtain
Figure BDA0003036025450000135
Wherein the content of the first and second substances,
Figure BDA0003036025450000136
for convenience, the above formula can also be written as
s＝Ax+z
Wherein s ═ s (t)0)，…，s(tN-1)]T，z＝[z(t0)，…，z(tN-1)]T，
Figure BDA0003036025450000137
Figure BDA0003036025450000138
(·)TDenotes transposition, A ═ A1，…，AM]Matrix AmThe p row and q column elements of (1) are:
(Am)pq＝exp(j2π(q-K-1)f0tp-1+jψm(tp-1))
p＝1，2，…，N；q＝1，2，…，2K+1
the coefficients are estimated taking into account the following optimization problem:
Figure BDA0003036025450000139
wherein | · | purple sweet2And | · | non-conducting phosphor1Respectively represent l2Norm and l1Norm, λ is the regularization parameter.
Figure BDA00030360254500001310
Is a data fidelity term, | | x | | non-woven phosphor1Is a sparse regularization term, λ is a constant, and these two terms are weighted. In fact, for the problem considered here, the solution of the above optimization problem is not sensitive to λ, and it is not difficult to select a suitable λ empirically, let λ be 1 in the present invention.
l1The norm minimization problem can be solved by a Fast Iterative Shrinkage Threshold Algorithm (FISTA) Algorithm. To findThe specific process of the solution is as follows: is provided with
Figure BDA0003036025450000141
The following approximate model was used:
Figure BDA0003036025450000142
wherein x is(n)Is the solution for the nth iteration; r (x)(n)) Is dependent on x(n)A relative constant of (d);<>representing an included angle; xi is the gradient
Figure BDA0003036025450000143
The Lipschitz (Lipschitz) constant of (A), the smallest Lipschitz constant beingHTwice the maximum eigenvalue of a.
Due to l1The norm is separable and the problem can be reduced to a series of one-dimensional minimization problems for each component of x, the solution being expressed as:
Figure BDA0003036025450000144
wherein, Tτ(.) is the shrink operator, i.e.
Tτ(xi)＝max(|xi|-τ，0)sign(xi)
Where sign () is a sign function and τ is a parameter of the shrink operator.
The following summarizes the FISTA-based envelope complex fourier coefficient estimation method.
First, s, A, xi, M are inputteditEpsilon is more than 0, xi and epsilon are constants respectively; and initializes a variable x(0)＝x(1)＝0，α(0)＝α(1)＝1。
Next, iteration is performed according to the following formula
Figure BDA0003036025450000145
x(n+1)＝pξ(y(n))
Wherein n ═ {1,2, …, Mit}，MitThe maximum number of iterations is indicated.
In that
Figure BDA0003036025450000151
The iteration is terminated when necessary, otherwise the order is
Figure BDA0003036025450000152
The iteration continues.
Outputting the result after the iteration is terminated
Figure BDA0003036025450000153
After envelope complex fourier coefficient estimation, each interference component may be reconstructed as
Figure BDA0003036025450000154
Wherein the content of the first and second substances,
Figure BDA0003036025450000155
is the sub-vector corresponding to the mth WBI component.
Figure BDA0003036025450000156
Can be directly from
Figure BDA0003036025450000157
Extracting.
The overall interference can be estimated as
Figure BDA0003036025450000158
The SAR echo after removing the interference can be obtained by the following method
Figure BDA0003036025450000159
The SAR echo after removing the interference is shown in fig. 6. This method is essentially a time-frequency filtering method, since the instantaneous frequency can be represented as a time-frequency curve on a time-frequency plane, and a sparse regularization term is used in the method. The filtering area is near the instantaneous frequency in the time-frequency plane and has the bandwidth of
B＝2Kf0
For relatively high signal-to-interference ratios, e.g., -3dB, the time-frequency representation of the SAR echo containing WBI is shown in fig. 7 (a). In this case, the time-frequency curve of the WBI component is mixed with the points where the SAR echo and noise are generated in certain areas, which makes the instantaneous frequency estimation more difficult. Fig. 7(b) shows the instantaneous frequency estimate obtained in this case, and fig. 7(c) shows the time-frequency representation of the recovered SAR echo. In this case, the recovery result of the SAR echo indicates the robustness of the multi-component joint reconstruction method proposed by the present invention to perform wideband interference.
Experimental verification is carried out on the multi-component combined reconstruction method in the embodiment, and different WBI inhibition methods are adopted to compare inhibition performance results; to quantitatively evaluate the performance of the proposed method, the Normalized Recovery Error (NRE) of the SAR echo is used, which can be expressed as:
Figure BDA0003036025450000161
where Z is a matrix containing the original SAR echo as its column,
Figure BDA0003036025450000162
is a matrix containing recovered SAR echoes, | · | | non-woven phosphorFThe Frobenius norm of the matrix is represented.
As shown in fig. 8(a), it is a true SAR echo imaging result graph without interference, and the four broadband interference suppression methods of the present invention are compared: transient spectrum notch, transient feature subspace filtering, interference suppression methods based on time-frequency transformation and wavelet transformation, and multi-component joint reconstruction methods provided by the invention.
WBI suppression contrast experiments of-12 dB signal-to-interference ratios were performed for each of the four methods in the above example, and the imaging results of the SAR echo corrupted by WBI are shown in fig. 8 (b).
Fig. 9(a) - (d) are WBI inhibition results for the four methods, respectively. The NRE of the transient frequency spectrum notch is-5.07 dB, the NRE of the transient characteristic subspace filtering is-6.85 dB, the NRE of the interference suppression method based on time-frequency transformation and wavelet transformation is-2.75 dB, and the NRE of the multi-component joint reconstruction method provided by the invention is-10.53 dB, so that the multi-component joint reconstruction method provided by the invention has the best suppression effect, and the effectiveness of the invention is verified.
Although the present invention has been described in detail in this specification with reference to specific embodiments and illustrative embodiments, it will be apparent to those skilled in the art that modifications and improvements can be made thereto based on the present invention. Accordingly, such modifications and improvements are intended to be within the scope of the invention as claimed.

Claims (10)

1. A multi-component combined reconstruction SAR echo broadband interference suppression method is characterized by comprising the following steps:
step 1, carrying out short-time Fourier transform of self-adaptive window width on an original SAR echo signal to obtain a corresponding two-dimensional time-frequency graph;
step 2, carrying out instantaneous frequency estimation on the interference component on the two-dimensional time-frequency diagram to obtain an instantaneous frequency estimation value of the interference component;
and 3, performing joint reconstruction on the instantaneous frequency estimation value of the interference component by using a multi-component joint reconstruction method to complete broadband interference suppression on the original SAR echo, and obtaining the SAR echo after the broadband interference suppression, namely the reconstructed SAR echo.
2. The method for suppressing the SAR echo broadband interference of the multi-component joint reconstruction according to claim 1, wherein the short-time Fourier transform of the adaptive window width is performed on the original SAR echo signal, and the specific steps are as follows:
(1.1) for a discrete-time signal, the short-time fourier transform is represented as:
Figure FDA0003036025440000011
wherein j represents an imaginary unit, N represents a time sampling number, k represents a frequency sampling number, N is a sample number, and s (N) and h (-) represent a discrete time signal and a window function, respectively;
(1.2) performing a short-time Fourier transform using a Gaussian window as a window function, the Gaussian window being defined as:
Figure FDA0003036025440000012
wherein σ is the standard deviation, and the range of the Gaussian window is [ -3 σ,3 σ ];
bringing the Gaussian window into the step (1.1) to obtain the short-time Fourier transform of the Gaussian window of the signal s;
(1.3) in order to obtain a time-frequency graph with high aggregation degree, optimizing the standard deviation of a Gaussian window, namely solving the following optimization problem:
Figure FDA0003036025440000021
wherein CM represents a measure of concentration, STFTs,σIs a short-time Fourier transform of a Gaussian window of the signal s whose standard deviation is σ, and σoptIs the standard deviation after optimization.
3. The method for suppressing the SAR echo broadband interference of the multi-component joint reconstruction according to claim 2, wherein the solving process of the optimization problem is as follows: the standard deviation optimization is carried out under a logarithmic scale, and the specific process is as follows:
(1.3a) setting the search interval [ sigma ] in advance1,σ2]Then the search interval on a logarithmic scale is converted to [ log ]uσ1,loguσ2]Wherein u is>1；
(1.3b) dividing the search space [ log ]uσ1,loguσ2]Dividing the cell into N-1 cells with the same length, and then calculating the aggregation degree measure of two end points of each cell; wherein N is the number of samples;
the calculation formula of the concentration measure is as follows;
Figure FDA0003036025440000022
wherein the content of the first and second substances,
Figure FDA0003036025440000023
is a time-frequency diagram WsNormalized representation of (n, k), i.e.
Figure FDA0003036025440000024
Wherein the content of the first and second substances,
Figure FDA0003036025440000025
Ws(n,k)＝STFTs(n, k), alpha and beta are respectively preset parameters and are constants;
(1.3c) constructing a new search interval by adopting two side points of the point corresponding to the maximum value of the concentration degree measure, wherein the new search interval is smaller than the search interval in the step (1.3 a);
repeating steps (1.3b) - (1.3c) until the length of the new search interval is less than a given threshold loguEta, using the finally obtained new search interval as the optimal standard deviation sigmaoptAnd taking the optimal standard deviation as the standard deviation of a Gaussian window, and carrying out short-time Fourier transform on the original SAR echo signal to obtain a corresponding high-aggregation two-dimensional time-frequency graph.
4. The method for suppressing the SAR echo broadband interference of the multi-component joint reconstruction according to claim 1, wherein the instantaneous frequency estimation of the interference component is performed on the two-dimensional time-frequency graph, specifically:
(2.1) time-frequency diagram W of original SAR echo signal by using Viterbi algorithms(n, k) performing ridge path detection to obtain an initial ridge path and a corresponding instantaneous frequency estimate;
(2.2) estimate for initial instantaneous frequency as
Figure FDA0003036025440000031
From the time-frequency diagram W of the corresponding raw SAR echo signali(n, k) constructing a new time-frequency representation Wi′(n,k)：
Figure FDA0003036025440000032
Where δ represents the maximum instantaneous frequency change threshold allowed between successive points;
(2.3) representing W for each new time-frequency using Viterbi algorithmi' (n, k) performing ridge path detection in the step (2.1a) to obtain a ridge path after the redetection, namely, the ridge path is an instantaneous frequency estimation after the redetection;
and (2.4) cutting off and recombining the redetected ridge paths to obtain new ridge paths, wherein each new ridge path corresponds to a final instantaneous frequency estimation of a WBI component.
5. The SAR echo broadband interference suppression method for multi-component joint reconstruction according to claim 4, characterized in that the time-frequency graph W of the original SAR echo signal by using Viterbi algorithms(n, k) performing ridge path detection, specifically comprising the following steps:
first, a search space is set to [ j '- δ, j' + δ](ii) a For a time interval n1,ni]The partial best path is represented as:
Figure FDA0003036025440000033
Figure FDA0003036025440000041
wherein L (p (n); n)1,n2) Represents a time interval [ n ]1,(ni,kj')]Penalty function of the medium path P (n), Pij'Is comprised of a time n1And point (n)1,kj') A set of all possible paths p (n) in between; f (-) and g (-) are penalty functions, W, respectivelys(n, p (n)) represents a time-frequency diagram;
the corresponding instantaneous frequency estimates are:
Figure FDA0003036025440000042
secondly, for the next instant ni+1The partial best path is represented as:
pi+1(n,kj')＝[pi(n；kl),(ni+1,kj')]
wherein l is obtained by the following formula
Figure FDA0003036025440000043
6. The method for suppressing the SAR echo broadband interference of multi-component joint reconstruction according to claim 5, wherein the construction process of the penalty function f (x) is as follows:
the time-frequency graph values at the time index n are ordered in a non-increasing order:
W(n,k1)≥W(n,k2)≥…≥W(n,kj')≥…≥W(n,kM)
kj'∈If,j'＝1,2,…,M
wherein M is IfThe number of elements in (b), j' is an index in the sorted sequence, and the function f (x) is expressed as:
f(Ws(n,kj'))＝j'-1；
the penalty function g (x, y) is expressed as:
Figure FDA0003036025440000044
where δ is a threshold representing the maximum instantaneous frequency variation allowed between successive points, and ρ is a penalty factor, both of which are constants.
7. The SAR echo broadband interference suppression method based on multi-component joint reconstruction according to claim 4, characterized in that the cutting and reconstruction are performed on the re-detected ridge path, specifically: firstly, cutting off the intersection points of the ridge paths after the re-inspection, and then, re-connecting the intersection points by using a re-assembly method based on the slope of the ridge paths to form new ridge paths.
8. The method for suppressing the SAR echo wideband interference through multi-component joint reconstruction according to claim 1, wherein the joint reconstruction is performed on the instantaneous frequency estimation value of the interference component by using a multi-component joint reconstruction method, specifically:
(3.1) will contain McThe discretely sampled SAR echo for each interference component is represented as:
Figure FDA0003036025440000051
wherein, IFm(u) an instantaneous frequency estimate corresponding to the mth interference component; n is the number of samples, z (t) is the sum of useful echo and noise,
Figure FDA0003036025440000052
is the complex envelope of the mth interference component, which is a band-limited signal, and is expanded as:
Figure FDA0003036025440000053
wherein the content of the first and second substances,
Figure FDA0003036025440000054
complex fourier coefficients that are the complex envelope of the mth component; f. of0＝fs/(QN) is the fundamental frequency, Q is a positive integer, fsIs the sampling frequency;
substituting the above formula into the expression of s (t), and then including McThe discretely sampled SAR echo for each interference component is represented as:
Figure FDA0003036025440000055
wherein the content of the first and second substances,
Figure FDA0003036025440000056
the above equation is written as:
s＝Ax+z
wherein s ═ s (t)0),…,s(tN-1)]T，z＝[z(t0),…,z(tN-1)]T，
Figure FDA0003036025440000061
Figure FDA0003036025440000062
(·)TDenotes transposition, A ═ A1,…,AM]Matrix AmThe p row and q column elements of (1) are:
(Am)pq＝exp(j2π(q-K-1)f0tp-1+jψm(tp-1))
p＝1,2,…,N；q＝1,2,…,2K+1
(3.2) estimating coefficients considering the following optimization problem:
Figure FDA0003036025440000063
wherein | · | purple sweet2And | · | non-conducting phosphor1Respectively represent l2Norm and l1Norm, λ is a regularization parameter, constant;
(3.3) solving the optimization problem by fast iterative shrinkage threshold algorithm, i.e./1Norm minimization problem, yielding an envelope complex Fourier coefficient estimate
Figure FDA0003036025440000064
(3.4) estimation by enveloping the complex Fourier coefficients
Figure FDA0003036025440000065
And reconstructing each interference component and obtaining the SAR echo after interference suppression.
9. The method for suppressing the SAR echo broadband interference of the multi-component joint reconstruction according to claim 8, wherein the optimization problem is solved by a fast iterative shrinkage threshold algorithm, and the specific steps are as follows:
first, s, A, xi, M are inputtedit,ε>0, xi and epsilon are constants respectively; and initializes a variable x(0)＝x(1)＝0,α(0)＝α(1)＝1。
Next, iteration is performed according to the following formula:
Figure FDA0003036025440000066
x(n+1)＝pξ(y(n))
wherein n ═ {1,2, …, Mit}，MitRepresenting the maximum number of iterations;
Figure FDA0003036025440000067
x(n)is the solution for the nth iteration;
Figure FDA0003036025440000068
the superscript H denotes the conjugate transpose, Tτ(. is) a shrink operator, i.e.
Tτ(xi)＝max(|xi|-τ,0)sign(xi)
Wherein sign () is a sign function, τ is a parameter of the shrink operator;
finally, whether the iteration termination condition is met or not is judged
Figure FDA0003036025440000071
If yes, outputting the result
Figure FDA0003036025440000072
Otherwise make
Figure FDA0003036025440000073
The iteration continues.
10. The method for suppressing the SAR echo wideband interference of multi-component joint reconstruction according to claim 8, characterized in that the estimation by enveloping complex Fourier coefficients
Figure FDA00030360254400000710
Reconstructing each interference component, and obtaining an SAR echo after interference suppression, specifically:
first, each interference component is reconstructed as:
Figure FDA0003036025440000074
wherein the content of the first and second substances,
Figure FDA0003036025440000075
is a sub-vector corresponding to the mth WBI component, directly from
Figure FDA0003036025440000076
Extracting;
the overall interference is estimated as:
Figure FDA0003036025440000077
SAR echo after interference suppression
Figure FDA0003036025440000078
Expressed as:
Figure FDA0003036025440000079
CN202110443842.4A 2021-04-23 2021-04-23 SAR echo broadband interference suppression method based on multi-component joint reconstruction Active CN113238193B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110443842.4A CN113238193B (en) 2021-04-23 2021-04-23 SAR echo broadband interference suppression method based on multi-component joint reconstruction

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110443842.4A CN113238193B (en) 2021-04-23 2021-04-23 SAR echo broadband interference suppression method based on multi-component joint reconstruction

Publications (2)

Publication Number Publication Date
CN113238193A true CN113238193A (en) 2021-08-10
CN113238193B CN113238193B (en) 2023-06-30

Family

ID=77129036

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110443842.4A Active CN113238193B (en) 2021-04-23 2021-04-23 SAR echo broadband interference suppression method based on multi-component joint reconstruction

Country Status (1)

Country Link
CN (1) CN113238193B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115015849A (en) * 2022-08-09 2022-09-06 西安电子科技大学 Interference suppression method and related device
CN116679277A (en) * 2023-07-26 2023-09-01 中国科学院空天信息创新研究院 Minimum L based on time-frequency domain 1 Double-channel SAR multi-jammer positioning method based on norm

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20110113926A (en) * 2010-04-12 2011-10-19 한국과학기술원 Synthetic aperture radar system for continuous wide swath high resolution imaging and method thereof
CN111239697A (en) * 2020-02-11 2020-06-05 西北工业大学 Multidimensional domain combined SAR broadband interference suppression method based on low-rank matrix decomposition
CN111273238A (en) * 2020-01-06 2020-06-12 中国航天科工集团八五一一研究所 SAR (synthetic aperture radar) wide-band and narrow-band interference simultaneous inhibition method based on low-rank recovery
CN112269168A (en) * 2020-10-14 2021-01-26 西安电子科技大学 SAR broadband interference suppression method based on Bayesian theory and low-rank decomposition

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20110113926A (en) * 2010-04-12 2011-10-19 한국과학기술원 Synthetic aperture radar system for continuous wide swath high resolution imaging and method thereof
CN111273238A (en) * 2020-01-06 2020-06-12 中国航天科工集团八五一一研究所 SAR (synthetic aperture radar) wide-band and narrow-band interference simultaneous inhibition method based on low-rank recovery
CN111239697A (en) * 2020-02-11 2020-06-05 西北工业大学 Multidimensional domain combined SAR broadband interference suppression method based on low-rank matrix decomposition
CN112269168A (en) * 2020-10-14 2021-01-26 西安电子科技大学 SAR broadband interference suppression method based on Bayesian theory and low-rank decomposition

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
YAN HUANG等: "Simultaneous Narrowband and Wideband Interference Suppression on Single-Channel SAR System via Low-Rank Recovery", 《IGARSS 2019》 *
孔舒亚等: "基于压缩感知的SAR宽带干扰抑制方法", 《电子测量技术》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115015849A (en) * 2022-08-09 2022-09-06 西安电子科技大学 Interference suppression method and related device
CN116679277A (en) * 2023-07-26 2023-09-01 中国科学院空天信息创新研究院 Minimum L based on time-frequency domain 1 Double-channel SAR multi-jammer positioning method based on norm
CN116679277B (en) * 2023-07-26 2023-10-17 中国科学院空天信息创新研究院 Minimum L based on time-frequency domain 1 Double-channel SAR multi-jammer positioning method based on norm

Also Published As

Publication number Publication date
CN113238193B (en) 2023-06-30

Similar Documents

Publication Publication Date Title
Hao et al. A joint framework for multivariate signal denoising using multivariate empirical mode decomposition
CN111273238B (en) SAR (synthetic aperture radar) wide-band and narrow-band interference simultaneous inhibition method based on low-rank recovery
CN110221256B (en) SAR interference suppression method based on deep residual error network
CN113238193B (en) SAR echo broadband interference suppression method based on multi-component joint reconstruction
CN112269168A (en) SAR broadband interference suppression method based on Bayesian theory and low-rank decomposition
CN111507047B (en) Inverse scattering imaging method based on SP-CUnet
Ren et al. RFI mitigation for UWB radar via hyperparameter-free sparse SPICE methods
CN111965615A (en) Radar target detection method based on estimation before detection
CN113534151B (en) Dual-band ISAR imaging method based on off-grid sparse Bayesian learning
CN113671498A (en) SAR radio frequency interference suppression method based on low-rank and dual sparse matrix decomposition
Han et al. Wideband interference suppression for SAR via instantaneous frequency estimation and regularized time-frequency filtering
CN109815849A (en) Chaotic signal Denoising Algorithm based on singular value decomposition
Gao et al. Image compressive sensing reconstruction based on z-score standardized group sparse representation
CN113887398A (en) GPR signal denoising method based on variational modal decomposition and singular spectrum analysis
CN116699526A (en) Vehicle millimeter wave radar interference suppression method based on sparse and low-rank model
Yang et al. Robust and efficient harmonics denoising in large dataset based on random SVD and soft thresholding
Chang et al. Parameter estimation for ultrasonic echo signals through improved matching pursuit and flower pollination algorithms
CN111046791A (en) Current signal filtering and denoising method based on generalized S transform containing variable factors
Hu et al. Compressive frequency hopping signal detection using spectral kurtosis and residual signals
Han et al. SAR wideband interference suppression method using second-order multisynchrosqueezing transform
Zhang et al. Seismic data denoising using double sparsity dictionary and alternating direction method of multipliers
Tao et al. Extraction and analysis of RFI signatures via deep convolutional RPCA
CN111398912B (en) Synthetic aperture radar interference suppression method based on tensor low-rank approximation
Mao et al. An Radio Frequency Interference Mitigation Approach for Spaceborne SAR System in Low SINR Condition
Zhang et al. Compressed sensing reconstruction of radar echo signal based on fractional fourier transform and improved fast iterative shrinkage-thresholding algorithm

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant