A Near-Field SAR Image Enhancement Method Based on Two-Dimensional Space-Variant Deconvolution

Abstract

The present invention discloses a near-field SAR image enhancement method based on two-dimensional spatially variant deconvolution. It calculates a rough response matrix, then constructs a rough convolutional degradation equation based on this matrix, and further uses the least absolute shrinkage and selection operator to solve the degradation equation to achieve deconvolution and obtain a rough enhanced image of the near-field SAR image. The K-means clustering algorithm is used to extract sub-rough enhanced images corresponding to different sub-class scattering points of the target, and the corresponding fine response matrices are calculated respectively, the corresponding fine spatially variant convolutional degradation equations are constructed, the degradation equations are solved to achieve spatially variant deconvolution, and the corresponding fine enhanced images are obtained. Finally, the overall SAR image enhancement result is obtained through the superposition of the fine enhanced images. Compared with other SAR image enhancement methods, this method can improve in terms of sidelobe suppression, clutter noise suppression, and target amplitude estimation retention during near-field imaging.

Description

Technical Field

[0001] The present invention belongs to the technical field of synthetic aperture radar, and particularly relates to the technical field of synthetic aperture radar (SAR) near-field imaging. Background Art

[0002] Synthetic Aperture Radar (SAR) is a radar system with all-weather and all-day imaging capabilities. It can image at any time during the day or night, in sunny or rainy and snowy weather, overcoming the disadvantages of optical and infrared systems that cannot image at night and in complex weather conditions. However, the traditional SAR has a long observation distance, a small relative viewing angle change with the observed target, and a small synthetic aperture angle, resulting in limited azimuth resolution of the image.

[0003] Different from traditional SAR, near-field SAR operates in the near-field region of the target. At this time, the relative viewing angle change with the observed target increases significantly, and the synthetic aperture angle is also large, so the azimuth resolution of the image is significantly improved, reaching centimeter-level resolution. In addition, the near-field SAR system usually can use stepped-frequency signals as the transmitted signal, so that the range resolution of the image can reach the same level as the azimuth resolution. Therefore, near-field SAR can achieve fine imaging of the target and has wide application value in many fields such as human body security inspection imaging, concealed object detection, building deformation detection, autonomous driving, and scattering characteristic measurement.

[0004] Currently, near-field SAR imaging usually adopts the back-projection algorithm based on matched filtering. Due to limited aperture sampling, there are inevitably sidelobes of the target in the near-field SAR imaging result. In addition, clutter from the ground and the surrounding environment also exists. These sidelobes and clutter have an adverse impact on the subsequent application of the image. For example, since they are easily misjudged as real targets, this will lead to a loss of accuracy in target detection and recognition. Therefore, it is necessary to enhance the imaging result through image enhancement methods to obtain a target image with lower target sidelobes and less clutter.

[0005] Currently, SAR image quality enhancement methods are mainly designed for far-field SAR. Therefore, the two-dimensional spatial variability of the response function of near-field SAR, which is different from that of far-field SAR, is not considered, resulting in limited applicability of the methods and limited enhancement of image quality. A representative method is the Iterative Adaptive Approach (IAA). Although some other methods, such as the CLEAN method, can be appropriately adjusted to consider the two-dimensional spatial variability. However, it still lacks additional consideration of the tightness of the target scattering distribution of near-field SAR, so it is prone to cumulative errors in forward and backward iterations, ultimately leading to deviations in target amplitude estimation. Generally speaking, the current enhancement methods applicable to near-field SAR images are still lacking. In order to better improve the quality of near-field SAR images, after fully considering the two-dimensional spatial variability of the near-field SAR system response function and the characteristics of the near-field SAR target scattering distribution, the present invention proposes a near-field SAR image enhancement method based on two-dimensional spatially variant deconvolution. Summary of the Invention

[0006] The present invention belongs to the field of synthetic aperture radar (SAR) imaging, and discloses a near-field SAR image enhancement method based on two-dimensional spatially variant deconvolution for enhancing the quality of near-field SAR images. The method first calculates a rough response matrix, then constructs a rough convolution degradation equation based on this matrix, and further uses the least absolute shrinkage and selection operator to solve the degradation equation to achieve deconvolution and obtain a rough enhanced image of the near-field SAR image. Further, the K-means clustering algorithm is used to extract sub-rough enhanced images corresponding to different sub-class scattering points of the target, and the corresponding fine response matrices are calculated respectively, and the corresponding fine spatially variant convolution degradation equations are constructed. Then, the least absolute shrinkage and selection operator is used to solve the degradation equation to achieve spatially variant deconvolution and obtain the corresponding fine enhanced images. Finally, the overall SAR image enhancement result is obtained by superimposing the fine enhanced images. Compared with other SAR image enhancement methods, this method can improve side lobe suppression, clutter and noise suppression, and target amplitude estimation retention in near-field imaging.

[0007] To facilitate the description of the content of the present invention, the following term definitions are made first:

[0008] Definition 1. Synthetic Aperture Radar

[0009] Synthetic Aperture Radar ((Synthetic Aperture Radar, SAR) is a high-resolution microwave imaging radar with the advantages of all-weather and all-day operation, and has been widely used in various fields, such as topographic mapping, guidance, environmental remote sensing, and resource exploration. For details of Synthetic Aperture Radar, see "Pi Yiming, Yang Jianyu, Fu Yusheng, Yang Xiaobo. Synthetic Aperture Radar Imaging Principle [M]. University of Electronic Science and Technology Press. 2007".

[0010] Definition 2. Stepped-frequency (SF) signal

[0011] The SF signal is a signal form that synthesizes a large equivalent bandwidth by transmitting a sub-pulse train with a stepped carrier frequency and has the ability of high-resolution range imaging. It can significantly reduce the instantaneous bandwidth of the system and the requirements for receiver hardware, and is widely used in the fields of microwave imaging and microwave measurement. For details of the stepped-frequency signal, see "Yang Ruliang. High-Resolution Microwave Imaging [M]. National Defense Industry Press. 2013".

[0012] Definition 3. Back-projection algorithm

[0013] The back-projection (BP) algorithm uses the trajectory information of the radar platform to obtain the range history between the radar platform and the scene pixels. Then, by traversing the range history, the matching echo data in the echo data is found, and then phase compensation and coherent accumulation are performed. Then, the complex-valued result is back-projected into the three-dimensional image space to complete the three-dimensional imaging process. For details of the back-projection algorithm, see "Shi Jun. Research on the Principles and Imaging Technologies of Bistatic SAR and Linear Array SAR [D]. PhD Thesis of University of Electronic Science and Technology of China. 2009".

[0014] Definition 4. Azimuth direction and range direction

[0015] The direction of the lateral movement of the radar platform is called the azimuth direction, and the direction perpendicular to the former is called the range direction.

[0016] Definition 5. Construction method of the N-point discrete Fourier transform matrix

[0017] The discrete Fourier transform matrix is an expression that represents the discrete Fourier transform in terms of matrix multiplication. The N-point discrete Fourier transform matrix can implement the N-point discrete Fourier transform. Specifically, the N-point discrete Fourier transform can be represented by an N×N matrix multiplication, that is, x F = Wx, where x is the original input signal, and x F is the output signal obtained after the discrete Fourier transform. The N-point discrete Fourier transform matrix W is composed of the following:

[0018]

[0019] where i is the imaginary unit. For details of the construction method of the N-point discrete Fourier transform matrix, see "He Zishu, Xia Wei. Modern Digital Signal Processing and Its Applications [M]. Beijing: Tsinghua University Press, 2009".

[0020] Definition 6. Matrix vectorization operator

[0021] The matrix vectorization operator vec(A) refers to the operation of stacking the input matrix A column by column to form a column vector. Specifically, the matrix vectorization operator arranges each column of a matrix A with dimensions m×n from left to right in a head-to-tail connection manner to form a column vector vec(A):

[0022] vec(A) = [a 1，1 , …, a m，1 , a 1，2 , …, a m，2 , …a 1，n , …, a m，n T

[0023] where a i，j represents A(i, j), the element in the i-th row and j-th column of matrix A; the superscript T represents the matrix transpose operation. For details of the matrix vectorization operator, see "Zhang Xianda. Matrix Analysis and Applications [M]. Tsinghua University Press Co., Ltd., 2004".

[0024] Definition 7. Vector-to-matrix diagonal operator

[0025] The vector-to-matrix diagonal operator diag(a) refers to the operation of generating a matrix A from the input vector a, and the diagonal elements of A are composed of the vector a. Specifically, the vector-to-matrix diagonal operator forms the diagonal elements of matrix A corresponding to the elements of a column vector a with dimensions n×1 in the order from top to bottom, and the other elements of the matrix are 0. diag(a) is expressed as:

[0026]

[0027] where a i , i = 1, 2, 3,...N represents the i-th element of a.

[0028] Definition 8. Standard least absolute value convergence and selection operator

[0029] The standard least absolute value convergence and selection operator l(y; A, γ) is used to solve the inverse problem based on the vector L1-norm regularization, and this problem can be expressed as the equation where γ is the regularization term weight parameter, is the square of the vector 2-norm. The operator is the solution of this equation, that is, l(y; A, γ) = x. For details of the standard least absolute value convergence and selection operator, see "Zhang Xianda. Matrix Analysis and Applications [M]. Tsinghua University Press Co., Ltd., 2004".

[0030] Definition 9. Traditional K-means clustering algorithm

[0031] ​The traditional K-means clustering algorithm is an iterative clustering analysis algorithm. The algorithm classifies the input data and outputs the classification results and the centers of each subclass. For the input data image matrix I, the algorithm output is the sub-images I of different subclass objects in the image matrix i and the class center c i =(c xi , c yi ), where i = 1, 2,..., K, K is the number of classes, c xi is the horizontal coordinate of the class center, and c yi is the vertical coordinate of the class center. The main steps of the algorithm include first randomly selecting K seed points as the initial class centers, then calculating the distances between the remaining points and each seed point, and assigning each point to the class center closest to it. This process will be repeated continuously until a certain termination condition is met. It is usually that the class centers no longer change. For details of the traditional K-means clustering algorithm, see "Wang Qian, Wang Cheng, Feng Zhenyuan, Ye Jinfeng. Research Review of K-means Clustering Algorithm [J]. Electronic Design Engineering, 2012, 20(07): 21-24".

[0032] Definition 10. Traditional Image Bilinear Interpolation Cropping and Rotation Algorithm

[0033] The traditional image bilinear interpolation rotation algorithm is a standard image rotation transformation algorithm. For the input image matrix I, the algorithm rotates it around the center of the image matrix according to the specified rotation angle θ i and the rotation direction (clockwise or counterclockwise), and calculates the missing values in the rotated image matrix using the bilinear interpolation method. Values outside the range of the original image matrix are cropped and discarded to keep the size of the output image matrix unchanged. For details of the traditional image bilinear interpolation cropping and rotation algorithm, see "Zhang Zheng, Wang Yanping, Xue Guixiang. Digital Image Processing and Machine Vision: Implementation with Visual C++ and Matlab [M]. Posts & Telecom Press, 2010."

[0034] Definition 11 Iterative Adaptive SAR Image Enhancement Method

[0035] The iterative adaptive SAR image enhancement method is a method with the core idea of maximum a posteriori probability estimation of two-dimensional spectrum. Specifically, this method uses the far-field space-invariant frequency response as the basic image model, and then takes the maximization of the target posterior probability as the optimization criterion to achieve high-resolution estimation of the target image spectrum. For details of the iterative adaptive SAR image enhancement method, see "Roberts W, Stoica P, Li J, et al. Iterative adaptive approaches to MIMO radar imaging[J]. IEEE Journal of Selected Topics in Signal Processing, 2010, 4(1): 5-20".

[0036] Definition of 12CLEAN SAR image enhancement method

[0037] The CLEAN SAR image enhancement method is an image enhancement method with the core idea of successively extracting target points and their system responses in the image. Specifically, the method iteratively extracts the point with the strongest energy in the image as the target point one by one, and at the same time subtracts the system response corresponding to the target point to achieve image enhancement. For details of the CLEAN SAR image enhancement method, see "Yun D J, Lee J I, Bae K U, et al. Improvement in computation time of 3-D scattering center extraction using the shooting and bouncing ray technique[J]. IEEE Transactions on Antennas and Propagation, 2017, 65(8): 4191-4199".

[0038] A near-field SAR image enhancement method based on two-dimensional space-variant deconvolution provided by the present invention includes the following steps:

[0039] Step 1. Initialize relevant parameters

[0040] Initialize the following parameters: the speed of light in air, denoted as c; take the natural exponential function, denoted as exp(·); the imaginary unit, denoted as j; pi, denoted as π; the system signal bandwidth, denoted as B; the wavelength of the system center frequency, denoted as λ; the distance from the center of the system imaging scene, denoted as r c ; the azimuth aperture length of the system, denoted as l a ; the number of pixels in the azimuth direction of the image, denoted as n a ; the number of pixels in the range direction of the image, denoted as n r; The pixel size in the azimuth direction of the image, denoted as i a ; The pixel size in the range direction of the image, denoted as i r ; The near-field SAR image matrix to be enhanced processed by the traditional back-projection algorithm, denoted as I bp ; The regularization weight parameter, denoted as γ; The number of clusters in the K-means clustering algorithm, denoted as K.

[0041] Step 2. Calculate the rough response matrix

[0042] For the speed of light c in air, the system signal bandwidth B, the wavelength λ of the system center frequency, and the distance r from the center of the system imaging scene in the initialization parameters in Step 1 c 、the azimuth aperture length l of the system a 、the number of pixels n in the range direction of the image r 、the number of pixels n in the azimuth direction of the image a 、the pixel size i in the range direction of the image r and the pixel size i in the azimuth direction of the image a , calculate the rough response function H of the system c :

[0043] Step 2.1. For the system signal bandwidth B and the speed of light c in air initialized in Step 1, use the formula d rc = c / (2B), calculate the range resolution of the center point of the imaging scene, denoted as d rc .

[0044] Step 2.2. For the wavelength λ of the system center frequency, the azimuth aperture length l of the system initialized in Step 1 a and the distance r from the center of the system imaging scene c , use the formula d ac = λr c / (2l a ), calculate the azimuth resolution of the center point of the imaging scene, denoted as d ac .

[0045] Step 2.3. For the number of pixels n in the range direction of the image initialized in Step 1 r , use the traditional method for constructing an N-point discrete Fourier transform matrix to construct an n r -point discrete Fourier transform matrix, denoted as W r .

[0046] Step 2.4 For the number of pixels n in the azimuth direction of the image initialized in Step 1 a , use the traditional method for constructing an N-point discrete Fourier transform matrix to construct an n a -point discrete Fourier transform matrix, denoted as W a .

[0047] Step 2.5 For the number of range pixels n of the image initialized in Step 1 r , the range pixel size i of the image r and the range resolution d of the center point of the imaging scene calculated in Step 2.1 rc , the size of the range support domain of the system's frequency-domain rough response function is calculated using the formula , denoted as n rfc , where represents rounding down.

[0048] Step 2.6 For the number of azimuth pixels n of the image initialized in Step 1 a , the azimuth pixel size i of the image a and the azimuth resolution d of the center point of the imaging scene calculated in Step 2.2 ac , the size of the azimuth support domain of the system's frequency-domain rough response function is calculated using the formula , denoted as m afc , where represents rounding down.

[0049] Step 2.7 For the number of range and azimuth pixels n of the image initialized in Step 1 r , the number of azimuth pixels n of the image a , the size of the range support domain n of the system's frequency-domain rough response function calculated in Step 2.5 rfc and the size of the azimuth support domain m of the system's frequency-domain rough response function calculated in Step 2.6 afc , the system's frequency-domain rough response matrix is calculated using the following formula, denoted as H cf .

[0050]

[0051] where 0 m×n represents a zero matrix of dimension m×n, and 1 m×n represents a matrix of all 1s of dimension m×n, represents rounding down.

[0052] Step 2.8 For the matrix W constructed in Step 2.3 r , the matrix W constructed in Step 2.4 a and the matrix H calculated in Step 2.7 cf , the system's rough response matrix is calculated using the formula , denoted as H c , where represents the Hadamard product of matrices, vec(·) is the matrix vectorization operator defined in Definition 6, and diag(·) is the vector matrix diagonalization operator defined in Definition 7.

[0053] Step 3. Construct a rough convolution degradation equation

[0054] For the system rough response matrix H calculated in Step 2 cf , the near-field SAR image matrix I to be enhanced initialized in Step 1 bp and the regularization weight parameter γ, denote I c as the rough enhanced image matrix, and construct the rough convolution degradation equation as follows:

[0055]

[0056] where represents the square of the vector 2-norm, |·| represents the vector 1-norm, and vec(·) is the matrix vectorization operator defined in Definition 6.

[0057] Step 4. Solve the rough convolution degradation equation

[0058] For the rough convolution degradation equation constructed in Step 3, use the standard least absolute value convergence and selection operator in Definition 8 to solve the equation to obtain the rough enhanced image matrix I c .

[0059] Step 5. Rough enhanced image clustering

[0060] For the number of classes K initialized in Step 1 and the rough enhanced image matrix I obtained in Step 4 c , use the traditional K-means clustering algorithm to obtain the rough enhanced image matrix of K different subclasses of target scattering points denoted as I i and the corresponding cluster centers denoted as c i =(c ai , c ri ), i = 1, 2,...K, where c ai is the azimuth coordinate of the center point of the i-th class of sub-scattering points, and c ri is the range coordinate of the center point of the i-th class of sub-scattering points.

[0061] Step 6. Calculate the fine response matrix

[0062] For the rough enhanced image matrix I of different subclasses of scattering points obtained in Step 5 i and the corresponding cluster centers c i =(c ai , c ri ), i = 1, 2,...K, the speed of light c in air, the system signal bandwidth B, the wavelength λ of the system center frequency, the center distance r of the system imaging scene c , the system azimuth aperture length l a , the number of pixels n in the range direction of the imager , the number of pixels in the azimuth direction of the image, n a , the pixel size in the range direction of the image, i r and the pixel size in the azimuth direction of the image, i a , calculate the system fine response function H i , i = 1, 2,... K:

[0063] Step 6.1. For the system center frequency wavelength λ and the system azimuth aperture length l initialized in Step 1 a and the cluster center c obtained in Step 5 i =(c ai , c ri ), use the formula to calculate the azimuth resolution of the center point of the image matrix of the i-th type of scattering point, denoted as d ai , i = 1, 2,... K.

[0064] Step 6.2 For the number of pixels in the azimuth direction of the image n initialized in Step 1 a , the pixel size in the azimuth direction of the image, i a and the azimuth resolution d of the center point of the image matrix of the i-th type of scattering point calculated in Step 6.1 ai , use the formula to calculate the azimuth support domain size of the frequency domain fine response function, denoted as n afi , where represents rounding down, i = 1, 2,... K.

[0065] Step 6.3 For the number of pixels in the range direction n initialized in Step 1 r , the number of pixels in the azimuth direction of the image, n a , the range support domain size n of the system frequency domain rough response function calculated in Step 2.5 rfc and the azimuth support domain size n of the frequency domain fine response function calculated in Step 6.2 afi , use the following formula to calculate the frequency domain fine response matrix, denoted as H f (i), i = 1, 2,... K.

[0066]

[0067] where 0 m×n represents a zero matrix of dimension m×n, and 1 m×n represents a matrix of all 1s of dimension m×n.

[0068] Step 6.4 For the cluster center c obtained in Step 5 i =(c ai , cri ) Calculate the fine response corrected rotation angle in the frequency domain using the formula and denote it as θ. i where arctan(·) is the arctangent function, and i = 1, 2,... K.

[0069] Step 6.5 For the fine response corrected rotation angle θ calculated in Step 6.4 i and the fine response matrix H obtained in Step 6.3 f (i), use the traditional bilinear interpolation rotation algorithm for images defined in Definition 10 to rotate the fine response matrix H f (i) counterclockwise by the corrected rotation angle θ i to obtain the fine response corrected matrix, denoted as H rf (i), where i = 1, 2,... K.

[0070] Step 6.6 For the matrix W constructed in Step 2.3 r , the matrix W constructed in Step 2.4 a and the matrix H calculated in Step 6.5 rf (i), use the formula to calculate the fine response matrix, denoted as H r (i), where represents the Hadamard product of matrices, vec(·) is the matrix vectorization operator defined in Definition 6, and diag(·) is the vector matrix diagonalization operator defined in Definition 7, and i = 1, 2,... K.

[0071] Step 7. Construct the fine convolution degradation equation

[0072] For the fine response matrix H calculated in Step 6 r (i), the rough enhanced image matrix I obtained in Step 5 i and the regularization weight parameter γ initialized in Step 1, denote Ir(i) as the fine enhanced image matrix of the i-th type of scattering point, and construct the following image rough-fine convolution equation:

[0073]

[0074] where represents the square of the vector 2-norm, |·| represents the vector 1-norm, vec(·) is the matrix vectorization operator defined in Definition 6, and i = 1, 2,... K.

[0075] Step 8. Solve the fine convolution degradation equation

[0076] For the refined convolutional degradation equation constructed in step 7, the standard least absolute value convergence and selection operator in Definition 8 is used to solve the convolutional degradation equation to achieve space-variant deconvolution, and the refined enhanced image matrix Ir(i) of the i-th type of scattering point is obtained, where i = 1, 2,... K.

[0077] Step 9. Superposition of refined enhanced images

[0078] For the refined enhanced image matrices Ir(i) of different subclasses of scattering points obtained in step 8, where i = 1, 2,... K, the formula is used to calculate the final enhanced image, denoted as I e 。

[0079] So far, the whole method ends.

[0080] The innovation and advantages of the present invention are as follows: A near-field SAR image enhancement method based on two-dimensional space-variant deconvolution is proposed. The advantage is that the degradation characteristics of near-field SAR images and the target scattering distribution characteristics are fully considered in the method design. The innovation lies in considering the two-dimensional space-variance of the system response and realizing image enhancement with the idea of space-variant deconvolution. Compared with other SAR image enhancement methods, in near-field imaging, it can improve in terms of sidelobe suppression, clutter and noise suppression, and target amplitude estimation retention. Description of the Drawings

[0081] Figure 1 is the flow chart of the present invention.

[0082] Figure 2 is the verification result of the simulation experiment of the present invention. Detailed Embodiment

[0083] The present invention is mainly verified by simulation experiments, and all steps and conclusions are verified to be correct on the mathematical calculation software Matlab2019b. The specific implementation steps are as follows:

[0084] Step 1. Initialize relevant parameters

[0085] Initialize the following parameters: The speed of light in air c = 3×10 8 m / s; Take the natural exponential function exp(·); Imaginary unit j; Pi π = 3.14; System signal bandwidth B = 2 GHz; Wavelength λ of the system center frequency = 0.03 m; Center distance r of the system imaging scene c = 25 m; Azimuth aperture length l of the system a = 5 m; Number of pixels n in the azimuth direction of the image a = 400; Number of pixels n in the range direction of the image r = 400; Pixel size i in the azimuth direction of the image a = 0.05 m; Pixel size i in the range direction of the imager = 0.05 m; the near-field SAR image matrix I to be enhanced processed by the back-projection algorithm defined in Definition 3 bp ; the regularization weight parameter γ = 1; the number of categories K of the K-means clustering algorithm = 9.

[0086] Step 2. Calculate the rough response matrix

[0087] According to the speed of light c = 3×10 8 m / s in air, system signal bandwidth B = 2 GHz, wavelength λ of the system center frequency = 0.03 m, center distance r of the system imaging scene c = 25 m, system azimuth aperture length l a = 5 m, number of pixels n in the azimuth direction of the image a = 400, number of pixels n in the range direction of the image r = 400, pixel size i in the range direction of the image r = 0.05 m, and pixel size i in the azimuth direction of the image a = 0.05 m, calculate the rough response function H of the system c :

[0088] Step 2.1. According to the system signal bandwidth B = 2 GHz and the speed of light c = 3×10 8 m / s initialized in Step 1, use the formula d rc = c / (2B) to calculate the range resolution d rc = 0.075 m of the center point of the imaging scene.

[0089] Step 2.2. According to the wavelength λ = 0.03 m of the system center frequency, system azimuth aperture length l a = 5 m, and system imaging scene center distance r c = 25 m initialized in Step 1, use the formula d ac = λr c / (2l a ) to calculate the azimuth resolution d ac = 0.075 m of the center point of the imaging scene.

[0090] Step 2.3. According to the number of pixels n in the range direction of the image r = 400 initialized in Step 1, construct a 400-point discrete Fourier transform matrix W according to the construction method of the N-point discrete Fourier transform matrix in Definition 5 r .

[0091] Step 2.4 According to the number of pixels n in the azimuth direction of the image a= 400. According to the method for constructing the N - point discrete Fourier transform matrix in Definition 5, construct a 400 - point discrete Fourier transform matrix W a .

[0092] Step 2.5. According to the number of pixels n in the range direction of the image initialized in Step 1 r = 400, the pixel size i in the range direction of the image r = 0.05 m, and the range - direction resolution d of the center point of the imaging scene calculated in Step 2.1 rc = 0.05 m, use the formula to calculate the size n of the range - direction support domain of the system's frequency - domain rough response function rfc = 266, where represents rounding down.

[0093] Step 2.6. According to the number of pixels n in the azimuth direction of the image initialized in Step 1 a = 400, the pixel size i in the azimuth direction of the image a = 0.05 m, and the azimuth - direction resolution d of the center point of the imaging scene calculated in Step 2.2 ac = 0.075 m, use the formula to calculate the size n of the azimuth - direction support domain of the system's frequency - domain rough response function afc = 266, where represents rounding down.

[0094] Step 2.7. According to the number of pixels n in the range direction and azimuth direction initialized in Step 1 r = 400, the number of pixels n in the azimuth direction of the image a = 400, the size n of the range - direction support domain of the system's frequency - domain rough response function calculated in Step 2.5 rfc = 266, and the size n of the azimuth - direction support domain of the system's frequency - domain rough response function calculated in Step 2.6 afc = 266, use the following formula to calculate the system's frequency - domain rough response matrix H cf .

[0095]

[0096] where 0 m×n represents a zero - matrix of dimension m×n, 1 m×n represents a one - matrix of dimension m×n,

[0097]

[0098]

[0099]

[0100]

[0101]

[0102] Step 2.8 According to the matrix W constructed in Step 2.3 r and the matrix W constructed in Step 2.4 a and the matrix H calculated in Step 2.7 cf , use the formula to calculate the system coarse response matrix H c , where represents the matrix Hadamard product, vec(·) is the matrix vectorization operator defined in Definition 6, and diag(·) is the vector matrix diagonalization operator defined in Definition 7.

[0103] Step 3. Construct the rough convolution degradation equation

[0104] According to the system coarse response matrix H calculated in Step 2 cf , the near-field SAR image matrix I to be enhanced initialized in Step 1 bp and the regularization weight parameter γ = 1, denote I c as the rough enhancement image matrix, and construct the rough convolution degradation equation as follows:

[0105]

[0106] where represents the square of the vector 2-norm, |·| represents the vector 1-norm, and vec(·) is the matrix vectorization operator defined in Definition 6.

[0107] Step 4. Solve the rough convolution degradation equation

[0108] According to the rough convolution degradation equation constructed in Step 3, use the standard least absolute value convergence and selection operator in Definition 8 to solve the equation, and obtain the rough enhancement image matrix I c .

[0109] Step 5. Rough enhancement map clustering

[0110] According to the number of categories K = 9 initialized in Step 1 and the rough enhancement image matrix I obtained in Step 4 c , use the traditional K-means clustering algorithm to obtain the rough enhancement image matrix I of 9 different subclasses of scattering points of the target iand the corresponding class centers c1(2.5, 32.5), c2(-2.5, 32.5), c3(7.5, 27.5), c4(7.5, 22.5), c5(2.5, 17.5), c6(-2.5, 17.5), c7(-7.5, 22.5), c8(-7.5, 27.5) and c9(0, 20).

[0111] Step 6. Calculate the fine response matrix

[0112] According to the rough enhanced image matrix I of different subclass scattering points obtained in Step 5 i and the corresponding class centers (c1(2.5, 32.5), c2(-2.5, 32.5), c3(7.5, 27.5), c4(7.5, 22.5), c5(2.5, 17.5), c6(-2.5, 17.5), c7(-7.5, 22.5), c8(-7.5, 27.5) and c9(0, 20)), the speed of light c = 3×10 8 m / s in air, the system signal bandwidth B = 2 GHz, the wavelength λ of the system center frequency = 0.03 m, the center distance r c of the system imaging scene = 25 m, the azimuth aperture length l a of the system = 5 m, the number of pixels n a in the azimuth direction of the image = 400, the number of pixels n r in the range direction of the image = 400, the pixel size i r in the range direction of the image = 0.05 m and the pixel size i a in the azimuth direction of the image = 0.05 m, calculate the system fine response function H i , f = 1, 2,... 9:

[0113] Step 6.1. According to the wavelength λ = 0.03 m of the system center frequency initialized in Step 1, the azimuth aperture length l a = 5 m of the system and the class centers c1(2.5, 32.5), c2(-2.5, 32.5), c3(7.5, 27.5), c4(7.5, 22.5), c5(2.5, 17.5), c6(-2.5, 17.5), c7(-7.5, 22.5), c8(-7.5, 27.5) and c9(0, 20) obtained in Step 5, use the formula and to calculate the azimuth resolution of the center points of the image matrices of the 1st to 9th subclass scattering points respectively as: d a1 = 0.0978 m, d a2 = 0.0978 m, d a3= 0.0855 m, d a4 = 0.0712 m, d a5 = 0.053 m, d a6 = 0.053 m, d a7 = 0.0712 m, d a8 = 0.0855 m and d a9 = 0.06 m.

[0114] Step 6.2 According to the number of pixels n in the azimuth direction of the image initialized in Step 1 a = 400, the pixel size i in the azimuth direction of the image a = 0.05 m and the azimuth resolution d of the center point of the image matrix calculated in Step 6.1 a1 = 0.0978 m, d a2 = 0.0978 m, d a3 = 0.0855 m, d a4 = 0.0712 m, d a5 = 0.053 m, d a6 = 0.053 m, d a7 = 0.0712 m, d a8 = 0.0855 m and d a9 = 0.06 m,

[0115] Adopt the formula and Calculate the size n of the azimuth support domain of the frequency-domain fine response function respectively af1 = 204, n af2 = 204, n af3 = 233, n af4 = 280, n af5 = 377, n af6 = 377, n af7 = 280, n af8 = 233 and n af9 = 333, where represents rounding down.

[0116] Step 6.3 According to the number of pixels n in the range and azimuth directions of the image initialized in Step 1 r = 400, the number of pixels n in the azimuth direction of the image a = 400, the size n of the range support domain of the system frequency-domain rough response function calculated in Step 2.5 rfc = 266 and the size n of the azimuth support domain of the image frequency-domain fine response function calculated in Step 6.2 af1 = 204, n af2 = 204, n af3 = 233, naf4 = 280, n af5 = 377, n af6 = 377, n af7 = 280, n af8 = 233 and n af9 = 333, calculate the fine response matrix H of the image frequency domain respectively using the following formula f (i), i = 1, 2,...9.

[0117]

[0118] where 0 m×n represents a zero matrix of dimension m×n, 1 m×n represents a matrix of all 1s of dimension m×n, n1(1) to n1(9) are respectively and n4(1) to n4(9) are respectively and n afc = 266, n rfc = 266.

[0119] Step 6.4 According to the cluster centers c1(2.5, 32.5), c2(-2.5, 32.5), c3(7.5, 27.5), c4(7.5, 22.5), c5(2.5, 17.5), c6(-2.5, 17.5), c7(-7.5, 22.5), c8(-7.5, 27.5) and c9(0, 20) obtained in Step 5, use the formula and calculate the fine response correction rotation angles θ1 = 4.3987°, θ2 = -4.3987°, θ3 = 15.2551°, θ4 = 18.4349°, θ6 = -8.1301°, θ7 = -18.4349°, θ ８ = -18.4349° and θ9 = 0°, where arctan(·) is the arctangent function.

[0120] Step 6.5 Correct the rotation angles (θ1 = 4.3987°, θ2 = -4.3987°, θ3 = 15.2551°, θ4 = 18.4349°, θ5 = 8.1301°, θ6 = -8.1301°, θ7 = -18.4349°, θ8 = -18.4349°, and θ9 = 0°) of the fine response in the image frequency domain according to the calculation in Step 6.4 and the fine response matrix H of the sub-image in the frequency domain obtained in Step 6.3 f (1)~H f (9) Use the traditional bilinear interpolation rotation algorithm for images defined in Definition 10 to the fine response matrix H of the sub-image in the frequency domain f (1)~H f (9) Rotate counterclockwise by the corrected rotation angles θ1 = 4.3987°, θ2 = -4.3987°, θ3 = 15.2551°, θ4 = 18.4349°, θ5 = 8.1301°, θ6 = -8.1301°, θ7 = -18.4349°, θ8 = -18.4349°, and θ9 = 0° respectively to obtain the corrected matrix H of the fine response in the image frequency domain rf (i), i = 1, 2,...K.

[0121] Step 6.6 According to the matrix W constructed in Step 2.3 r 、the matrix W constructed in Step 2.4 a and the matrix H calculated in Step 6.5 rf (i), use the formula to calculate the fine response matrix H of the image r (i), where represents the Hadamard product of matrices, vec(·) is the matrix vectorization operator defined in Definition 6, diag(·) is the vector matrix diagonal operator defined in Definition 7, and i = 1, 2,...9.

[0122] Step 7. Construct the fine convolution degradation equation

[0123] According to the fine response matrix H of the image calculated in Step 6 r (i), the rough enhanced image matrix I obtained in Step 5 i and the regularization weight parameter γ = 1 initialized in Step 1, denote Ir(i) as the fine enhanced image matrix of the i-th type of scattering point, and construct the following image rough-fine convolution equation

[0124]

[0125] where represents the square of the vector 2-norm, |·| represents the vector 1-norm, vec(·) is the matrix vectorization operator defined in Definition 6, and i = 1, 2,...9.

[0126] Step 8. Solve the fine convolution degradation equation

[0127] According to the fine convolution degradation equation of the image constructed in Step 7, use the standard least absolute value convergence and selection operator in Definition 8 to solve the convolution degradation equation, realize space-variant deconvolution, and obtain the fine enhanced image matrix Ir(i) of different subclass scattering points, where i = 1, 2,... 9.

[0128] Step 9. Superimpose the fine enhanced images

[0129] According to the fine enhanced image matrix Ir(i) of different subclass scattering points obtained in Step 8, where i = 1, 2,... 9, use the formula to calculate and obtain the final enhanced image I e .

[0130] So far, the whole method ends.

[0131] As Figure 2 shown, the computer simulation results show that the present invention fully considers the degradation characteristics of near-field SAR images and the target scattering distribution characteristics, realizes image enhancement with the idea of space-variant deconvolution, and has improvements in sidelobe suppression, clutter and noise suppression, and target amplitude estimation retention compared with other SAR image enhancement methods in near-field imaging.

Claims

1. A near-field SAR image enhancement method based on two-dimensional spatially variant deconvolution, characterized in that it Including the following steps: Step 1. Initialize relevant parameters Initialize the following parameters: the speed of light in air, denoted as c; the natural exponential function, denoted as exp(·); the imaginary unit, denoted as j; the pi, denoted as π; the system signal bandwidth, denoted as B; the center frequency wavelength of the system, denoted as λ; the center distance of the system imaging scene, denoted as r c ; the azimuth aperture length of the system, denoted as l a ; the number of azimuth pixels of the image, denoted as n a ; the number of range pixels of the image, denoted as n r ; the azimuth pixel size of the image, denoted as i a ; the range pixel size of the image, denoted as i r ; the near-field SAR image matrix to be enhanced processed by the back-projection algorithm, denoted as I bp ; the regularization weight parameter, denoted as γ; the number of classes of the K-means clustering algorithm, denoted as K; Step 2. Calculate the rough response matrix For the speed of light c in air, the system signal bandwidth B, the wavelength λ of the system center frequency, the center distance r of the system imaging scene in the initialization parameters in step 1 c 、the azimuth aperture length l of the system a 、the number of pixels n in the range direction of the image r 、the number of pixels n in the azimuth direction of the image a 、the pixel size i in the range direction of the image r and the pixel size i in the azimuth direction of the image a ,calculate the rough response function H of the system c : Step 2.

1. For the system signal bandwidth B initialized in Step 1 and the speed of light c in air, use the formula d rc = c / (2B) to calculate the range resolution of the center point of the imaging scene, denoted as d rc ; Step 2.

2. For the system center frequency wavelength λ, the system azimuth aperture length l a and the system imaging scene center distance r c initialized in Step 1, use the formula d ac = λr c / (2l a ) to calculate the azimuth resolution of the imaging scene center point, denoted as d ac ; Step 2.

3. For the number of pixels n in the range direction of the image initialized in Step 1 r , use the traditional method for constructing an N-point discrete Fourier transform matrix to construct an r n-point discrete Fourier transform matrix, denoted as W r ; Step 2.4 For the number n of azimuth pixels of the image initialized in step 1 a , a traditional method for constructing an N-point discrete Fourier transform matrix is used to construct an a n-point discrete Fourier transform matrix, denoted as W a ; Step 2.5 For the number of range pixels n of the image initialized in Step 1 r , the range pixel size i of the image r and the range resolution d of the center point of the imaging scene calculated in Step 2.1 rc , the formula is used to calculate the size of the range support domain of the system's frequency-domain rough response function, denoted as n rfc , where represents rounding down; Step 2.6 For the number of azimuth pixels n of the image initialized in Step 1 a , the azimuth pixel size i of the image a and the azimuth resolution d of the center point of the imaging scene calculated in Step 2.2 ac , the formula is used to calculate the azimuth support domain size of the system frequency domain rough response function, denoted as n afc , where represents rounding down; Step 2.7 For the number of range and azimuth pixels n initialized in Step 1 r , the number of azimuth pixels n of the image a , the size n of the range support domain of the system frequency-domain rough response function calculated in Step 2.5 rfc and the size n of the azimuth support domain of the system frequency-domain rough response function calculated in Step 2.6 afc , use the following formula to calculate the system frequency-domain rough response matrix, denoted as H cf ; where 0 m×n represents a zero matrix of dimension m×n, 1 m×n represents a matrix of all ones of dimension m×n, represents rounding down; Step 2.8 For the matrix W constructed in Step 2.3 r and the matrix W constructed in Step 2.4 a and the matrix H calculated in Step 2.7 cf , use the formula to calculate the system's rough response matrix, denoted as H c , where represents the matrix Hadamard product, vec(·) is the matrix vectorization operator, and diag(·) is the vector matrix diagonalization operator; Step 3. Construct the rough convolution degradation equation For the system rough response matrix H calculated in step 2 cf and the near-field SAR image matrix I to be enhanced initialized in step 1 bp and the regularization weight parameter γ, denote I c as the rough enhanced image matrix, and construct the rough convolution degradation equation as follows: where denotes the square of the vector 2-norm, |·| denotes the vector 1-norm, and vec(·) is the matrix vectorization operator; Step 4. Solve the rough convolution degradation equation For the rough convolutional degradation equation constructed in step 3, use the standard least absolute shrinkage and selection operator to solve the equation to obtain the rough enhanced image matrix I c ; Step 5. Rough enhanced image clustering For the number of categories K obtained by initialization in step 1 and the roughly enhanced image matrix I obtained in step 4 c , the traditional K-means clustering algorithm is used to obtain the roughly enhanced image matrix of different sub-class scattering points of K classes of targets, denoted as I i and the corresponding clustering centers, denoted as c i =(c ai , c ri ), i = 1, 2,...K, where c ai is the azimuth coordinate of the center point of the i-th type of sub-scattering point, and c ri is the range coordinate of the center point of the i-th type of sub-scattering point; Step 6. Calculate the fine response matrix For the rough enhanced image matrix I of different subclass scattering points obtained in step 5 i and the corresponding clustering center c i =(c ai , c ri ), i = 1, 2,...K, the speed of light c in air, system signal bandwidth B, system center frequency wavelength λ, system imaging scene center distance r c , system azimuth aperture length l a , number of pixels in the image range direction n r , number of pixels in the image azimuth direction n a , pixel size in the image range direction i r and pixel size in the image azimuth direction i a , calculate the system fine response function H i , i = 1, 2,...K: Step 6.

1. For the system center frequency wavelength λ and the system azimuth aperture length l initialized in Step 1 a and the clustering center c obtained in Step 5 i =(c ai , c ri ), use the formula to calculate the azimuth resolution of the image matrix center point of the i-th type of scattering point, denoted as d ai , i = 1, 2,... K; Step 6.2 For the number of azimuth pixels n of the image initialized in Step 1 a , the azimuth pixel size i of the image a and the azimuth resolution d of the center point of the image matrix of the i-th type of scattering point calculated in Step 6.1 ai , the formula is used to calculate the azimuth support domain size of the fine response function in the frequency domain, denoted as n afi , where represents rounding down, i = 1, 2,...K; Step 6.3 For the number of range pixels n initialized in Step 1 r , the number of azimuth pixels n of the image a , the size n of the range support domain of the system frequency-domain rough response function calculated in Step 2.5 rfc and the size n of the azimuth support domain of the frequency-domain fine response function calculated in Step 6.2 afi , use the following formula to calculate the frequency-domain fine response matrix, denoted as H f (i), where i = 1, 2,...K; where 0 m×n represents a zero matrix of dimension m×n, 1 m×n represents a matrix of all ones of dimension m×n; Step 6.4 For the clustering center c obtained in Step 5 i =(c ai , c ri ), the formula is used to calculate the fine response correction rotation angle in the frequency domain, denoted as θ i , where arctan(·) is the arctangent function, and i = 1, 2,... K; Step 6.5 Correct the rotation angle θ for the fine response in the frequency domain calculated in Step 6.4 i and the fine response matrix H in the frequency domain obtained in Step 6.3 f (i), use the traditional bilinear interpolation rotation algorithm for images to rotate the fine response matrix H f (i) counterclockwise by the corrected rotation angle θ i to obtain the corrected fine response matrix in the frequency domain, denoted as H rf (i), where i = 1, 2,...K; Step 6.6 For the matrix W constructed in Step 2.3 r and the matrix W constructed in Step 2.4 a and the matrix H calculated in Step 6.5 rf (i), use the formula to calculate the fine response matrix, denoted as H r (i), where represents the matrix Hadamard product, vec(·) is the matrix vectorization operator, diag(·) is the vector matrix diagonalization operator, and i = 1, 2,... K; Step 7. Construct the fine convolution degradation equation For the fine response matrix H obtained in step 6 r (i), the rough enhanced image matrix I obtained in step 5 i and the regularization weight parameter γ initialized in step 1, let Ir(i) be the fine enhanced image matrix of the i-th type of scattering point, and construct the following image rough-fine convolution equation: where denotes the square of the vector 2-norm, |·| denotes the vector 1-norm, vec(·) is the matrix vectorization operator, and i = 1, 2,... K; Step 8. Solve the fine convolution degradation equation For the fine convolution degradation equation constructed in Step 7, use the standard least absolute value convergence and selection operator to solve the convolution degradation equation, realize space-variant deconvolution, and obtain the fine enhanced image matrix Ir(i) of the i-th type of scattering point, where i = 1, 2,... K; Step 9. Superimpose the fine enhanced images For the fine-enhanced image matrix Ir(i) of different subclass scattering points obtained in step 8, where i = 1, 2,... K, the formula is used to calculate the final enhanced image, denoted as I e ; So far, the whole method ends.

Images