CN102788977B - 基于l1/2正则化的合成孔径雷达成像方法 - Google Patents

Abstract

本发明公开了一种基于L_1/2正则化的合成孔径雷达成像方法，涉及雷达二维成像技术，包括步骤：a、根据合成孔径雷达观测模型，建立基于L_1/2正则化的合成孔径雷达成像模型；b、采用迭代Half阈值算法实现对观测场景后向散射系数的重建。本发明方法相比于传统的合成孔径雷达成像方法，可降低正确重建目标场景所需的采样量，实现合成孔径雷达数据有效成像。

Description

基于L1/2正则化的合成孔径雷达成像方法

技术领域

本发明涉及雷达二维成像技术领域，是一种基于L_1/2正则化的合成孔径雷达成像方法，基于现代稀疏信号处理方法的成像处理技术。 

背景技术

合成孔径雷达(Synthetic Aperture Radar，SAR)作为一种主动式微波成像系统，具有全天时、全天候和高分辨率成像特点。它已经被广泛应用于军事及国民经济的许多领域，如军事侦察、环境监测、土地资源管理等方面。随着SAR技术的发展，要求雷达系统的分辨率和测绘带宽不断提高，大数据量的瓶颈效应也越发明显。 

传统SAR成像的采样数据是按照Nyquist率进行采样。Nyquist采样定理指出，对于有限带宽的实信号，为了保证无失真的恢复信号所需要的采样率应不低于两倍的信号带宽。基于Nyquist率的采样不可避免地产生大量数据，从而难以存储、处理和传输。然而，许多SAR场景在某些特定的基底下具有稀疏性。对于稀疏场景的成像技术已经有了许多的研究。Cetin提出了一种基于非二次正则化SAR成像方法用于稀疏场景的重建。该方法利用整个相位历史进行重建，且正则化的模型参数需要人工事先选取。最近兴起的压缩感知技术获得了大量的关注。它主要通过随机测量，采用非线性重建算法得到线性反问题的稀疏解。不同于传统的采样方式，压缩感知利用信号的稀疏性，将压缩与采样合并，从而以低于传统Nyquist率进行采样，实现信号的恢复。近年来，压缩感知被广泛应用于SAR成像。近期中国科学院电子学研究所提出稀疏微波成像新体制。目前压缩感知成像算法主要是基于L₁正则化实现对观测场景后向散射系数的恢复。和传统的SAR成像方法相比，基于L₁正则化的压缩感知的成像方法能大大降低SAR系统所需的数据量。然而基于L₁正则化的成像方法所需的采样量仍然相对较高。 

近期，Xu建立了L_1/2正则化框架，给出了L_1/2正则化模型解的解析表达式，提出了一种求解L_1/2正则化模型的快速算法。通过相变图研究，揭示了L_1/2正则化在L_q(0＜q＜1)正则化中具有代表性。相比于L₁正则化，L_1/2正则化具有更强的稀疏能力。 

对背景技术的了解可以参考下面三篇文章以及它们的引用文章。 

[1]M.Cetin and W.Karl，Feature-enhanced synthetic aperture radar image formation based on non-quadratic regularization，IEEE Trans.on Image Processing，vol.10，no.4，pp.623-631，2001. 

[2]D.L.Donoho，“Compressed sensing，”IEEE Trans.Inform.Theory，2006，52，(4)，pp.1289-1306. 

[3]Z.B.Xu，“Data modeling：Visual psychology approach and L1/2regularization theory，”Proceedings of International Congress of Mathematicians，Switzerland：European Mathematical Society Publishing House，2010. 

发明内容

本发明提供了一种基于L_1/2正则化的合成孔径雷达成像方法，根据雷达观测模型，建立基于L_1/2正则化的合成孔径雷达成像模型；采用迭代Half阈值算法实现对观测区域目标场景后向散射系数的重建。本发明方法相比于传统的合成孔径雷达成像方法，可降低正确重建目标场景所需的采样量，实现合成孔径雷达数据有效成像。 

为实现上述目的，本发明的技术解决方案是： 

一种基于L_1/2正则化的合成孔径雷达成像方法，其包括步骤： 

A)根据合成孔径雷达观测模型，建立基于L_1/2正则化的SAR成像模型； 

B)采用迭代Half阈值算法重建观测场景后向散射系数。 

所述的合成孔径雷达成像方法，其所述A)步，包括： 

a、根据雷达观测模型建立基于L_1/2正则化的合成孔径雷达成像模型： 

$x^{*}=\arg \min ({\left|\right|y-\mathrm{\Φx}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|x\left|\right|}_{1/2}^{1/2})---\left(1\right)$

其中，x^*是目标最优值，y是雷达观测回波数据，x是目标场景后向散射系数，Ф是雷达系统观测矩阵，λ是正则化参数，arg min是最小化计算式； 

b、建立的基于L_1/2正则化的合成孔径雷达成像模型的解，具有对角非线性阈值形式，换言之，其解可表示为x^*＝H_λμ，1/2(B_μ(x^*))，其中H_λμ，1/2(·)为Half阈值迭代算子，B_μ(x^*)＝x^*+μФ^H(y-Фx^*)，Ф^H为矩阵Ф的共轭转置矩阵，μ∈(0，1)为一常数。 

具体地，对任意 n为向量x的维数，x_i为向量x第i个分量，i＝1，2，…，n， 为复数域，H_λμ，1/2(x)＝(h_λμ，1/2(x₁)，…，h_λμ，1/2(x_n))^T，h_λμ，1/2(·)为Half阈值函数，其中： 

和 

为一带有参数λμ的函数，具体定义见上式，cos(·)为余弦函数，arccos(·)为反余弦函数。 

所述的合成孔径雷达成像方法，其所述B)步，根据所得到的L_1/2正则化的雷达成像解的阈值表达形式，设计出一种求解L_1/2正则化模型的快速算法，称为迭代Half阈值算法(Iterative Half Thresholding Algorithm)，主要步骤如下： 

步骤1：初始化目标场景的后向散射系数x₁，给定梯度下降步长μ、目标场景稀疏度预估值K和迭代终止准则，令迭代步数t:＝1； 

步骤2：更新梯度下降序列B_t： 

B_t＝x_t+μФ^H(y-Фx_t)(2) 

其中，Ф^H表示Ф的共轭转置矩阵； 

步骤3：更新正则化参数λ_t； 

步骤4：更新目标场景的后向散射系数： 

$x_{t+1}=H_{{\mathrm{\λ}}_t\mathrm{\μ},1/2}\left(B_t\right)---\left(3\right)$

其中， 为Half阈值算子； 

步骤5：判断是否满足算法迭代终止准则？满足，迭代终止，算法输出为目标场景重建后向散射系数；不满足，则更新迭代步数，令t:＝t+1，转步骤2。 

所述的合成孔径雷达成像方法，其根据L_1/2正则化理论，给出基于L_1/2正则化的合成孔径雷达成像方法的正则化参数最优选取策略： 

${\mathrm{\λ}}^{*}\∈[\frac{\sqrt{96}}{9\mathrm{\μ}}{\left({\left|B_{\mathrm{\μ}}\left(x^{*}\right)\right|}_{k+1}\right)}^{3/2},\frac{\sqrt{96}}{9\mathrm{\μ}}{\left({\left|B_{\mathrm{\μ}}\left(x^{*}\right)\right|}_k\right)}^{3/2}),$

其中 

B_μ(x^*)＝x^*+μФ^H(y-Фx^*)， 

x^*为目标场景，Ф^H记为Ф的共轭转置矩阵，|B_μ(x^*)|_k和|B_μ(x^*)|_k+1分别记为B_μ(x^*)对应的模值向量第k和k+1个最大的分量，k是目标场景的稀疏度，即非零元素的个数。 

所述的合成孔径雷达成像方法，其在算法迭代过程中，取迭代序列 来逼近最优的正则化参数λ^*，t为迭代步数。 

所述的合成孔径雷达成像方法，其为保证迭代Half阈值算法的收敛性，梯度下降步长μ应满足  表示矩阵Ф的2-范数平方；算法迭代初始值x₁取为零；迭代Half阈值算法对目标场景稀疏度预估值K具有稳健性，在实际中，K取为Ф^Hy相应40dB输出分量的数目，即所有幅度大于Ф^Hy中绝对值最大分量幅度的百分之一的分量数目。 

本发明方法相比于传统的合成孔径雷达成像方法，可降低正确重建目标场景所需的采样量，实现合成孔径雷达数据有效成像。 

附图说明

图1为采用迭代Half阈值算法实现基于L_1/2正则化合成孔径雷达成像 的处理流程图； 

图2a为传统成像结果； 

图2b为基于L_1/2正则化SAR成像方法成像结果； 

图3a为传统成像方法满采样成像结果； 

图3b为基于L_1/2正则化SAR成像方法10％采样成像结果； 

图3c为基于L₁正则化SAR成像方法10％采样成像结果； 

图3d为基于L_1/2正则化SAR成像方法5％采样成像结果。 

具体实施方式

本发明的一种基于L_1/2正则化合成孔径雷达成像方法，其雷达观测回波模型如下： 

y＝Φx+n    (4) 

其中y是雷达观测回波数据，x是目标场景后向散射系数，Ф是雷达系统观测矩阵，n是接收热噪声。 

根据雷达观测回波模型，基于L_1/2正则化，建立如下的SAR成像模型 

$x^{*}=\arg \min ({\left|\right|y-\mathrm{\Φx}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|x\left|\right|}_{1/2}^{1/2})---\left(5\right)$

其中，x^*是目标最优值，y是雷达观测回波数据，x是目标场景后向散射系数，Ф是雷达系统观测矩阵，λ是正则化参数，arg min是最小化计算式。 

根据L_1/2正则化理论，所建立的基于L_1/2正则化的成像模型的解具有对角非线性阈值形式，换言之，其解可表示为x^*＝H_λμ，1/2(B_μ(x^*))。其中H_λμ，1/2(·)为Half阈值迭代算子，B_μ(x^*)＝x^*+μΦ^H(y-Фx^*)，Ф^H为矩阵Ф的共轭转置矩阵，μ∈(0，1)为一常数。具体地，对任意 n为向量x的维数，x_i为向量x第i个分量，i＝1，2，…，n， 为复数域，H_λμ，1/2(x)＝(h_λμ，1/2(x₁)，…，h_λμ，1/2(x_n))^T，h_λμ，1/2(·)为Half阈值函数，其中 

和 

为一带有参数λμ的函数，具体定义见上式，cos(·)为余弦函数，arccos(·)为反余弦函数。 

根据所得到的L_1/2正则化模型解的阈值表达形式，可以设计出一种求解L_1/2正则化模型的快速算法，即迭代Half阈值算法(Iterative Half Thresholding Algorithm)。采用迭代Half阈值算法重建目标场景后向散射系数，算法步骤如下： 

步骤1：初始化目标场景的后向散射系数x₁，给定梯度下降步长μ，目标场景稀疏度预估值K和迭代终止准则，令迭代步数t:＝1； 

步骤2：更新梯度下降序列B_t

B_t＝x_t+μФ^H(y-Фx_t)(6) 

其中Ф^H表示Ф的共轭转置矩阵； 

步骤3：更新正则化参数λ_t

${\mathrm{\λ}}_t=\frac{\sqrt{96}{\left({\left|B_t\right|}_{K+1}\right)}^{3/2}}{9\mathrm{\μ}}---\left(7\right)$

其中|B_t|_K+1表示B_t对应的模值向量第K+1个最大的分量； 

步骤4：更新目标场景的后向散射系数 

$x_{t+1}=H_{{\mathrm{\λ}}_t\mathrm{\μ},1/2}\left(B_t\right)---\left(8\right)$

其中 为Half阈值算子； 

步骤5：判断是否满足算法迭代终止准则？满足，迭代终止，算法输出为目标场景的后向散射系数；否则，更新迭代步数，令t:＝t+1，转步骤2。 

下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。 

见图1，采用迭代Half阈值算法实现基于L_1/2正则化的合成孔径雷达 成像处理流程的具体实施步骤如下： 

步骤1：后向散射系数初始值x₁取为零，给定梯度下降步长  表示矩阵Ф的2-范数平方，目标场景稀疏度预估值K取为Ф^Hy相应40dB输出分量的数目，给定误差容限ε及迭代终止准则||x_t+1-x_t||₂＜ε(也可以选取其他的迭代终止准则，如||x_t+1-x_t||₂/||x_t||₂＜ε或根据需要设定最大迭代步数MaxIter等)，令初始迭代步数t:＝1； 

步骤2：更新梯度下降序列B_t

B_t＝x_t+μФ^H(y-Фx_t)(9) 

其中Ф^H表示Ф的共轭转置矩阵； 

步骤3：更新正则化参数λ_t

${\mathrm{\λ}}_t=\frac{\sqrt{96}{\left({\left|B_t\right|}_{K+1}\right)}^{3/2}}{9\mathrm{\μ}}---\left(10\right)$

其中|B_t|_K+1表示B_t对应的模值向量第K+1个最大的分量，K是目标场景的稀疏度预估值； 

步骤4：更新观测目标场景后向散射系数 

$x_{t+1}=H_{{\mathrm{\λ}}_t\mathrm{\μ},1/2}\left(B_t\right)---\left(11\right)$

其中 为Half阈值算子； 

步骤5：判断是否满足算法迭代终止准则？满足，迭代终止，算法输出为目标场景后向散射系数；否则，更新迭代步数，令t:＝t+1，转步骤2。 

见图2，传统雷达成像及基于L_1/2正则化SAR成像方法一维成像结果如下： 

图2a为传统成像结果，图2b为基于L_1/2正则化SAR成像方法成像结果。 

见图3，分别采用传统成像、基于L₁正则化SAR成像与基于L_1/2正则化SAR成像方法处理RADARSAT-1实验数据成像结果如下： 

图3a为传统成像方法满采样成像结果，图3b为基于L_1/2正则化SAR成像方法10％采样成像结果，图3c为基于L₁正则化SAR成像方法10％采样成像结果，图3d为基于L_1/2正则化SAR成像方法5％采样成像结果。 

Claims (5)

1.一种基于L_1/2正则化的合成孔径雷达成像方法，其特征在于：包括步骤：

A)根据合成孔径雷达观测模型，建立基于L_1/2正则化的SAR成像模型；

B)采用迭代Half阈值算法重建观测场景后向散射系数；

其中，所述B)步，根据所得到的基于L_1/2正则化的SAR成像模型解的阈值表达形式，设计出一种SAR成像快速算法，称为迭代Half阈值算法(Iterative Half Thresholding Algorithm)，主要步骤如下：

步骤1：初始化目标场景的后向散射系数x₁，给定梯度下降步长μ、目标场景稀疏度预估值K和迭代终止准则，令迭代步数t：＝1；

步骤2：更新梯度下降序列B_t

B_t＝x_t+μΦ^H(y-Φx_t)    (2)

其中，y是雷达观测回波数据，Φ是雷达系统观测矩阵，Φ^H表示Φ的共轭转置矩阵；

步骤3：更新正则化参数λ_t；

步骤4：更新目标场景的后向散射系数：

x_t+1＝H_λtμ，1/2(B_t)    (3)

其中，H_λtμ，1/2(·)为Half阈值算子；

步骤5：判断是否满足算法迭代终止准则？满足，迭代终止，算法输出为目标场景重建后向散射系数；不满足，则更新迭代步数t，令t：＝t+1，转步骤2。

2.根据权利要求1所述的合成孔径雷达成像方法，其特征在于：所述A)步，包括：

a、根据雷达观测模型建立基于L_1/2正则化的合成孔径雷达成像模型：

$x^{*}=\arg \min ({\left|\right|y-\mathrm{\Φx}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|x\left|\right|}_{1/2}^{1/2})---\left(1\right)$

其中，x^*是目标最优值，y是雷达观测回波数据，x是目标场景后向散射系数，Φ是雷达系统观测矩阵，λ是正则化参数，arg min是最小化计算式；

b、建立的基于L_1/2正则化的合成孔径雷达成像模型的解，具有对角非线性阈值形式，换言之，其解可表示为x^*＝H_λμ，1/2(B_μ(x^*))，其中H_λμ，1/2(·)为Half阈值迭代算子，B_μ(x^*)＝x^*+μΦ^H(y-Φx^*)，Φ^H为矩阵Φ的共轭转置矩阵，μ∈(0，1)为一常数；

具体地，对任意n为向量x的维数，x_i为向量x第i个分量，i＝1，2，…，n，为复数域，H_λμ，1/2(x)＝(h_λμ，1/2(x₁)，…，h_λμ，1/2(x_n))^T，h_λμ，1/2(·)为Half阈值函数，其中：

和

为一带有参数λμ的函数，具体定义见上式，cos(·)为余弦函数，arccos(·)为反余弦函数。

3.根据权利要求1所述的合成孔径雷达成像方法，其特征在于：

根据L_1/2正则化理论，给出基于L_1/2正则化合成孔径雷达成像模型的正则化参数最优选取策略：

${\mathrm{\λ}}^{*}\∈(\frac{\sqrt{96}}{9\mathrm{\μ}}{\left({\left|B_{\mathrm{\μ}}\left(x^{*}\right)\right|}_{k+1}\right)}^{3/2},\frac{\sqrt{96}}{9\mathrm{\μ}}{\left({\left|B_{\mathrm{\μ}}\left(x^{*}\right)\right|}_k\right)}^{3/2}),$

其中

B_μ(x^*)＝x^*+_μΦ^H(y-Φx^*)，

λ^*为最优的正则化参数，x^*为目标场景，Φ^H记为Φ的共轭转置矩阵，|B_μ(x^*)|_k和|B_μ(x^*)|_k+1分别记为B_μ(x^*)对应的模值向量第k和k+1个最大的分量，k是目标场景的稀疏度，即非零元素的个数。

4.根据权利要求1或3所述的合成孔径雷达成像方法，其特征在于：在算法迭代过程中，取迭代序列来逼近最优的正则化参数λ^*，t为迭代步数。

5.根据权利要求1所述的合成孔径雷达成像方法，其特征在于：

为保证迭代Half阈值算法的收敛性，梯度下降步长μ应满足 表示矩阵Φ的2-范数平方；算法迭代初始值x₁取为零；迭代Half阈值算法对目标场景稀疏度预估值K具有稳健性，在实际中，K取为Φ^Hy相应40dB输出分量的数目，即所有幅度大于Φ^Hy中绝对值最大分量幅度的百分之一的分量数目。