一种基于柯西正则化的阵列三维SAR稀疏成像方法

一种基于柯西正则化的阵列三维sar稀疏成像方法
技术领域
1.本发明属于合成孔径雷达技术领域，它特别涉及到了合成孔径雷达(synthetic aperture radar,sar)稀疏成像技术领域。


背景技术：

2.合成孔径雷达(synthetic aperture radar,sar)是一种具有全天时、全天候的成像能力，即能够在白天或是黑夜、晴天或是雨雪天气随时成像，克服了光学和红外系统不能在晚上和复杂天气条件进行成像的缺点。传统的sar成像一般只具有二维成像分辨率，因此存在空间的一些重要信息(比如高度)丢失，为了获取完整的空间信息，实现三维成像是非常有必要的。目前常见的三维成像技术有圆周sar三维成像、层析sar三维成像、线阵sar三维成像。
3.阵列三维sar是一种将宽带雷达信号处理与天线阵列处理相结合，在切方位向添加阵列天线，平台沿方位向运动形成虚拟的面阵进而获得二维分辨率，距离向再通过脉冲压缩获得第三维的分辨率。相比于圆周sar三维成像，线阵sar三维成像不需要圆周运动的轨迹；相比于层析sar三维成像需要航过多次，线阵sar三维成像只需一次航过，所以线阵sar三维成像相对于层析sar和圆周sar三维成像有更强的灵活性。因此，在目标检测与识别、反恐安检、智能汽车雷达等领域具有重要的应用价值。
4.目前，经典的基于匹配滤波的成像方法，如后向投影，受到分辨率的限制，获得的三维图像中包含了较高的目标旁瓣与杂波，容易被错误判别为真实存在的目标，从而降低了图像的质量，对后续图像应用，如目标检测与识别等，造成不利影响。因此有必要提出能实现更高成像质量的成像方法，即能够获得更高成像分辨率、更低目标旁瓣与杂波的图像。其中，主要的改进方法便是稀疏成像方法。
5.稀疏成像方法通过引入关于目标场景的稀疏先验信息，构建了带稀疏正则项的数学优化问题，利用数学优化求解方法，获得成像结果。现有的阵列三维sar稀疏成像方法主要采用基于l1范数的稀疏正则化项，这会导致目标能量存在部分损失，尤其是对于能量较弱的目标，这会导致目标的视觉信息丢失，如目标的细节轮廓，这对于后续目标的识别、分类等应用带来了不利影响。因此，有必要提出能够更好的保留目标能量，完善目标细节轮廓的稀疏成像方法。为了满足这一需求，本发明提出了一种基于柯西正则化的阵列三维sar稀疏成像方法。


技术实现要素：

6.本发明提出了一种基于柯西正则化的阵列三维sar稀疏成像方法。首先，它利用阵列三维sar获取被测的目标场景回波数据矩阵；利用脉冲压缩和频率升采样技术完成信号的初步处理；利用三维bp算法得到三维初始图像；然后利用基于柯西正则化项的近端求解算子实现三维增强成像，最终获得目标场景三维稀疏成像结果。与基于匹配滤波的成像方法比，本发明方法能够提高分辨率、抑制杂波与噪声，提高成像质量，与基于l1范数正则化
项的稀疏成像方法相比，本发明方法能够在提高图像质量的同时，更好得保留目标能量，保留目标的细节轮廓。
7.为了方便描述本发明的内容，首先作以下术语定义：
8.定义1.合成孔径雷达
9.合成孔径雷达((synthetic aperture radar,sar)是一种高分辨率微波成像雷达，具有全天时和全天候工作的优点，已被广泛应用各个领域，如地形测绘、制导、环境遥感和资源勘探等。sar应用的重要前提和信号处理的主要目标是通过成像算法获取高分辨、高精度的微波图像。详见“皮亦鸣,杨建宇,付毓生,杨晓波.合成孔径雷达成像原理[m].电子科技大学出版社.2007”。
[0010]
定义2.步进频率(stepped-frequency,sf)信号
[0011]
sf信号是一种通过发射载频步进变化的子脉冲串来合成大的等效带宽、具备高分辨距离像能力的信号形式。它能够显著降低系统的瞬时带宽和对接收机硬件的需求，在微波成像与微波测量领域被广泛使用。详见“杨汝良.高分辨率微波成像[m].国防工业出版社.2013”。
[0012]
定义3.脉冲压缩
[0013]
脉冲压缩是一种现代雷达信号处理技术，简单来说就是雷达发射宽脉冲，然后再接收端“压缩”为窄脉冲，从而改善雷达的两种性能：作用距离和距离分辨率。详见“皮亦鸣,杨建宇,付毓生,杨晓波.合成孔径雷达成像原理[m].电子科技大学出版社.2007”。
[0014]
定义4.升采样
[0015]
升采样是一种在离散信号域提高信号采样率的方法，有时域升采样和频域升采样两种实现方式。
[0016]
定义5.快速傅里叶变换对算法
[0017]
快速傅里叶变换对算法是计算离散傅里叶变换对的一种快速算法，可分为快速傅里叶变换(fast fourier transform,fft)和逆快速傅里叶变换(inverse fast fourier transform,ifft)。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数越多，fft/ifft算法计算量的节省就越显著。详见“程乾生.数字信号处理[m].北京大学出版社.2003”。
[0018]
定义6.三维后向投影算法
[0019]
三维后向投影算法bp(back-projection,bp)利用雷达平台的轨迹信息求出雷达平台与场景像素点的距离历史，然后通过遍历距离历史找出回波数据中匹配的回波数据，再进行相位补偿并相干累加，然后把复值结果反向投影到三维图像空间中完成三维成像处理。详见“师君.双基地sar与线阵sar原理及成像技术研究[d].电子科技大学博士论文.2009”。
[0020]
定义7.方位向、高度向、距离向
[0021]
将雷达平台横向运动的方向叫做方位向，纵向运动的方向叫做高度向，垂直于前两者的方向叫做距离向。
[0022]
定义8.l1范数正则化项
[0023]
l1范数正则化项是一种目标先验信息约束项，它由目标场景图像x的l1范数构成，可表示为||x||1＝∑|x|，其中x为x的元素，∑(
·
)为求和累加函数。具体通过对图像像素x
的绝对值累加计算。该正则项常用于sar稀疏成像方法，用于增强所获得目标场景图像的稀疏性。详见“wang y,zhang x,zhan x,et al.an rcs measurement method using sparse imaging based 3d sar complex image[j].ieee antennas and wireless propagation letters,2021”。
[0024]
定义9.柯西正则化项
[0025]
柯西正则化项是一种目标先验信息约束项，可表示为p(x)＝∑ln(γ/γ2 |x|)，其中γ为尺度参数，ln(
·
)为取自然对数函数，∑(
·
)为求和累加函数。相较为定义8中的l1范数正则化项，柯西正则化项在零值处呈现出较温和的稀疏性。详见“o,mayo p,achim a.convergence guarantees for non-convex optimisation with cauchy-based penalties[j].ieee transactions on signal processing,2020,68:6159-6170”。
[0026]
定义10.标准的柯西正则化近端求解算子方法
[0027]
柯西正则化近端求解算子用于求解基于柯西正则化的求逆问题，该问题形式可表示为其中λ为正则项权重参数，为l2范数的平方，通过取绝对值平方计算。求解算子为该问题的最优解，即该算子具体通过求解一个对应的三次方程(u
3-xu2 (γ2 2u)u-xγ2＝0)实现。标准的柯西正则化近端求解算子方法详见“o,mayo p,achim a.convergence guarantees for non-convex optimisation with cauchy-based penalties[j].ieee transactions on signal processing,2020,68:6159-6170”。
[0028]
本发明提供的一种基于柯西正则化的阵列三维sar稀疏成像方法，其特征是它包括以下步骤：
[0029]
步骤1.初始化相关参数
[0030]
采用标准步进频率(sf)信号形式以及阵列三维sar模型，并初始化如下参数：光在空气中传播速度，记为c；取自然指数函数，记为exp(
·
)；虚数单位，记为j；圆周率，记为π；步进频率信号中心频率，记为fc；步进频率序号，记为n＝1，2，...，n，其中n表示频点总数；方位向阵元序号，记为l＝1，2，...，l，其中l表示方位向阵元总数；方位向阵元间隔，记为d
l
；方位向合成孔径，记为da；高度向阵元序号，记为m＝1，2，...，m，其中m表示高度向阵元总数；高度向阵元间隔，记为dm；高度向合成孔径，记为de；频率升采样倍数，记为k；参考距离，记为r0；三维目标场景空间，记为v
′
；三维图像像素位置向量，记为r
′
＝[x y z]，且满足r
′
∈v
′
；目标原始回波数据矩阵为s
l
×m×n，其中第一维表示方位向响应，第二维表示高度向响应，第三维表示频率响应；γ为柯西正则项尺度参数；λ为柯西正则项权重平衡参数。
[0031]
步骤2.对目标原始回波数据矩阵进行脉冲压缩处理，得到距离向数据矩阵
[0032]
以步骤1中的目标原始回波数据矩阵s
l
×m×n作为输入，采用标准的脉冲压缩方法对s
l
×m×n中的第三维频率响应进行压缩，得到距离向数据矩阵p
l
×m×n。
[0033]
步骤3.对脉冲压缩后的数据矩阵进行频率升采样
[0034]
以步骤2得到的距离向数据矩阵p
l
×m×n和步骤1中初始化的频率升采样倍数k作为输入，进行k倍频域升采样处理：
[0035]
步骤3.1.取出步骤2中距离向数据矩阵p
l
×m×n的第l行和第m列，得到1
×1×
n大小的向量p
lm
，其中n为步骤1中初始化的频点总数。
[0036]
步骤3.2.采用标准的快速傅里叶变换(fft)处理向量p
lm
，得到向量f
lm
。
[0037]
步骤3.3.从向量f
lm
的n/2 1位置开始插入(k-1)
·
n个零元素，得到f
′
lm
＝[f
lm
(1，2，
…
n/2 1) o
(k-1)
·
n f
lm
(n/2 2，
…
，n)]，其中f
lm
(1，2，
…
n/2 1)表示f
′
lm
中的前n/2 1个元素，f
lm
(n/2 2，
…
，n)表示f
′
lm
中的后n/2 1个元素，o
(k-1)
·n表示插入的(k-1)
·
n个零元素。
[0038]
步骤3.4.采用标准的逆快速傅里叶变换(ifft)处理向量f
′
lm
，得到向量p
′
lm
。
[0039]
步骤3.5.将向量p
′
lm
存放至一个l
×m×
(kn)矩阵的第l行和第m列，得到升采样后的数据矩阵p
′
l
×m×
(kn)
。
[0040]
步骤4.计算成像相位补偿因子
[0041]
以步骤1中初始化的光在空气中传播速度c、自然指数函数exp(
·
)、虚数单位j、圆周率π、方位向阵元序号l、方位向阵元间隔d
l
、方位向最大孔径da、高度向阵元序号m、高度向阵元间隔dm、高度向最大孔径de、信号中心频率fc、三维目标场景空间v
′
、三维图像像素位置向量r
′
和参考距离r0作为输入，计算成像相位补偿因子：
[0042]
步骤4.1.采用公式r＝[r
0 l
·dl-da/2 m
·dm-de/2]计算得到当前阵元位置向量，记为r；然后，采用公式rr＝||r-r
′
||
２
计算得到三维目标场景空间v
′
中的距离历史，记为rr，其中||
·
||2表示l2范数，通过元素取绝对值的平方累加后开根号计算，r0为参考距离，l为方位向阵元序号，d
l
为方位向阵元间隔，da为方位向最大孔径，m为高度向阵元序号，dm为高度向阵元间隔，de为高度向最大孔径。
[0043]
步骤4.2.遍历三维目标场景空间v
′
，计算成像相位补偿因子k
nn
＝exp(-j4πrrfc/c)，记为k
nn
，其中exp(
·
)为取自然指数函数，j为虚数单位，π为圆周率，rr为步骤4.1中所得距离历史，fc为信号中心频率，c为光在空气中传播速度。
[0044]
步骤5.利用三维bp算法获取三维初始图像
[0045]
以步骤3中得到的数据矩阵p
′
l
×m×n和步骤4中得到的成像相位补偿因子k
nn
作为输入，采用标准的三维bp算法进行成像处理，得到被测目标场景的三维初始图像，记为y
mf
。
[0046]
步骤6.利用基于柯西正则化的近端求解算子实现三维增强成像
[0047]
步骤6.1.以步骤5中得到三维初始图像y
mf
、步骤1中定义的柯西正则项尺度参数γ和权重平衡参数λ作为输入，构建基于柯西正则化的三维稀疏成像方程，基于柯西正则化的三维稀疏成像方程如下：
[0048][0049]
其中x为待求三维稀疏图像，x为图像中的单个像素值，y
mf
为初始图像y
mf
中的单个像素值。为取斐波那契范数的平方，通过对每个像素值取平方，累加后计算。ln(
·
)为取自然对数函数。argmin
x
{
·
}为以x为自变量的取最小值函数。∑(
·
)为求和累加函数。表示l2范数的平方，通过像素值取平方计算。
[0050]
步骤6.2.采用标准的柯西正则化近端求解算子方法求解步骤6.1中基于柯西正则化的三维稀疏成像方程，得到三维稀疏图像x。
[0051]
x即为所提基于柯西正则化的阵列三维sar稀疏成像方法的最终成像结果。
[0052]
本发明的创新点与优点在于：采用了一种不同于基于l1范数正则化的稀疏成像方
法，通过在成像模型中加入更合适的柯西正则化项，来实现阵列三维sar成像。与基于匹配滤波的成像方法比，成像结果分辨率更高、杂波与噪声更少，成像质量更高，与基于l1范数正则化项的稀疏成像方法相比，本发明方法能够在提高图像质量的同时，更好得保留目标能量，保留目标的细节轮廓。成像质量更优，有利于后续的目标分类、识别等应用。
附图说明：
[0053]
图1是本发明流程图。
[0054]
图2是本发明仿真实验验证结果。
具体实施方式
[0055]
本发明主要采用仿真实验的方法进行验证，所有步骤、结论都在数学计算软件matlab2019b上验证正确。具体实施步骤如下：
[0056]
步骤1.初始化相关参数
[0057]
采用标准步进频率(sf)信号形式以及阵列三维sar模型，并初始化如下参数：光在空气中传播速度，记为c＝3
×
108m/s；取自然指数函数，记为exp(
·
)；虚数单位，记为j；圆周率，记为π＝3.14；步进频率信号中心频率，记为fc＝37.5ghz；步进频率序号，记为n＝1，2，...，256，频点总数为256；方位向阵元序号，记为l＝1，2，...，76，方位向阵元总数为76；方位向阵元间隔d
l
＝0.04m；方位向合成孔径为da＝3m；高度向阵元序号，记为m＝1，2，...，76，高度向阵元总数为76；高度向阵元间隔dm＝0.04m；高度向合成孔径为de＝3m；频率升采样倍数为k＝8；参考距离为r0＝30m，目标场景空间，记为v
′
＝[
±
10
ꢀ±
10
ꢀ±
20]m；三维图像像素位置向量，记为r
′
＝[x y z]∈[
±
10
ꢀ±
10
ꢀ±
20]；仿真产生的目标原始回波数据矩阵为s
76
×
76
×
101
，其中第一维表示方位向响应，第二维表示高度向响应，第三维表示频率响应；柯西正则项尺度参数为γ＝4.4；柯西正则项权重平衡参数λ＝1。
[0058]
步骤2.对目标原始回波数据矩阵进行脉冲压缩处理，得到距离向数据矩阵
[0059]
以步骤1中的目标原始回波数据矩阵s
76
×
76
×
101
作为输入，采用标准的脉冲压缩方法对s
76
×
76
×
101
中的第三维频率响应进行压缩，得到距离向数据矩阵p
76
×
76
×
101
。
[0060]
步骤3.对脉冲压缩后的数据矩阵进行频率升采样
[0061]
以步骤2得到的距离向数据矩阵p
76
×
76
×
101
和步骤1中初始化的频率升采样倍数k＝8作为输入，进行8倍频域升采样处理：
[0062]
步骤3.1.取出步骤2中距离向数据矩阵p
76
×
76
×
256
的第l行和第m列，得到1
×1×
256大小的向量p
lm
。
[0063]
步骤3.2.采用标准的快速傅里叶变换(fft)处理向量p
lm
，得到向量f
lm
。
[0064]
步骤3.3.从向量f
lm
的129位置开始插入1792个零元素，得到f
′
lm
＝[f
lm
(1，2，
…
129) o
1792 f
lm
(130，
…
，256)]，其中f
lm
(1，2，
…
129)表示f
′
lm
中的前129个元素，f
lm
(130，
…
，256)表示f
′
lm
中的后127个元素，o
1792
表示插入的1792个零元素。
[0065]
步骤3.4.采用标准的逆快速傅里叶变换(ifft)处理向量f
′
lm
，得到向量p
′
lm
。
[0066]
步骤3.5.将向量p
′
lm
存放至一个101
×
101
×
1792矩阵的第l行和第m列，得到升采样后的数据矩阵p
′
101
×
101
×
1792
。
[0067]
步骤4.计算成像相位补偿因子
[0068]
以步骤1中初始化的光在空气中传播速度c＝3
×
108m/s、自然指数函数exp(
·
)、虚数单位j、圆周率π＝3.14、方位向阵元序号l＝1，2，...，76、方位向阵元间隔d
l
＝0.04m、方位向最大孔径da＝3m、高度向阵元序号m＝1，2，...，76、高度向阵元间隔d
m＝0.04m
、高度向最大孔径de＝3m、信号中心频率fc＝37.5ghz、三维目标场景空间v
′
＝[
±
10
ꢀ±
10
ꢀ±
20]m、三维图像像素位置向量r
′
和参考距离r0＝30m作为输入，计算成像相位补偿因子：
[0069]
步骤4.1.采用公式r＝[30 0.04l-1.5 0.04m-1.5]计算得到当前阵元位置向量，记为r；然后，采用公式rr＝||r-r
′
||2计算得到三维目标场景空间v
′
中的距离历史，记为rr，其中||
·
||2表示l2范数，通过元素取绝对值的平方累加后开根号计算。
[0070]
步骤4.2.遍历三维目标场景空间v
′
，计算成像相位补偿因子k
nn
＝exp(-j4πrrfc/c)，记为k
nn
。其中exp(
·
)为取自然指数函数，j为虚数单位，π＝3.14，rr为步骤4.1中所得距离历史，fc＝37.5ghz，c＝3
×
108m/s。
[0071]
步骤5.利用三维bp算法获取三维初始图像
[0072]
以步骤3中得到的数据矩阵p
′
101
×
101
×
1792
和步骤4中得到的成像相位补偿因子k
nn
作为输入，采用标准的三维bp算法进行成像处理，得到被测目标场景的三维初始图像，记为y
mf
。
[0073]
步骤6.利用基于柯西正则化的近端求解算子实现三维增强成像
[0074]
步骤6.1.以步骤5中得到三维初始图像y
mf
、步骤1中定义的柯西正则项尺度参数γ＝4.4和权重平衡参数λ＝1作为输入，构建基于柯西正则化的三维稀疏成像方程，基于柯西正则化的三维稀疏成像方程如下：
[0075][0076]
其中x为待求三维稀疏图像，x为图像中的单个像素值，y
mf
为初始图像y
mf
中的单个像素值。为取斐波那契范数的平方，通过对每个像素值取平方，累加后计算。ln(
·
)为取自然对数函数。argmin
x
{
·
}为以x为自变量的取最小值函数。∑(
·
)为求和累加函数。表示l2范数的平方，通过像素值取平方计算。
[0077]
步骤6.2.采用标准的柯西正则化近端求解算子求解步骤6.1中构建的基于柯西正则化的三维稀疏成像方程，得到三维稀疏图像x。x即为所提基于柯西正则化的阵列三维sar稀疏成像方法的最终成像结果。计算机仿真结果为图2所示。
[0078]
经过计算机仿真结果表明，本发明将柯西正则化引入阵列三维sar稀疏成像中，可获得高质量的成像结果，与基于匹配滤波的成像方法比，本发明方法能够提高分辨率、抑制杂波与噪声，提高成像质量，与基于l1范数正则化项的稀疏成像方法相比，本发明方法能够在提高图像质量的同时，更好得保留目标能量，保留目标的细节轮廓。