一种太赫兹雷达超高分辨成像方法与流程

本发明属于雷达信号处理技术领域，具体涉及一种太赫兹雷达超高分辨成像方法。

背景技术：

对于散射能量具有各向异性的目标成像时，如果直接利用宽孔径数据相干积累成像，由于散射体的强散射能量仅集中在小部分孔径内会导致目标的信噪比降低进而使目标的棱角模糊。对散射能量具有各向异性目标成像的传统成像方法是：将宽孔径划分为多个固定中心等宽度的子孔径，并基于glrt成像方法非相关联合这些子孔径图像形成融合图像，并且子孔径宽度越小孔径内细节越清晰。由于太赫兹雷达具有目标散射能量各向异性更加剧烈的缺点，因此太赫兹雷达尤其是用于探测多棱角、多边沿的人造目标的太赫兹雷达应当采用划分子孔径的方法来进行成像。为了更好地应用划分子孔径的成像方法，研究者针对传统划分子孔径的成像方法划分的子孔径没有合适的子孔径中心以及适应不同sar图像的最佳子孔径宽度，无法保证同一个强散射体的强散射能量回波数据被划分到同一个子孔径中的缺点做出了改进。

2009年，中国微波成像技术国家级重点实验室的q.li等人采用了对全孔径数据分别采用一个滤波器生成子孔径数据以及采用与散射体的散射特性相匹配的滤波器组生成子孔径数据的方法来提高子孔径图像包含的探测目标细节，并采用基于glrt成像方法非相关联合成像的方式融合子孔径图像成像。

2018年，电子科技大学的刘通采用基于子孔径能量的自适应子孔径划分方法将全孔径划分为多个自适应中心的等宽度子孔径，并采用基于glrt成像方法非相关联合成像的方式融合子孔径图像成像完成了太赫兹雷达的孔径自适应高分辨成像。

然而该技术以子孔径能量为依据划分子孔径仅能实现太赫兹雷达对人造目标的高分辨成像，且成像效果不稳定。原因在于该技术时常能高分辨成像的关键在于将自适应得到的子孔径统一为相同宽度，然而统一子孔径宽度不是太赫兹雷达对人造目标高分辨成像真正的关键原理，该技术能实现高分辨成像是因为该技术统一子孔径宽度的过程中时常偶然实现关键孔径的子孔径相邻重叠。该技术实现依靠一定的偶然性因此成像效果不稳定。同时该技术只对能量较高的子孔径成像其实丢失了挺多的成像细节，无法实现超高分辨成像。并且子孔径能量来划分子孔径，之后再统一子孔径宽度没有将子孔径边界设置在最合理位置，导致该方法没有完全解决小孔径子孔径划分成像的旁瓣散焦问题，成像效果有一定模糊。

技术实现要素：

本发明的目的，就是针对上述存在的问题及不足，提供一种尤其适用于对散射能量具有各向异性目标太赫兹雷达csar成像的超高分辨成像方法。

本发明的目的是通过以下技术方案来实现的：一种太赫兹雷达超高分辨成像方法，其特征在于划分子孔径的过程中采用了基于子孔径互相关系数结合子孔径能量的算法来确定子孔径边界，并且在划分子孔径的算法中采用了将子孔径依次重叠一部分的技术，包括如下步骤：

(1)使用太赫兹雷达系统对需要成像的目标进行csar数据采集，并将原始数据处理为包含一个圆周孔径回波信息的矩阵文件，其中行表示快时间，列表示慢时间也就是角度：

雷达系统数据采集后产生一个矩阵，该矩阵的行数为快时间采样数的两倍，其中前一半行存储快时间回波信号的实部，后一半行存储快时间回波信号的虚部。该矩阵的列数为慢时间采样数。

对矩阵进行处理：从该矩阵中截取一个圆周孔径对应的慢时间数，该慢时间数m可以根据雷达系统旋转的角速度ω和雷达系统的脉冲重复频率prf计算得到。从截取到的矩阵中截取表示快时间回波信号实部的行，并对每一矩阵单元加上对应回波信号的虚部，从而使截取到的矩阵行表示快时间回波复信号。将截取到的矩阵作为用于雷达成像或目标识别的通用回波数据。

(2)开始利用回波数据划分子孔径，首先找到圆周孔径中包含强散射能量的角度：

首先计算目标散射能量与角度的函数：i(θ,n)和q(θ,n)分别表示回波信号的实部与虚部，θ表示慢时间方位角度，n，1<n<n表示快时间采样，n为快时间采样数即行数，对于成像场景中仅存在唯一目标的实测数据，目标散射能量与角度的函数表示为：

然后计算目标散射能量函数的包络：

phe(θ)＝interp{lmv[ph(θ)]}(公式3)

其中lmv[·]表示提取局部极大值操作，interp[·]表示插值操作。

接着计算强散射能量门限：

其中a表示强散射能量门限与平均散射能量的倍率，m为.mat文件的慢时间数。

最后挑出具有强散射能量的角度：创建一个可变长数组indivisible来存储强散射能量的角度。遍历phe(θ)，如果phe(θ)>et，则：

indivisible(ii)＝θ

其中ii表示该慢时间在数组中的下标。

(3)根据慢时间的相邻互相关系数划分子孔径区间：

首先计算目标互相关系数与角度的函数：依次计算每个角度与相邻角度的互相关系数

其中，x、y分别为.mat文件矩阵相邻的两个列，cov(x,y)为x、y的协方差，var[x]为x的方差，var[y]为y的方差。然后计算每个角度与相邻角度互相关系数的平均值作为函数值

其中由于雷达系统采集的是一个圆周孔径的csar回波数据，第一个慢时间和最后一个慢时间视为相邻。

然后对目标互相关系数函数求相邻平均：

其中2×long+1表示求相邻平均的区间长度。

接着计算目标互相关系数函数的包络：

mue(θ)＝interp{lmv[mu(θ)]}(公式8)

其中lmv[·]表示提取局部极大值操作，interp[·]表示插值操作。

搜索适合作为子孔径边界点的角度：

[xpeakmin1,ypeakmin1]＝getpeakmin(mue)；

[xpeakmin2,ypeakmin2]＝getpeakmin(ypeakmin1)；

thetamin＝xpeakmin1(xpeakmin2)；

其中[a,b]＝getpeakmin(c)表示一个从函数c中分别找到极小值的横坐标数组a和纵坐标数组b的方法。xpeakmin1,ypeakmin1分别表示函数mue极小值的横纵坐标数组，xpeakmin2,ypeakmin2分别表示纵坐标数组ypeakmin1中极小值的序号数组和纵坐标数组，thetamin表示对函数mue极小值再取极小值的横坐标数组。

划分子孔径区间:以数组thetamin的元素作为划分子孔径的边界点，将圆周孔径划分为多个子孔径区间。得到子孔径区间左边界数组thetal以及子孔径区间右边界数组thetar；

(4)根据步骤s2中获得的圆周孔径中包含强散射能量的角度，从步骤s3划分出的子孔径区间中标记出三类子孔径区间：如果子孔径区间中包含强散射能量角度则该区间为第一类子孔径区间，如果子孔径区间位于两个第一类子孔径区间之间并且左右第一类子孔径区间之间的间隔小于子孔径宽度门限则中间的这些子孔径区间为第二类窄子孔径区间，其余子孔径区间为第三类子孔径区间；所述子孔径宽度门限为预设值：

创建一个数组flag，元素分别对应子孔径区间，初始值为0。创建一个可变长数组flagindex用于存储第一类数组序号。如果子孔径区间属于第一类子孔径区间则将数组flag的对应元素标记为1并且将该子孔径区间序号添加到数组flagindex。如果子孔径区间属于第二类子孔径区间则将数组flag的对应元素标记为2。

(5)对第一类子孔径区间拓宽得到第一类成像子孔径：

遍历数组flagindex，对该类子孔径区间依次执行拓宽算法

用变量thetaltemp、thetartemp存储拓宽中的子孔径左右边界，初始值为第ii个子孔径的左右边界：

forjj＝ii+1:lengthof(flagindex)

if(第jj个子孔径与拓宽中的子孔径不相邻)

结束子孔径拓宽；

elseif(子孔径边界回波能量大于门限值||拓宽中的子孔径区间宽度小于门限值)

将第jj个子孔径所在的子孔径区间添加到拓宽中的子孔径中；

end

其中thetaltemp、thetartemp在拓宽算法结束后表示拓宽后子孔径的左右边界，lengthof(·)表示求数组的长度的方法。

(6)对第三类子孔径区间拓宽得到第二类成像子孔径：

遍历数组flag，对其中的第一类子孔径区间依次执行拓宽算法

用变量thetaltemp、thetartemp存储拓宽中的子孔径左右边界，初始值为第ii个子孔径的左右边界

forjj＝ii+1:lengthof(flag)

if(拓宽中的子孔径区间宽度小于门限值&&第jj个子孔径是第一类子孔径区间)

将第jj个子孔径所在的子孔径区间添加到拓宽中的子孔径中；

else

结束子孔径拓宽；

end

其中thetaltemp、thetartemp在拓宽算法结束后表示拓宽后子孔径的左右边界，lengthof(·)表示求数组的长度的方法。

(7)对得到的两类成像子孔径采用基于glrt成像方法非相关联合成像的方式融合子孔径图像成像：

分别对两类成像子孔径使用pfa算法成像，并根据子图像所使用的子孔径角度将子图像旋转至同一角度对齐。然后根据子图像形成融合图像，融合图像的每一个像素值均为所有子图像该像素值的最大值。

本发明与现有技术相比，具有的有益效果是：

(1)采用了基于子孔径互相关系数结合子孔径能量的算法来确定子孔径边界角度，更好地消除了小孔径子孔径联合成像时的旁瓣散焦问题，成像效果更加清晰。并且相比现有技术成像效果消除旁瓣散焦的能力更加稳定。

(2)采用了子孔径依次重叠一部分的子孔径划分技术，解决了复杂人造目标进行小孔径的子孔径联合成像时细节回波被子孔径边界截断严重，导致太赫兹雷达对该类目标csar成像细节损失严重的问题。实现了对该类目标成像细节接近完整的效果。

(3)采用了寻找相对优秀的子孔径边界点而不是强散射能量孔径的方式划分子孔径，基本覆盖了全孔径。在增加很少运算量的同时，子孔径覆盖了目标的更多细节，对复杂人造目标的成像更加完整。

(4)前二者的影响下，对复杂人造目标可以进行孔径更小的子孔径联合成像，而子孔径联合成像的子孔径越小孔径内的细节就越清晰。因此成像的细节可以更清晰。

附图说明

图1是本发明的方法的流程图；

图2是成像实验场景示意图；

图3是对成像实验场景不划分子孔径，对圆周孔径回波数据pfa算法成像效果图；

图4是对成像实验场景使用现有成像方法，利用子孔径能量自适应划分子孔径并成像的效果图；

图5是对成像实验场景使用本发明的成像方法的成像效果图；

图6是本发明的成像方法与现有成像方法成像效果细节对比图。

具体实施方式

下面结合附图对本发明进行详细的描述

参照附图1，本发明的具体实施步骤如下：

(1)实测数据，将雷达系统得到的原始数据处理为通用.mat矩阵文件

雷达系统数据采集后产生一个用于计算机处理的.dat文件，该文件的内容是一个矩阵。该矩阵的行数为快时间采样数的两倍，其中前一半行存储快时间回波信号的实部，后一半行存储快时间回波信号的虚部。该矩阵的列数为慢时间采样数。

对该.dat文件进行处理：从该文件矩阵中截取一个圆周孔径对应的慢时间数，该慢时间数m可以根据雷达系统旋转的角速度ω和雷达系统的脉冲重复频率prf计算得到。从截取到的矩阵中截取表示快时间回波信号实部的行，并对每一矩阵单元加上对应回波信号的虚部，从而使截取到的矩阵行表示快时间回波复信号。将截取到的矩阵作为.mat文件保存，形成可以用于雷达成像或目标识别的通用回波数据。

(2)找到圆周孔径中包含强散射能量的角度

首先计算目标散射能量与角度的函数：i(θ,n)和q(θ,n)分别表示回波信号的实部与虚部，θ表示慢时间方位角度，n，1<n<n表示快时间采样，对于成像场景中仅存在唯一目标的实测数据，目标散射能量与角度的函数可以近似表示为：

然后计算目标散射能量函数的包络：

phe(θ)＝interp{lmv[ph(θ)]}(公式3)

其中lmv[·]表示提取局部极大值操作，interp[·]表示插值操作。

接着计算强散射能量门限：

其中a表示强散射能量门限与平均散射能量的倍率，m为.mat文件的慢时间数。

最后挑出具有强散射能量的角度：创建一个可变长数组indivisible来存储强散射能量的角度。遍历phe(θ)，如果phe(θ)>et，则：

indivisible(ii)＝θ

其中ii表示该慢时间在数组中的下标。

(3)根据慢时间的相邻互相关系数划分子孔径区间

首先计算目标互相关系数与角度的函数：依次计算每个角度与相邻角度的互相关系数

其中，x、y分别为.mat文件矩阵相邻的两个列，cov(x,y)为x、y的协方差，var[x]为x的方差，var[y]为y的方差。然后计算每个角度与相邻角度互相关系数的平均值作为函数值

其中由于雷达系统采集的是一个圆周孔径的csar回波数据，第一个慢时间和最后一个慢时间视为相邻。

然后对目标互相关系数函数求相邻平均：

其中2×long+1表示求相邻平均的区间长度。

接着计算目标互相关系数函数的包络：

mue(θ)＝interp{lmv[mu(θ)]}(公式8)

其中lmv[·]表示提取局部极大值操作，interp[·]表示插值操作。

搜索适合作为子孔径边界点的角度：

[xpeakmin1,ypeakmin1]＝getpeakmin(mue)；

[xpeakmin2,ypeakmin2]＝getpeakmin(ypeakmin1)；

thetamin＝xpeakmin1(xpeakmin2)；

其中[a,b]＝getpeakmin(c)表示一个从函数c中分别找到极小值的横坐标数组a和纵坐标数组b的方法。

划分子孔径区间:以数组thetamin的元素作为划分子孔径的边界点，将圆周孔径划分为多个子孔径区间。得到子孔径区间左边界数组thetal以及子孔径区间右边界数组thetar。

(4)从划分出的子孔径区间中标记出三类子孔径区间：包含总回波能量高于门限的角度的子孔径区间，第一类子孔径区间之间的间隙的窄子孔径区间，其他不受总回波能量高于门限的角度影响的子孔径区间

创建一个数组flag，元素分别对应子孔径区间，初始值为0。创建一个可变长数组flaginde用于存储第一类数组序号。如果子孔径区间属于第一类子孔径区间则将数组flag的对应元素标记为1并且将该子孔径区间序号添加到数组flagindex。如果子孔径区间属于第二类子孔径区间则将数组flag的对应元素标记为2。

(5)对第一类子孔径区间拓宽得到第一类成像子孔径

遍历数组flagindex，对该类子孔径区间依次执行拓宽算法

用变量thetaltemp、thetartemp存储拓宽中的子孔径左右边界，初始值为第ii个子孔径的左右边界：

forjj＝ii+1:lengthof(flagindex)

if(第jj个子孔径与拓宽中的子孔径不相邻)

结束子孔径拓宽；

elseif(子孔径边界回波能量大于门限值||拓宽中的子孔径区间宽度小于门限值)

将第jj个子孔径所在的子孔径区间添加到拓宽中的子孔径中；

end

其中thetaltemp、thetartemp在拓宽算法结束后表示拓宽后子孔径的左右边界，

lengthof(·)表示求数组的长度的方法。

(6)对第三类子孔径区间拓宽得到第二类成像子孔径：

遍历数组flag，对其中的第一类子孔径区间依次执行拓宽算法

用变量thetaltemp、thetaltemp存储拓宽中的子孔径左右边界，初始值为第ii个子孔径的左右边界

forjj＝ii+1:lengthof(flag)

if(拓宽中的子孔径区间宽度小于门限值&&第jj个子孔径是第一类子孔径区间)

将第jj个子孔径所在的子孔径区间添加到拓宽中的子孔径中；

else

结束子孔径拓宽；

end

其中thetaltemp、thetartemp在拓宽算法结束后表示拓宽后子孔径的左右边界，lengthof(·)表示求数组的长度的方。

(7)对得到的两类成像子孔径采用基于glrt成像方法非相关联合成像的方式融合子孔径图像成像

分别对两类成像子孔径使用pfa算法成像，并根据子图像所使用的子孔径角度将子图像旋转至同一角度对齐。然后根据子图像形成融合图像，融合图像的每一个像素值均为所有子图像该像素值的最大值。

该成像方法中存在一些最佳值根据回波特性变化的参数：

其中局部最大值间隔最佳值区间为代表孔径0.02°-0.08°对应的慢时间，以能去除函数毛刺为准。

强散射能量门限与平均散射能量的倍率a最佳值区间为0.5-5倍，以能挑出作为子孔径边界会产生较为严重旁瓣散焦的角度为准。

求相邻平均的区间长度2×long+1最佳值区间为代表孔径0.06°-0.20°对应的慢时间，以保证后续步骤两次求极小值的结果为函数包络总体趋势的极小值为准。

拓宽子孔径时使用的能量门限最佳值应该小于0.15倍最大能量大于0.05倍最大能量，以能消除强散射能量的旁瓣散焦而又不至于导致大部分该类子孔径区间都被合并形成非常大的孔径为准。

拓宽子孔径时使用的角度门限最佳值区间为3°-12°，以单个子孔径宽度足够容纳单一细节回波为准，太赫兹低频段最佳值在4°附近。

下面结合实测数据对本发明的效果作进一步说明。

利用工作带宽28.8ghz的0.34thz雷达系统采集实测数据。将j15飞机模型作为csar成像目标，其尺寸为44.5cm×30cm。飞机模型与真实飞机的尺寸比例为1：50。目标放置在距离雷达天线4.02米处的泡沫平台上，与雷达天线的相对高度为0.42m。雷达天线到泡沫平台中心的斜视角近视为6°。雷达角速度为4°/s，载波322.02ghz，带宽28.8ghz，脉冲宽度0.3ms，快时间采样频率1.5625mhz，脉冲重复频率1000hz。实验场景示意图见附图2。处理后的回波数据.mat文件矩阵包含512行、90000列。其中行是快时间，列是慢时间也就是角度。

对该回波数据不划分子孔径，使用圆周孔径回波数据pfa算法成像，成像效果见附图3。

对该回波数据使用基于子孔径能量的自适应子孔径划分方法成像，其中能量门限设置为最大能量的0.01，子孔径宽度门限设置为1°，成像效果见附图4。

对该回波数据使用本发明的成像方法进行成像，其中局部最大值间隔设置为0.04°，相邻平均区间宽度设置为0.124°，强散射能量门限设置为1倍目标散射能量的平均值，能量门限设置为最大能量的0.08倍，子孔径宽度门限设置为4°(本成像方法门限的用途和现有成像方法不一样)，成像效果见附图5。

对本发明的成像方法效果图与现有自适应子孔径划分方法成像效果图进行细节对比分析，见附图6，其中左边是本发明的成像方法效果图的细节，右边是现有自适应子孔径划分方法成像效果图的对应位置。