合成孔径雷达稀疏成像方法
【专利摘要】一种SAR稀疏成像方法，包括：构造全孔径回波模拟算子M，并得到模拟回波数据；基于所述模拟回波数据来构建雷达观测方程；根据构建的雷达观测方程，建立基于SAR回波模拟算子的Lq正则化成像模型；以及采用阈值迭代方法求解建立的基于SAR回波模拟算子的Lq正则化成像模型，以便重建观测场景散射强度X*。根据本发明实施例的SAR成像方法，可以在低于奈奎斯特率采样下实现成像，并可抑制旁瓣，从而获得更清晰的SAR图像。
【专利说明】
合成孔径雷达稀疏成像方法
技术领域
[000?] 本发明实施例涉及合成孔径雷达(Synthetic Aperture Radar，SAR)成像，更具体 地，涉及一种SAR稀疏成像方法。
【背景技术】
[0002] SAR作为一种主动式微波成像系统，具有全天时、全天候和高分辨率成像等特点。 传统的全孔径滑动聚束成像方法通过解斜降低方位向带宽，然后对全孔径数据利用传统条 带聚焦算法进行成像，最后再通过方位向基带变标，从而解决方位向图像混叠问题，最终得 到聚焦图像。稀疏微波成像因其可以降低方位采样率，并得到良好重构图像而在雷达成像 中得以应用。但是，传统的稀疏微波成像框架都是基于雷达二维精确观测得到的雷达观测 矩阵，求解该模型的计算代价过于庞大，因此难以用于大场景的成像。
[0003] 对已有技术的了解可参考以下文献及其中的相关引文。
[0004] [ I ]Donoho D L Compressed sensing[J]. Information Theory ?IEEE Transactions on，2006,52(4):1289-1306.
[0005] [2]Zhang B C?Hong ff?ffu Y R.Sparse microwave imaging:Principles and applications[J].Science China Information Sciences?2012?55(8)：1722-1754.
[0006] [3]Fang J，Xu Z?Zhang B，et al.Fast compressed sensing SAR imaging based on approximated observation[J].Selected Topics in Applied Earth Observations and Remote Sensing，IEEE Journal of，2014,7(1):352-363.
[0007] [4]Khwaja A S，Ferro_Famil L，Pottier E.SAR Raw Data Simulation-Using High Precision Focusing Methods[J].EUSAR 2006,2006.
[0008] [5]Patel V M，Easley G R，Healy Jr D M，et al·Compressed synthetic aperture radar[J].Selected Topics in Signal Processing，IEEE Journal of，2010,4 (2):244-254.
[0009] [6]Prats P? Scheiber R，Mittermayer J，et al. Processing of sliding spotlight and TOPS SAR data using baseband azimuth scaling[J].Geoscience and Remote Sensing，IEEE Transactions on，2010,48(2):770-780.
[0010] [7]Sun G?Xing M?ffang Y?et al.Sliding spotlight and TOPS SAR data processing without subaperture[J].Geoscience and Remote Sensing Letters，IEEE， 2011,8(6):1036-1040.
【发明内容】
[0011] 根据本发明实施例的一个方面，提供了一种SAR稀疏成像方法，包括：
[0012] 步骤Sl，构造全孔径回波模拟算子M，并得到模拟回波数据：
[0013] f^M(X)
[0014] 其中X表示目标场景的散射强度矩阵，f是模拟回波数据，M是全孔径回波模拟算 子，其中由脉冲压缩算法的逆过程来构造 Μ。
[0015] 步骤S2,基于所述模拟回波数据来构建雷达观测方程：
[0016]
[0017] 其中Ys表示经二维采样后的回波数据，N是噪声，Θη，Θτ分别表示方位采样矩阵和 距离随机降采样矩阵，M是构造的全孔径回波模拟算子。
[0018] 步骤S3,根据构建的雷达观测方程，建立基于SAR回波模拟算子的Lq正则化成像模 型：
[0019]
[0020] 其中，X#是重建的目标场景散射强度，arg min是最小化计算式，I |X| |q为X的q(q = 1)范数。
[0021] 步骤S4,采用阈值迭代方法求解建立的基于SAR回波模拟算子的Lq正则化成像模 型，以便重建目标场景的散射强度浐。
[0022] 根据本发明实施例的技术方案，利用SAR回波特性及观测场景的稀疏性，建立基于 滑动聚束SAR回波模拟算子的稀疏正则化模型。此外，利用融合回波模拟算子的阈值迭代方 法实现对观测区域目标场景雷达成像。相比于已有的基于二维观测模型的稀疏SAR成像算 法，根据本发明实施例的成像方法能够提高运行效率并降低计算成本。此外，相比于已有的 匹配滤波成像方法，根据本发明实施例的方法可以在低于奈奎斯特率采样下实现成像，并 可抑制旁瓣，从而获得更清晰的SAR图像。
【附图说明】
[0023] 仅作为示例参照附图对本发明实施例进行详细描述，在附图中：
[0024] 图1示出了根据本发明实施例的基于全孔径回波模拟算子的滑动聚束稀疏SAR成 像方法的流程图；
[0025] 图2示出了根据本发明实施例的构造全孔径回波模拟算子的流程图；
[0026] 图3示出了根据本发明实施例的利用阈值迭代求解Lq正则化成像模型的流程图； 以及
[0027] 图4示出了采用传统雷达成像与采用根据本发明实施例的基于模拟回波算子的稀 疏SAR成像方法得到的二维仿真成像结果的比较。
【具体实施方式】
[0028] 下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是， 所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。
[0029] 图1示出了根据本发明实施例的SAR稀疏成像方法的流程图。如图1所示，在步骤 Sl，构建回波模拟算子M，并得到模拟回波数据：
[0030] Y =M(X) ( I )
[0031] 其中X代表目标场景的回波散射强度矩阵，:f为模拟回波数据;M为回波模拟算子， 由脉冲压缩算法的逆过程构造。
[0032]优选地，根据本发明实施例，还可以包括步骤:获取目标场景的回波散射强度矩阵 X。
[0035]
[0033] 优选地，如图2所示，步骤Sl具体可以包括，[0034]对目标散射强度矩阵X进行逆去斜相位补偿相乘得到S1:
[0036]
[0037] 其中，H1是逆方位向基带变标的补偿相位，V是雷达平台速度，λ是发射信号波长， r r。t为滑动聚束旋转中心最短斜距，r Q为场景中心最短斜距。
[0038] 然后，对S1进行方位向傅里叶变换得到S2a:
[0039] S2a = FFTn(Si)
[0040]其中:FFTn表示方位向傅里叶变换。
[0041 ]然后，对S2a进行逆去斜操作，得到S2b:
[0042]
[0043]
[0044] 其中，H2是逆方位向基带变标相位。
[0045] 然后，对S2b进行方位向相位补偿得到&。：
[0046]
[0047]
[0048]
[0049]
[0050]其中：Φ4为相位补偿因子，c为光速，Kr为发射信号线性调频率1为雷达脉冲和SRC 滤波器的综合调频率，Cs(f)为变标因子。
[0051] 然后，对&。进行方位向解压缩得到S2d:
[0052]
[0053]
[0054] 其中，Φ3为方位解压缩因子，r为最短斜距。
[0055] 然后，对S2d进行距离向傅里叶变换得到S2e:
[0056] S2e = FFTx (S2d)
[0057] FFlMf表距离向快速傅里叶变换。
[0058] 对S2e进行距离徙动、距离向解压缩及二次距离徙动得到S2f:
[0059]
[0060]
[0061]其中，Φ2为距离徙动、距离向解压缩及二次距离徙动因子，为距离向频率。
[0062]然后，对S2f进行逆距离向傅里叶变换得到S2g:
[0063] S2g=IFFTx(S2f)
[0064] IFFIV表示距离向逆傅里叶变换。
[0065]然后，对S2g进行逆线性调频变标得到S2h:
[0068] Φι = θχρ( j3iKs(f ,r〇)Cs(f) [τ-τ0(?)])[0069] 其中，Φ:为逆线性调频变标因子，T〇(f)为距离向压缩位置。[0070] 然后，对S2h进行逆方位向傅里叶变换得到模拟S2:[0071] S2 = IFFTn(S2h)[0072] IFFTn代表逆方位向快速傅里叶变换[0073] 然后，对S2进行逆去斜相位相乘操作，得到S3a:[0074]
[0066]
[0067]
[0075]
[0076] 然后，对S3a进行方位向逆傅里叶变换，得到S3b:
[0077] Ssb = IFFTn(S3a)
[0078] 然后，对S3b进行逆去斜操作，得到模拟回波数据f :
[0079]
[0080] 如图1所示，在步骤S2,根据得到的模拟回波数据构建雷达观测模型： _
⑴
[0082]其中Ys表示经二维采样后的回波数据，N为噪声，Θη，Θτ分别代表方位采样矩阵和 距离随机降采样矩阵。
[0083]然后，在步骤S3,根据观测模型（1)建立基于滑动聚束SAR回波模拟算子的Lq正则 化成像模型：
[0084]
(2)
[0085] 然后，在步骤S4,采用阈值迭代方法求解基于滑动聚束SAR回波模拟算子的Lq正则 化成像模型(2)，从而重建目标场景的散射强度X'
[0086] 如图3所示，步骤S4可以包括：
[0087] 初始化目标场景的散射强度XQ，并设定目标场景的稀疏度预估值K和迭代终止准 贝 1J，令 n = l;
[0088] 然后，按照下式更新梯度下降序列Bn:
[0089]
[0090]其中μ表示梯度下降的步长，
[0091]
[0092] %，//1冲,3沖#1，4分别为!11，!12,〇4,〇3,〇2,〇1，8 1)的复共辄。
[0093] 然后，更新正则化参数：
[0094]
[0095] 其中|Bn|k+1表示将序列Bn的模值按降序排列后的第k+Ι个元素；
[0096]然后，更新目标场景的散射强度Xn+1
[0097] Xn+i = Hq (Bn)
[0098] 其中Hq( ·)为阈值算子：
[0099] Hq(X) = (hq(Xl)，…，hq(Xn) )τ
[0100] 其中当q = l时
[0101]
[0102] 然后，若I |χη+1-χη| If/I Ιχη| |「<￡迭代终止，输出是目标场景的回
[0103] 波强度;否则，令η = η+1，返回更新梯度下降序列仏的步骤。
[0104]其中，ε是终止迭代的门限。
[0105] 因此，结合图3,模拟回波数据可以表示为：
[0106]
[0107] 图4示出了采用传统雷达成像与采用根据本发明实施例的SAR稀疏成像方法得到 的二维仿真成像结果的比较。其中，主要雷达参数包括:景中心距离r〇 = 671km，旋转中心距 离:Trot = 1061.5km飞行器速度V=7391m/s，信号带宽Br = I50MHz，脉冲持续时间Tr = IOys，载 频fo = 9.65GHz，脉冲发射频率 PRF = 3798Hz。
[0108] 图4(a)示出了已有的全孔径匹配滤波方法的仿真结果。图4(b)、图4(c)、图4(d)分
【主权项】
1. 一种SAR稀疏成像方法，包括： 步骤S1，构造全孔径回波模拟算子M，并得到模拟回波数据； 步骤S2,基于所述模拟回波数据来构建雷达观测方程； 步骤S3,根据构建的雷达观测方程，建立基于SAR回波模拟算子的Lq正则化成像模型；W 及 步骤S4,采用阔值迭代方法求解建立的基于SAR回波模拟算子的Lq正则化成像模型，W 便重建目标场景的回波散射强度X*。2. 根据权利要求1所述的方法，其中，所述步骤S1包括利用W下公式来得到模拟回波数 据：其中X表示目标场景的回波散射强度矩阵，:f是模拟回波数据，Μ是全孔径回波模拟算 子，其中由脉冲压缩算法的逆过程来构造 Μ。3. 根据权利要求1所述的方法，其中，所述步骤S2包括利用W下公式来构建雷达观测方 程：其中Ys表示经二维采样后的回波数据，Ν是噪声，Θη，Θτ分别表示方位采样矩阵和距离 随机降采样矩阵，Μ是构造的全孔径回波模拟算子。4. 根据权利要求1所述的方法，其中，所述步骤S3包括利用W下公式来建立Lq正则化成 像模型：其中，X*是重建的回波散射强度矩阵，Ao为一常数，arg min是最小化计算式，IIX II q为X 的q(q = l)范数。5. 根据权利要求2所述的方法，其中，所述步骤S1还包括： 对目标场景的回波散射强度数据X进行逆去斜相位补偿相乘得到Si; 对Si进行方位向傅里叶变换得到S2a; 对S2a进行逆去斜操作，得到S2b ; 对S2b进行方位向相位补偿得到S2。； 对S2。进行方位向解压缩得到S2d ; 对S2d进行距离向傅里叶变换得到S2e; 对S2e进行距离徙动、距离向解压缩及二次距离徙动得到S2f ; 对S2f进行逆距离向傅里叶变换得到S2g; 对S2g进行逆线性调频变标得到S2h，W及对S2h进行逆方位向傅里叶变换得到模拟S2; 对S2进行逆去斜相位相乘操作，得到S3a ; 对S3a进行方位向逆傅里叶变换，得到S3b; 对S3b进行逆去斜操作，得到模拟回波数据f。 6 .根据权利要求1所述的方法，其中，所述步骤S4还包括： 初始化目标场景的回波散射强度Xo，设定目标场景的稀疏度预估值Κ和迭代终止准则； 更新梯度下降序列Bn; 更新正则化参数λ; 更新目标场景的散射强度Χη+1; 若|帖+片。|^/|1乂。|^<6迭代终止，输出是目标场景的回波强度;否则，令11 = 11+1，并 返回所述更新梯度下降序列Bn的步骤。 7 .根据权利要求6所述的方法，其中，所述更新梯度下降序列Bn包括： 其中μ表示梯度下降的步长，幻W分别为出，也，巫4，巫3，巫2，巫1，Sp的复共辆。8. 根据权利要求6所述的方法，其中，所述更新正则化参数λ包括 义=打巧,L;其中，q=i 其中|Bn|k+l表示将序列Bn的模值按降序排列后的第k+1个元素。9. 根据权利要求6所述的方法，其中，所述更新目标场景的散射强度Xn+1包括 Xn+l = Hq 化η) 其中Hq( ·)为阔值算子： Hq(X)=化q(Xl)，···，hq(Xn))T 其中当q=l时10. 根据权利要求1所述的方法，还包括:在构造全孔径回波模拟算子Μ之前，获取目标 场景的散射强度数据X。
【文档编号】G01S13/90GK105842699SQ201610384461
【公开日】2016年8月10日
【申请日】2016年6月2日
【发明人】张冰尘, 魏中浩, 毕辉, 吴戎, 吴一戎
【申请人】中国科学院电子学研究所