CS-SAR Imaging Method Based on Inverse Omega-K Algorithm

Hu Jingqiu Liu Falin Zhou Chongbin Li Bo Wang Dongjin

引用本文:
Citation:

一种新的基于Omega-K算法的稀疏场景压缩感知SAR成像方法

    通讯作者:
  • 基金项目:

    The National Natural Science Foundation of China 61431016

  • 中图分类号: TN957

CS-SAR Imaging Method Based on Inverse Omega-K Algorithm

    Author Bio: Hu Jingqiu (1992-) was born in Anhui, China. She is currently a master student in University of Science and Technology of China. She received her Bachelor degree in electronic engineering and information science from University of Science and Technology of China. Her current research interests include SAR imaging, compressed sensing, and sparse signal reconstruction;Liu Falin (1963-) was born in Xingtai, Hebei. He received the B.E. degree from Tsinghua University, Beijing, China, in 1985, and the M.E. and Ph.D. degrees both from University of Science and Technology of China (USTC), Hefei, China, in 1988 and 2004, respectively. Since 1988, he has been with the Department of Electronic Engineering and Information Science, USTC, where he is now a professor. His research interests include mm-wave devices, computational electromagnetics, microwave communications, and radar imaging. E-mail:liufl@ustc.edu.cn;Zhou Chongbin (1988-) was born in Luoyang, Henan Province, China. He received the B.Eng. degree from University of Science and Technology of China (USTC), Hefei, China, in 2011. He is currently working toward the Ph.D. degree at USTC. His research interests include compressive sensing, signal processing, and radar imaging;Li Bo (1991-) was born in Baoji, Shaanxi Province, China. He received the Bachelor degree from Xidian University, Xi'an, China, in 2013. He is currently working toward the Ph.D. degree at USTC. His research interests include compressive sensing, signal processing, and radar imaging;Wang Dongjin (1955-) was born in Huainan, Anhui Province. He received the Bachelor's degree from University of Science and Technology of China (USTC), Hefei, China, in 1982, and the M.E. degree from Nanjing Institute of Electronic Technology, Nanjing, China, in 1985. He has been a full professor since 1998. Prof. Wang had been the vice president of USTC since 2003 and is now the director of the Key Laboratory of Electromagnetic Space Information, Chinese Academy of Sciences. His research interests include electromagnetic theory, mm-wave communications and radar systems, and applications.
    Corresponding author: Liu Falin. E-mail:liufl@ustc.edu.cn
  • Fund Project: The National Natural Science Foundation of China 61431016

    CLC number: TN957

  • 摘要: 很多文献已经证明压缩感知应用在SAR成像中的有效性.现有的CS-SAR成像算法非常耗时, 尤其是对于高分辨率的图像来说更甚.该文针对稀疏场景提出了一种基于omega-K算法, 精确且高效的CS-SAR成像算法——CS-OKA算法.我们首先推导出了omega-K算法的逆算子, 可不通过发射信号和场景的卷积来直接得到回波信号.在此基础上我们将SAR成像问题建立为一个稀疏优化问题, 并用迭代阈值的方法来求解.仿真结果表明, 当场景稀疏时该文的方法可以在远低于Nyquist采样率的前提下有效地恢复出原始场景, 并且时耗和存储量都显著降低.
  • Figure 1.  Reconstruction results with full sampled data but different SNRs. The left column is recovered by omega-K algorithm, and the right column is recovered by CS-OKA. From top to bottom, SNRs are 20, 10, and 5 dB, respectively

    Figure 2.  Contours of magnitude of lower left target in Fig. 1. The left column is recovered by omega-K algorithm, and the right column is recovered by CS-OKA. From top to bottom, SNRs are 20, 10, and 5 dB, respectively

    Figure 3.  Reconstruction results with 50% samples in both range and azimuth. From left to right, SNRs are 20, 10, and 5 dB, respectively

    Figure 4.  Reconstruction results with 10% samples in both range and azimuth. From left to right, SNRs are 20, 10, and 5 dB, respectively

    Figure 5.  Reconstruction results of two close targets

    Algorithm 1: Iterative thresholding algorithm for proposed CS-SAR imaging
    Input: SAR raw echoes SS,
        omega-K algorithm M and inverse omega-K algorithm T,
        sampling operator Θτ and Θη
    Initial: G(0), λ, μ, and max iteration Imax
        1: for i = 1 to Imax do
        2: residue: $ {{\mathit{\boldsymbol{R}}}^{\left( i-1 \right)}}={{\mathit{\boldsymbol{S}}}_{S}}-{\mathit{\Theta }_\tau }T\left( {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}} \right){{\mathit{\Theta }}_{\eta }}$
        3: omega-K on residue: $\Delta {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}}=M\left( \mathit{\Theta } _{\tau }^{\rm{T}}{{\mathit{\boldsymbol{R}}}^{\left( i-1 \right)}}\mathit{\Theta } _{\eta }^{\rm{T}} \right)$
        4: Thresholding: ${{\mathit{\boldsymbol{G}}}^{\left( i \right)}}={{E}_{1, \lambda \mu }}\left( {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}}+\mu \Delta {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}} \right)$
        5: end for
    Output: The recovery image G*=G(i)
    下载: 导出CSV

    Table 1.  Parameters used in the simulation

    Parameter Value
    Pulse duration (μs) 1.33
    Bandwidth in range (MHz) 150
    Carrier frequency (MHz) 600
    Sampling rate (MHz) 225
    Slant range of scene center (m) 1200
    Length of synthetic aperture (m) 300
    Pulse repeat frequency (Hz) 150 000 000
    Radar velocity in azimuth (m/s) 150
    下载: 导出CSV
  • [1] Cumming I G and Wong F H. Digital Processing of Synthetic Aperture Radar Data:Algorithms and Implementation[M]. Norwood, MA, USA, Artech House, 2004:225-367.
    [2] Donoho D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4):1289-1306. doi: 10.1109/TIT.2006.871582
    [3] Baraniuk R G. Compressive sensing[J]. IEEE Signal Processing Magazine, 2007, 24(4):118-120. doi: 10.1109/MSP.2007.4286571
    [4] Candes E J, Romberg J K, and Tao T. Stable signal recovery from incomplete and inaccurate measurements[J]. Communications on Pure and Applied Mathematics, 2006, 59(8):1207-1223. doi: 10.1002/(ISSN)1097-0312
    [5] Davenport M A, Duarte M F, Eldar Y C, et al.. Compressed Sensing:Theory and Applications[M]. Cambridge, U.K., Cambridge University Press, 2012:1-55.
    [6] 吴一戎, 洪文, 张冰尘, 等.稀疏微波成像研究进展[J].雷达学报, 2014, 3(4):383-395.Wu Yi-rong, Hong Wen, Zhang Bing-chen, et al.. Current developments of sparse microwave imaging[J]. Journal of Radars, 2014, 3(4):383-395.
    [7] Baraniuk R and Steeghs P. Compressive radar imaging[C]. IEEE Radar Conference, Waltham, MA, USA, 2007: 128-133.
    [8] Alonso Mariví Tello, López-Dekker Paco, and Mallorquí Jordi J. A novel strategy for radar imaging based on compressive sensing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(12):4285-4295. doi: 10.1109/TGRS.2010.2051231
    [9] Yang Jungang, Thompson J, Huang Xiaotao, et al.. Segmented reconstruction for compressed sensing SAR imaging[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(7):4214-4225. doi: 10.1109/TGRS.2012.2227060
    [10] Fang Jian, Xu Zongben, Zhang Bingchen, et al.. Fast compressed sensing SAR imaging based on approximated observation[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(1):352-363. doi: 10.1109/JSTARS.2013.2263309
    [11] Dong Xiao and Zhang Yunhua. A novel compressive sensing algorithm for SAR imaging[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(2):708-720. doi: 10.1109/JSTARS.2013.2291578
    [12] Bu Hongxia, Tao Ran, Bai Xia, et al.. A novel SAR imaging algorithm based on compressed sensing[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12(5):1003-1007. doi: 10.1109/LGRS.2014.2372319
    [13] Yang Dong, Liao Guisheng, Zhu Shengqi, et al.. SAR imaging with undersampled data via matrix completion[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(9):1539-1543. doi: 10.1109/LGRS.2014.2300170
    [14] Dong Xiao and Zhang Yunhua. A MAP approach for 1-bit compressive sensing in synthetic aperture radar imaging[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12(6):1237-1241. doi: 10.1109/LGRS.2015.2390623
    [15] Zhang Siqian, Zhu Yutao, Dong Ganggang, et al.. Truncated SVD-based compressive sensing for downward-looking three-dimensional SAR imaging with uniform nonuniform linear array[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12(9):1853-1857. doi: 10.1109/LGRS.2015.2431254
    [16] Bi Hui, Jiang Chenglong, Zhang Bingchen, et al.. Radar change imaging with undersampled data based on matrix completion and Bayesian compressive sensing[J].IEEE Geoscience and Remote Sensing Letters, 2015, 12(7):1546-1550. doi: 10.1109/LGRS.2015.2412677
    [17] Khwaja Ahmed Shaharyar, Ferro-Famil Laurent, and Pottier Eric. Efficient SAR raw data generation for anisotropic urban scenes based on inverse processing[J]. IEEE Geoscience and Remote Sensing Letters, 2009, 6(4):757-761. doi: 10.1109/LGRS.2009.2024559
    [18] Blumensath T and Davies M E. Iterative hard thresholding for compressed sensing[J]. Applied and Computational Harmonic Analysis, 2009, 27(3):265-274. doi: 10.1016/j.acha.2009.04.002
    [19] Daubechies Ingrid, Defrise Michel, and De Mol Christine. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[J]. Communications on Pure and Applied Mathematics, 2004, 57(11):1413-1457. doi: 10.1002/(ISSN)1097-0312
  • [1] 洪文 . 圆迹SAR 成像技术研究进展. 雷达学报, 2012, 1(2): 124-135. doi: 10.3724/SP.J.1300.2012.20046
    [2] Alejandro C. Frery . Stochastic Contrast Measures for SAR Data: A Survey (in English). 雷达学报, 2019, 8(6): 758-781. doi: 10.12000/JR19108
    [3] 顾福飞张群杨秋霍文俊王敏 . 基于NCS算子的大斜视SAR压缩感知成像方法. 雷达学报, 2016, 5(1): 16-24. doi: 10.12000/JR15035
    [4] 韦顺军田博坤张晓玲师君 . 基于半正定规划的压缩感知线阵三维SAR自聚焦成像算法. 雷达学报, 2018, 7(6): 664-675. doi: 10.12000/JR17103
    [5] 闫敏韦顺军田博坤张晓玲师君 . 基于稀疏贝叶斯正则化的阵列SAR高分辨三维成像算法. 雷达学报, 2018, 7(6): 705-716. doi: 10.12000/JR18067
    [6] 李烈辰李道京黄平平 . 基于变换域稀疏压缩感知的艇载稀疏阵列天线雷达实孔径成像. 雷达学报, 2016, 5(1): 109-117. doi: 10.12000/JR14159
    [7] 刘振魏玺章黎湘 . 一种新的随机PRI脉冲多普勒雷达无模糊MTD算法. 雷达学报, 2012, 1(1): 28-35. doi: 10.3724/SP.J.1300.2013.10063
    [8] 赵娟白霞 . 一种适用于TDOMP算法的测量矩阵优化方法. 雷达学报, 2016, 5(1): 8-15. doi: 10.12000/JR15131
    [9] 杨军张群罗迎邓冬虎 . 基于压缩感知的认知雷达多目标跟踪方法. 雷达学报, 2016, 5(1): 90-98. doi: 10.12000/JR14107
    [10] 朱敏慧 . SAR 的海洋动力探测研究及应用浅析. 雷达学报, 2012, 1(4): 342-352. doi: 10.3724/SP.J.1300.2012.20088
    [11] 张问一胡东辉丁赤飚 . 基于FABEMD 和Goldstein 滤波器的SAR 舰船尾迹图像增强方法. 雷达学报, 2012, 1(4): 426-435. doi: 10.3724/SP.J.1300.2012.20059
    [12] 王志豪李刚蒋骁 . 基于光学和SAR遥感图像融合的洪灾区域检测方法. 雷达学报, 2020, 9(): 1-15. doi: 10.12000/JR19095
    [13] 洪文林赟谭维贤王彦平向茂生 . 地球同步轨道圆迹SAR研究. 雷达学报, 2015, 4(3): 241-253. doi: 10.12000/JR15062
    [14] 刘宁赵博黄磊 . 单通道SAR抗欺骗干扰方法. 雷达学报, 2019, 8(1): 73-81. doi: 10.12000/JR18072
    [15] 周伟刘永祥黎湘凌永顺 . MIMO-SAR 技术发展概况及应用浅析. 雷达学报, 2014, 3(1): 10-18. doi: 10.3724/SP.J.1300.2013.13074
    [16] 杨震杨汝良 . HJ-1-C 卫星SAR 系统的内定标. 雷达学报, 2014, 3(3): 314-319. doi: 10.3724/SP.J.1300.2014.14028
    [17] 张增辉郁文贤 . 稀疏微波SAR图像特征分析与目标检测研究. 雷达学报, 2016, 5(1): 42-56. doi: 10.12000/JR15097
    [18] 贾丽贾鑫许小剑何永华 . 机场场景SAR原始数据模拟. 雷达学报, 2014, 3(5): 565-573. doi: 10.3724/SP.J.1300.2014.14071
    [19] 毛永飞汪小洁向茂生 . 机载干涉SAR 区域网三维定位算法. 雷达学报, 2013, 2(1): 60-67. doi: 10.3724/SP.J.1300.2012.20107
    [20] 陈功伯李勇陶满意 . 基于信号的环视SAR 成像参数估计方法. 雷达学报, 2013, 2(2): 203-209. doi: 10.3724/SP.J.1300.2013.20073
  • 加载中
图(5)表(2)
计量
  • 文章访问数:  1276
  • HTML浏览量:  59
  • PDF下载量:  791
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-01-30
  • 录用日期:  2016-03-29
  • 网络出版日期:  2016-06-03
  • 刊出日期:  2017-02-28

CS-SAR Imaging Method Based on Inverse Omega-K Algorithm

    通讯作者:
    作者简介: Hu Jingqiu (1992-) was born in Anhui, China. She is currently a master student in University of Science and Technology of China. She received her Bachelor degree in electronic engineering and information science from University of Science and Technology of China. Her current research interests include SAR imaging, compressed sensing, and sparse signal reconstruction;Liu Falin (1963-) was born in Xingtai, Hebei. He received the B.E. degree from Tsinghua University, Beijing, China, in 1985, and the M.E. and Ph.D. degrees both from University of Science and Technology of China (USTC), Hefei, China, in 1988 and 2004, respectively. Since 1988, he has been with the Department of Electronic Engineering and Information Science, USTC, where he is now a professor. His research interests include mm-wave devices, computational electromagnetics, microwave communications, and radar imaging. E-mail:liufl@ustc.edu.cn;Zhou Chongbin (1988-) was born in Luoyang, Henan Province, China. He received the B.Eng. degree from University of Science and Technology of China (USTC), Hefei, China, in 2011. He is currently working toward the Ph.D. degree at USTC. His research interests include compressive sensing, signal processing, and radar imaging;Li Bo (1991-) was born in Baoji, Shaanxi Province, China. He received the Bachelor degree from Xidian University, Xi'an, China, in 2013. He is currently working toward the Ph.D. degree at USTC. His research interests include compressive sensing, signal processing, and radar imaging;Wang Dongjin (1955-) was born in Huainan, Anhui Province. He received the Bachelor's degree from University of Science and Technology of China (USTC), Hefei, China, in 1982, and the M.E. degree from Nanjing Institute of Electronic Technology, Nanjing, China, in 1985. He has been a full professor since 1998. Prof. Wang had been the vice president of USTC since 2003 and is now the director of the Key Laboratory of Electromagnetic Space Information, Chinese Academy of Sciences. His research interests include electromagnetic theory, mm-wave communications and radar systems, and applications
  • ①. Department of EEIS, University of Science and Technology of China, Hefei 230027, China
  • ②. Key Laboratory of Electromagnetic Space Information, Chinese Academy of Sciences, Hefei 230027, China
基金项目:  The National Natural Science Foundation of China 61431016

摘要: 很多文献已经证明压缩感知应用在SAR成像中的有效性.现有的CS-SAR成像算法非常耗时, 尤其是对于高分辨率的图像来说更甚.该文针对稀疏场景提出了一种基于omega-K算法, 精确且高效的CS-SAR成像算法——CS-OKA算法.我们首先推导出了omega-K算法的逆算子, 可不通过发射信号和场景的卷积来直接得到回波信号.在此基础上我们将SAR成像问题建立为一个稀疏优化问题, 并用迭代阈值的方法来求解.仿真结果表明, 当场景稀疏时该文的方法可以在远低于Nyquist采样率的前提下有效地恢复出原始场景, 并且时耗和存储量都显著降低.

English Abstract

    • Synthetic Aperture Radar (SAR) is widely used to obtain images under all weather conditions at any time of a day[1]. With the growing demand for SAR images and the development of SAR imaging technology, the image resolution has increased significantly. Traditionally, the scene is reconstructed by Matched Filter (MF)-based focusing algorithms. However, these methods are efficient above Nyquist sample rate, which would bring a lot of burden to current SAR system hardware with too much data.

      The recent theory of Compressed Sensing (CS) appears to be a powerful method to solve the problem above. CS theory proves that a sparse signal can be reconstructed with much fewer measurements than Nyquist sample rate needs[2-5]. In the past several years, CS has been used in many radar applications and shown very good results. The theory was applied to eliminate the need of the pulse compression matched filter and reduce the required receiver's analog-to-digital conversion bandwidth[6, 7]. In Ref. [8], after range compression, CS was then used on azimuth, but the signal was operated in One-Dimensional (1-D) form. For Two-Dimensional (2-D) operations, some work has been done in Refs. [9-12]. Many different algorithms have been used to solve SAR imaging problem, including matrix completion[13], 1-bit compressive sensing[14], truncated Singular Value Decomposition (SVD)[15], and Bayesian compressive sensing[16].

      Traditional MF-based methods have relatively lower computational complexity and memory cost compared to CS-SAR methods, but they suffer from high sidelobes. Range Doppler Algorithm (RDA), a traditional MF-based method of SAR imaging, has some drawback in correcting the coupling between range and azimuth. Another traditional MF-based method, the Chirp Scaling Algorithm (CSA), uses an approximated equation of signal when operating range IFFT. The omega-K algorithm adopts a special operation to correct the coupling of range and azimuth to overcome the drawback of RDA. Besides, omega-K algorithm uses a precise expression of signal compared to CSA[1].

      In this study, we propose a new CS-SAR method for sparse scene based on omega-K algorithm, which is named as Compressed Sensing omega-K Algorithm (CS-OKA). This new method could reduce computational complexity and memory cost while achieving high resolution and eliminating sidelobes with reduced sampling rate. First, we introduce the procedure of omega-K algorithm and represent it in the form of multiplications of matrices. Then we try to derive its inverse in order to get raw signal from SAR image. Similar practice has been explored before finding the inverse of MF-based method[17]. Finally, we combine the inverse of omega-K algorithm and CS framework to formulate a method from which undersampled data can recover the original SAR image fast and precisely. Different from previous CS-SAR algorithm, we get the echoes from the inverse of omega-K algorithm, not from the original observation. This results in the reduction of computational complexity. The whole procedure is performed two-dimensionally. To get the expected sparse scene, we use Iterative Thresholding Algorithm[18, 19] (ITA) to solve the sparse problem. The CS-OKA is much faster than previous CS-SAR algorithms since we substitute the large matrix-vector multiplication with omega-K algorithm and its inverse which can be realized efficiently. Also, memory cost is significantly reduced. Simulations show that CS-OKA can reconstruct the sparse SAR images effectively below Nyquist rate. Refs. [10] and [11] use similar process but with RDA and CSA, respectively. As mentioned above, omega-K has some advantages compared to RDA and CSA, and that is why we chose omega-K in our algorithm.

      This paper is organized as follows: In Section 2, the CS theory is introduced. In Section 3, the inverse of omega-K algorithm is described. In Section 4, our new CS-SAR method CS-OKA is presented. In Section 5, simulation results are shown to evaluate the performance of CS-OKA. Finally, Section 6 concludes the paper.

    • Suppose x$\in {{\mathbb{C}}^{N}}$ is the signal to be estimated. We call it K-sparse if its nonzero elements are no more than K[2, 5]. Let $x\in {{\mathbb{C}}^{M}}$ be the collected data, and matrix $\mathit{\Theta }\in {{\mathbb{R}}^{M\times N}}(M < N)$ be the observation matrix. According to CS theory, it is possible to recover x if x is sparse:

      $$ \mathit{\boldsymbol{y}}=\mathit{\Theta }\mathit{\boldsymbol{x}}+\mathit{\boldsymbol{n}} $$

      where n $\in {{\mathbb{C}}^{M}} $ represents the additive noise. For those signals which are not sparse themselves, we sometimes can find a transform or basis so that their coefficients under the basis are sparse:

      $$ \mathit{\boldsymbol{f}}=\mathit{\Psi }\mathit{\boldsymbol{x}} $$

      f is the non-sparse signal to be estimated; Ψ is the basis; x and is the sparse coefficient vector under the basis. In this case, the expression would be

      $$ \mathit{\boldsymbol{y=}}\mathit{\Theta \Psi }\mathit{\boldsymbol{x}}+\mathit{\boldsymbol{n=Ax}}+\mathit{\boldsymbol{n}} $$

      where A is sensing matrix.

      To recover the sparse signal x, CS theory puts forward some requirements on the sensing matrix. One sufficient condition is the Restricted Isometry Property (RIP). If there exists a constant ${{\delta }_{K}}\in \left( 0, 1 \right) $ such that

      $$ \left( 1-{{\delta }_{K}} \right)\left\| \mathit{\boldsymbol{x}} \right\|_{2}^{2}\le \left\| \mathit{\boldsymbol{Ax}} \right\|_{2}^{2}\le \left( 1+{{\delta }_{K}} \right)\left\| \mathit{\boldsymbol{x}} \right\|_{2}^{2} $$

      holds for all K-sparse signal x, the matrix A is said to satisfy the RIP of order K. When the RIP constant satisfies some conditions, the recovery can be achieved by solving a convex optimization problem:

      $$ \mathop {\min }\limits_\mathit{\boldsymbol{x}} {\mkern 1mu} {\left\| \mathit{\boldsymbol{x}} \right\|_1}\;\;{\rm{subject}}\;\;{\rm{to}}{\left\| {\mathit{\boldsymbol{Ax}} - \mathit{\boldsymbol{y}}} \right\|_{\rm{2}}} \leqslant \varepsilon $$

      where ‖·‖1 denotes the l1 norm and e bounds the energy of noise (${{\left\| \mathit{\boldsymbol{n}} \right\|}_{2}}\le \varepsilon $). An equivalent regularization scheme has the form of

      $$ \mathop {\min }\limits_\mathit{\boldsymbol{x}} {\mkern 1mu} \left\| {\mathit{\boldsymbol{Ax}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + \lambda {\left\| \mathit{\boldsymbol{x}} \right\|_{\rm{1}}} $$

      where λ is a regularization parameter. ITA is an general algorithm to solve the optimization problem. It uses the iteration

      $$ {\mathit{\boldsymbol{x}}^{\left( {i + 1} \right)}} = {E_{q,\lambda \mu }}\left( {{\mathit{\boldsymbol{x}}^{\left( i \right)}} + \mu {\mathit{\boldsymbol{A}}^{\rm{H}}}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{x}}^{\left( {\bf{i}} \right)}}} \right)} \right) $$

      to implement the procedure. In Eq. (7), μ is the step size that controls the convergence. Let σ=λμ, then Eq, σ is the thresholding operator which is defined in terms of components by

      $$ {{E}_{q, \sigma }}\left( \mathit{\boldsymbol{x}} \right)={{\left( {{\rm{e}}_{q, \sigma }}\left( {{\mathit{\boldsymbol{x}}}_{1}} \right), {{\rm{e}}_{q, \sigma }}\left( {{\mathit{\boldsymbol{x}}}_{2}} \right), \cdots, {{\rm{e}}_{q, \sigma }}\left( {{\mathit{\boldsymbol{x}}}_{n}} \right) \right)}^{\rm{T}}} $$

      where eq, σ can be analytically specified when q=0, 1/2, 2/3, 1. If q=1, then it corresponds to soft-thresholding, and Eq, σ would have the following expression:

      $$ {{\rm{e}}_{1, \sigma }}\left( \mathit{\boldsymbol{x}} \right)=\left\{ \begin{align} & \left( \mathit{\boldsymbol{x/}}{{\left\| \mathit{\boldsymbol{x}} \right\|}_{2}} \right)\left( \left| \mathit{\boldsymbol{x}} \right|-\sigma \right), \left| \mathit{\boldsymbol{x}} \right|\ge \sigma \\ & 0, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \rm{otherwise} \\ \end{align} \right. $$

      Solving Eq. (5) or Eq. (6) with ITA is very time consuming, especially for large scale SAR imaging. Because it involves large matrix-vector multiplications. In Section 3, we will derive a new way to form the sensing matrix, which would accelerate the process significantly.

    • As mentioned in Section 2, when using Eq. (7) to solve the optimization problem, we have to conduct multiplications of very large scale matrices. However, if we can substitute A and AH with much simpler operators, we will avoid such time consuming multiplications. In this study, we derive such operators based on omega-K algorithm. Detailed information is further discussed later.

      The omega-K algorithm[1] has a sufficient accuracy with a moderate computational cost. In this section, after a brief review of omega-K algorithm, we will derive the inverse of omega-K algorithm.

    • Let s0(τ, η) be the received data, where τ and η represent fast range time and slow azimuth time, respectively:

      $$ \begin{align} &{{s}_{0}}\left( \tau, \eta \right)={{\text{A}}_{0}}{{w}_{\text{r}}}\left( \tau-2R\left( \eta \right)/c \right){{w}_{\text{a}}}\left( \eta-{{\eta }_{c}} \right) \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdot \text{exp}\left\{-\text{j4 }\!\!\pi\!\!\text{ }{{f}_{0}}R\left( \eta \right)/c \right\} \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdot \text{exp}\left\{ \text{j }\!\!\pi\!\!\text{ }{{K}_{\text{r}}}{{\left( \tau -2R\left( \eta \right)/c \right)}^{2}} \right\} \\ \end{align} $$

      where A0 is a complex constant; wr(τ)=rect (τ/Tr) and Tr is the transmitted signal duration of each pulse; wa(η)=rect (η/Tη) and Tη is the integration time; R(η) is the slant range; f0 is the carrier frequency; Kr is chirp rate. To perform omega-K algorithm, we need to go through four steps.

      The first step is 2-D Fourier transform. Let ${{\cal F}_\tau }\left\{ \cdot \right\}$ and ${{\cal F}_\eta }\left\{ \cdot \right\}$ be the Fourier transform operator for range and azimuth, respectively. The transform would be:

      $$ {S_{2{\rm{df}}}}\left( {{f_\tau }, {f_\eta }} \right) = {{\cal F}_\eta }\left\{ {{{\cal F}_\tau }\left\{ {{s_0}} \right\}} \right\} $$

      The second step is multiplying Eq. (11) with a reference function ${H_{{\rm{ref}}}}\left( {{f_\tau }, {f_\eta }} \right)$. This step is to compensate the phase for the target in the scene center:

      $$ {{H}_{\text{ref}}}\left( {{f}_{\tau }}, {{f}_{\eta }} \right)=\text{exp}\left[\text{j}\frac{4\text{ }\!\!\pi\!\!\text{ }{{R}_{\text{ref}}}}{c}\sqrt{{{\left( {{f}_{0}}+{{f}_{\tau }} \right)}^{2}}-\frac{{{c}^{2}}f_{\eta }^{2}}{4V_{{{\text{r}}_{\text{ref}}}}^{2}}+\text{j}\frac{\text{ }\!\!\pi\!\!\text{ }f_{\tau }^{2}}{{{K}_{\text{r}}}}} \right] $$

      where Rref is the slant range of reference target, and ${V_{{{\rm{r}}_{{\rm{ref}}}}}}$ is radar velocity at reference target position. We usually assume that radar velocity remains the same (Vr) during the procedure. After reference function multiplication, Eq. (11) becomes

      $$ {S_{{\rm{RFM}}}}\left( {{f_\tau }, {f_\eta }} \right) = {S_{2{\rm{df}}}}\left( {{f_\tau }, {f_\eta }} \right){H_{{\rm{ref}}}} $$

      The phase of ${S_{{\rm{RFM}}}}\left( {{f_\tau }, {f_\eta }} \right)$ is

      $$ {{\theta }_{\text{RFM}}}\approx-\frac{4\text{ }\!\!\pi\!\!\text{ }\left( {{R}_{0}}-{{R}_{\text{ref}}} \right)}{c}\sqrt{{{\left( {{f}_{0}}+{{f}_{\tau }} \right)}^{2}}-\frac{{{c}^{2}}{{f}_{\eta }}^{2}}{4V_{{{\text{r}}_{\text{ref}}}}^{2}}} $$

      The third step is Stolt interpolation. If the square root term in Eq. (14) is substituted by

      $$ \sqrt{{{\left( {{f}_{0}}+{{f}_{\tau }} \right)}^{2}}-\frac{{{c}^{2}}{{f}_{\eta }}{{}^{2}}}{4V_{{{\text{r}}_{\text{ref}}}}^{2}}}={{f}_{0}}+f_{\tau }^{'} $$

      then the coupling effect between range and azimuth is thus removed. After the substitution, the frequency spectrum of ${S_{{\rm{RFM}}}}\left( {{f_\tau }, {f_\eta }} \right) $ is changed. For the convenience of the last step, a truncated sinc-kernel interpolation is carried out. We use ${\rm{C}}\left\{ \cdot \right\}$ to represent the Stolt interpolation operator, then Eq. (13) becomes

      $$ {S_{{\rm{Stolt}}}}\left( {{f_\tau }, {f_\eta }} \right) = {\rm{C}}\left\{ {{S_{{\rm{RFM}}}}\left( {{f_\tau }, {f_\eta }} \right)} \right\} $$

      The last step is 2-D inverse Fourier transform to get the final image $g\left( {\tau, \eta } \right)$:

      $$ g\left( {\tau, \eta } \right) = {\cal F}_\eta ^{ - 1}\left\{ {{\cal F}_\tau ^{ - 1}\left\{ {{S_{{\rm{Stolt}}}}\left( {{f_\tau }, {f_\eta }} \right)} \right\}} \right\} $$

      We can sample the continuous-time analog signals, ${s_0}\left( {\tau, \eta } \right)$, $g\left( {\tau, \eta } \right)$, into two-dimensional arrays, S $\in {{\mathbb{C}}^{{{N}_{\tau }}\times {{N}_{\eta }}}}$, G $\in {{\mathbb{C}}^{{{N}_{\tau }}\times {{N}_{\eta }}}}$, then the omega-K algorithm can be written as

      $$ \mathit{\boldsymbol{G=}}\left\{ \mathit{\boldsymbol{F}}_{\tau }^{\rm{H}}\left\langle C\left[\left( {{\mathit{\boldsymbol{F}}}_{\tau }}S \right){{\mathit{\boldsymbol{F}}}_{\eta }}{{\mathit{\boldsymbol{H}}}_{\rm{ref}}} \right] \right\rangle \right\}\mathit{\boldsymbol{F}}_{\eta }^{\rm{H}} $$

      where Fτ and Fη represent Fourier transform matrices in range and azimuth, respectively; ${\left\{ \cdot \right\}^{\rm{H}}}$ denotes the conjugate transposition.

    • At the beginning of this section, we have introduced that if there exists a substitution for A, then the efficiency of the algorithm can be improved. Inspired by this idea, we derive the inverse of omega-K algorithm. Since omega-K algorithm is used to get SAR image from echoes, then it is intuitive that we can get echoes from SAR image using the inverse of omega-K algorithm. From Eq. (3), we are able to say that the inverse of omega-K algorithm has the same effect as A, thus can substitute A in Eq. (7).

      We can easily conclude that omega-K is a linear algorithm that only contains Fourier transform, complex multiplication, and Stolt interpolation. The inverse of these three operations are also linear and the general procedures are quite obvious:

      $$ \mathit{\boldsymbol{S=F}}_{\tau }^{\rm{H}}\left\{ \mathit{\boldsymbol{F}}_{\eta }^{\rm{H}}\left[<D\left[{{\mathit{\boldsymbol{F}}}_{\tau }}\left( G{{\mathit{\boldsymbol{F}}}_{\eta }} \right) \right]>\mathit{\boldsymbol{H}}_{\rm{ref}}^{\rm{H}} \right] \right\} $$

      where D represents the inverse of Stolt interpolation. It takes a similar form as Eq. (15):

      $$ \sqrt{{{\left( {{f}_{0}}+f_{\tau }^{'} \right)}^{2}}+\frac{{{c}^{\text{2}}}{{f}_{\eta }}{{}^{\text{2}}}}{\text{4}V_{{{\text{r}}_{\text{ref}}}}^{\text{2}}}}={{f}_{0}}+{{f}_{\tau }} $$

      We substitute the square root with $ {f_0} + {f_\tau }$. The inverse of truncated sinc-kernel interpolation is hard to achieve directly. In our method, we also use a truncated sinc-kernel interpolation to approximate its inverse.

      Note that matrices in Eq. (19) could be very large, but we can perform Fourier transform and inverse Fourier transform directly instead of con-ducting matrices multiplications. It can be a signi-ficant save of time consumption. Meanwhile, we do not have to restore the large matrices either, which is also a save of memory cost.

      The inverse of omega-K algorithm can be concluded accordingly: Fourier transform in range and azimuth first; then the inverse of Stolt interpolation; followed by is multiplication with the conjugate transposition of the reference function; at last IFFT in range and azimuth.

      We use T to represent the procedure of inverse omega-K algorithm above, then Eq. (19) becomes

      $$ \mathit{\boldsymbol{S=}}T\left( \mathit{\boldsymbol{G}} \right) $$

      Our goal is to get G from S using T(·)·.

    • We can conclude from Eq. (21) that it is possible to get echoes from the scene directly and precisely if no noise exists. If there is noise in the signal, we can draw the following objective function for sparse scene with the help of CS theory:

      $$ \mathop {\min }\limits_\mathit{\boldsymbol{G}} {\mkern 1mu} \left\| {\mathit{\boldsymbol{S}} - {\rm{T}}\left( \mathit{\boldsymbol{G}} \right)} \right\|_2^2 + \lambda {\left\| \mathit{\boldsymbol{G}} \right\|_1} $$

      We know that there is a lot of redundancy in SAR signals and that CS theory has proved to be capable of reconstructing a sparse signal with much fewer measurements than Nyquist theory needs. So we sample the signals to reduce redundancy:

      $$ {{\mathit{\boldsymbol{S}}}_{\rm{S}}}={{\mathit{\Theta }}_{\tau }}\mathit{\boldsymbol{S}}{{\mathit{\Theta }}_{\eta }} $$

      where $ {{\mathit{\Theta }}_{\tau }}\in {{\mathbb{C}}^{N_{\tau }^{'}\times {{N}_{\tau }}}}, \left( N_{\tau }^{'} < {{N}_{\tau }} \right)$ and ${{\mathit{\Theta }}_{\eta }}\in {{\mathbb{C}}^{{{N}_{\eta }}\times N_{\eta }^{'}}}, \left( N_{\eta }^{'}\times {{N}_{\eta }} \right)$ represent random sampling operators in range and azimuth, respectively. Then the objective function becomes:

      $$ \mathop {\min }\limits_\mathit{\boldsymbol{G}} {\mkern 1mu} \left\| {{\mathit{\boldsymbol{S}}_{\rm{S}}} - {\mathit{\Theta }_\tau }T\left( \mathit{\boldsymbol{G}} \right){\mathit{\Theta }_\eta }} \right\|_2^2 + \lambda {\left\| \mathit{\boldsymbol{G}} \right\|_1} $$

      Since T is linear, Eq. (24) is a convex problem that can be solved with many convex optimization algorithms. In this study, we choose ITA to get G. In Section 2, we have introduced that Eq. (7) can be a solution to Eq. (6). Then in Section 3, we derived the inverse of omega-K, T, to substitute A. In order to solve the objective function Eq. (24) effectively without multiplications of large matrices, we still need to find a simpler operator to substitute AH. Since AH can be viewed as the inverse of A, it is clear that omega-K algorithm can be a substitution for AH. We use M to represent omega-K algorithm, then from Eq. (7) we can get the following solution to Eq. (24):

      $$ {{\mathit{\boldsymbol{G}}}^{\left( i+1 \right)}}={{E}_{1, \lambda \mu }}\left( {{\mathit{\boldsymbol{G}}}^{\left( i \right)}}+\mu M\left( {\mathit{\Theta }_\tau }^{\rm{T}}\left( {{\mathit{\boldsymbol{S}}}_{s}}-{\mathit{\Theta }_\tau }T\left( {{\mathit{\boldsymbol{G}}}^{\left( i \right)}} \right){{\mathit{\Theta }}_{\eta }} \right){{\mathit{\Theta }}_{\eta }}^{\rm{T}} \right) \right) $$

      M has a similar form as T, so it can save time and memory cost too. From Eq. (25), we can give the following algorithm procedure:

      Algorithm 1: Iterative thresholding algorithm for proposed CS-SAR imaging
      Input: SAR raw echoes SS,
          omega-K algorithm M and inverse omega-K algorithm T,
          sampling operator Θτ and Θη
      Initial: G(0), λ, μ, and max iteration Imax
          1: for i = 1 to Imax do
          2: residue: $ {{\mathit{\boldsymbol{R}}}^{\left( i-1 \right)}}={{\mathit{\boldsymbol{S}}}_{S}}-{\mathit{\Theta }_\tau }T\left( {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}} \right){{\mathit{\Theta }}_{\eta }}$
          3: omega-K on residue: $\Delta {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}}=M\left( \mathit{\Theta } _{\tau }^{\rm{T}}{{\mathit{\boldsymbol{R}}}^{\left( i-1 \right)}}\mathit{\Theta } _{\eta }^{\rm{T}} \right)$
          4: Thresholding: ${{\mathit{\boldsymbol{G}}}^{\left( i \right)}}={{E}_{1, \lambda \mu }}\left( {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}}+\mu \Delta {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}} \right)$
          5: end for
      Output: The recovery image G*=G(i)

      Let's generally compare the computational complexity of CS-OKA and other CS-SAR methods which also use ITA to solve the problem. First, we denote $N = {N_\tau } \times {N_\eta }$. In each iteration, the main computational complexity of CS-OKA results from omega-K, inverse omega-K, which have commonly the computational complexity of $ {\mathcal O}\left( {N{\rm{log}}N} \right)$, and the thresholding operator with complexity of ${\mathcal O}\left( N \right)$. The total complexity is ${\mathcal O}\left( {{I_{\max }}N{\rm{log}}N} \right)$. For other CS-SAR methods that use Eq. (7) to solve the problem, multiplication of two large matrices brings the complexity of ${\mathcal O}\left( {NsU} \right)$ in each iteration, where represents the sampling rate, and U is the product of samples of sent pulses and samples of azimuth time. Then the total complexity of other CS-SAR methods is ${\mathcal O}\left( {{I_{\max }}NsU} \right) $. In practice, U is usually very large to acquire high Signal-to-Noise Ratio (SNR). As a result, ${\mathcal O}\left( {{I_{\max }}NsU} \right)$ is much higher than ${\mathcal O}\left( {{I_{\max }}N{\rm{log}}N} \right) $. CS-OKA obviously has lower computational complexity.

      As for memory cost, since other CS-SAR methods that use ITA have to restore the big matrices A and AH, the memory cost would be ${\mathcal O}\left( {{N^2}} \right)$. For CS-OKA, we only need to restore the reference function Href. That leads to the memory cost of ${\mathcal O}\left( N \right)$. Apparently, CS-OKA has much lower memory cost than other CS-SAR methods.

      CS-OKA uses the redundancy of SAR signal to reduce the data rate compared with traditional MF-based methods. Furthermore, sparse optimization method is adopted to explore the sparse property of the signal, which could significantly suppress the sideslobes and improve the resolution. Also, it applies the inverse of omega-K algorithm to substitute the large matrix operations in other CS-SAR algorithms, which significantly contributes to the improvement of time complexity.

    • In this section, we will carry out several simulations to demonstrate the effectiveness of our method. The simulations were taken under different sampling rates and SNRs. In our simulations, signal in SNR refers to the echo. All simulations are performed under side-looking SAR. Tab. 1 lists the parameters of our simulations.

      Parameter Value
      Pulse duration (μs) 1.33
      Bandwidth in range (MHz) 150
      Carrier frequency (MHz) 600
      Sampling rate (MHz) 225
      Slant range of scene center (m) 1200
      Length of synthetic aperture (m) 300
      Pulse repeat frequency (Hz) 150 000 000
      Radar velocity in azimuth (m/s) 150

      Table 1.  Parameters used in the simulation

      In our simulations, we put 4 targets in the scene, all with the same unit amplitude. We first recover the scene using omega-K algorithm and CS-OKA with full sampled data. As shown in Fig. 1, both algorithms can recover the scene under different SNRs from 20 dB to 5 dB. Fig. 2 shows the contours of magnitude of the lower left target in Fig. 1. It is clear that results of CS-OKA has much lower sidelobes than omega-K algorithm. Then we undersample the signal (echo) and try to recover the scene with CS-OKA. In Fig. 3, we randomly sample 50% data in range and azimuth, respectively (i.e., for ${\mathit{\Theta }_{\tau }}\in {{\mathbb{C}}^{N_{\tau }^{'}\times {{N}_{\tau }}}}, N_{\tau }^{'}=0.5{{N}_{\tau }};\text{for}\ {\mathit{\Theta }_{\eta }}\in {{\mathbb{C}}^{{{N}_{\eta }}\times N_{\eta }^{'}}}, N_{\eta }^{'}=0.5{{N}_{\eta }} $). The results proved the validity of our method. We further undersample the signal to 10% data in both range and azimuth. Results are shown in Fig. 4. The scene can still be recovered perfectly even with 1% sampled data. In Fig. 5, we put two targets that are very close to each other. omega-K algorithm failed to separate the two targets because of high sidelobes [Figs. 5(a) and 5(b)], while CS-OKA can clearly distinguish the two targets [Figs. 5(c) and 5(d)].

      Figure 1.  Reconstruction results with full sampled data but different SNRs. The left column is recovered by omega-K algorithm, and the right column is recovered by CS-OKA. From top to bottom, SNRs are 20, 10, and 5 dB, respectively

      Figure 2.  Contours of magnitude of lower left target in Fig. 1. The left column is recovered by omega-K algorithm, and the right column is recovered by CS-OKA. From top to bottom, SNRs are 20, 10, and 5 dB, respectively

      Figure 3.  Reconstruction results with 50% samples in both range and azimuth. From left to right, SNRs are 20, 10, and 5 dB, respectively

      Figure 4.  Reconstruction results with 10% samples in both range and azimuth. From left to right, SNRs are 20, 10, and 5 dB, respectively

      Figure 5.  Reconstruction results of two close targets

      It is obvious that CS-OKA works very well even when the sample rate is low. Besides, CS-OKA maintains a high quality of the recovered images under different SNRs. It proves that our method is robust against noise because of the intrinsic property of Eq. (6). Also, CS-OKA has a higher resolution compared to omega-K algorithm.

    • In this study, we have proposed a new method for SAR imaging called CS-OKA. Different from previous CS-SAR methods, we combine the omega-K algorithm and CS to recover the scene.

      We first derived the inverse of omega-K algorithm to get echoes from the scene. It can be a substitution for a large matrix in the original ITA. The inverse of omega-K algorithm can be realized efficiently with the help of several linear operations, which is also a save for memory cost. Then we proved that omega-K algorithm can substitute the other large matrix for convenience. It has a similar form as its inverse. These two substitutions together have significantly improved the computational complexity since the multiplications of large matrices are avoided. We have analyzed the detailed computational complexity and proved the efficiency of CS-OKA.

      Simulations have shown that our method per-formed well under low sampling rate and is robust under different SNRs when the scene is sparse.

参考文献 (19)

目录

    /

    返回文章
    返回