一种噪声环境下的雷达目标高分辨率距离像鲁棒识别方法

李玮杰; 杨威; 黎湘; 刘永祥

doi:10.12000/JR19093

一种噪声环境下的雷达目标高分辨率距离像鲁棒识别方法

DOI: 10.12000/JR19093

国防科技大学电子科学学院长沙 410073

基金项目: 国家自然科学基金(61871384, 61401486)，湖南省自然科学基金(2017JJ3367)，上海航天科技创新基金(SAST2017-048)

详细信息

作者简介:
李玮杰(1998–)，男，湖南株洲人，国防科技大学电子科学学院在读博士生，研究方向为雷达目标识别，机器学习，深度学习。E-mail: lwj2150508321@sina.cn

杨　威(1985–)，男，江西丰城人，博士，国防科技大学电子科学学院讲师，研究方向包括信息融合、多目标跟踪、融合目标识别。E-mail: yw850716@sina.com

黎　湘(1967–)，男，湖南浏阳人，教授，长江学者特聘教授，研究方向为目标探测识别与雷达成像等

刘永祥(1976–)，男，河北唐山人，博士，国防科技大学电子科学学院教授，博士生导师，研究方向为目标微动特性分析与识别。E-mail: lyx_bible@sina.com

通讯作者:
杨威 yw850716@sina.com

责任主编：陈渤 Corresponding Editor: CHEN Bo
中图分类号: TN957.51
计量
- 文章访问数: 3503
- HTML全文浏览量: 1303
- PDF下载量: 408
- 被引次数: 13
出版历程
- 收稿日期: 2019-10-23
- 修回日期: 2019-12-07
- 网络出版日期: 2020-08-28

Robust High Resolution Range Profile Recognition Method for Radar Targets in Noisy Environments

College of Electronic Science and Engineering, National University of Defense Technology, Changsha 410073, China

Funds: The National Natural Science Foundation of China (61871384, 61401486), The Natural Science Foundation of Hunan Province (2017JJ3367), The Innovation Foundation of Aerospace Science and Technology of Shanghai (SAST2017-048)

More Information

Corresponding author: YANG Wei, yw850716@sina.com

摘要

摘要: 随着深度学习技术被应用于雷达目标识别领域，其自动提取目标特征的特性大大提高了识别的准确率和鲁棒性，但噪声环境下的鲁棒性有待进一步研究。该文提出了一种在噪声环境下基于卷积神经网络(CNN)的雷达高分辨率距离像(HRRP)数据识别方法，通过增强训练集和使用残差块、inception结构和降噪自编码层增强网络结构，实现了在较宽信噪比范围下的较高识别率，其中在信噪比为0 dB的瑞利噪声条件下，识别率达到96.14%，并分析了网络结构和噪声类型对结果的影响。
- 目标识别 /
- 高分辨距离像 /
- 噪声环境 /
- 卷积神经网络
Abstract: With the application of deep learning technology in the radar target recognition field, the automatic extraction of the target feature greatly improves the accuracy and robustness of the recognition, but its robustness in noisy environments needs to be further investigated. This paper proposes a robust target recognition method for radar High Resolution Range Profile (HRRP) data based on Convolutional Neural Networks (CNN). By enhancing training set and using the residual block, inception structure, and denoising sparse autoencoder layer to enhance the network structure, a higher recognition rate is achieved in a wider SNR range, under the condition of 0 dB Rayleigh noise, the recognition rate reaches 96.14%, and the influence of the network structure and noise type on results is analyzed.
- Target recognition /
- High-Resolution Range Profile (HRRP) /
- Noise environment /
- Convolutional Neural Network (CNN)

HTML全文

1. 引言

作为一种主动式成像传感器，相比于红外、光学等被动式传感器，合成孔径雷达(Synthetic Aperture Radar, SAR)可以全天时、全天候地获取高分辨率图像^[1-4]。其具有成像模式丰富、覆盖范围广、分辨率高等特点，能够长期、动态、宏观地对海洋进行监测，目前已成为世界各国监测海洋目标的主流方法之一^[5-11]。

单通道SAR图像在舰船目标检测系统中应用极为广泛，一直是海洋监视领域的研究重点^[12-17]。随着SAR技术的快速发展，先进的SAR传感器如TerraSAR-X, COSMOS-SkyMed, RADARSAT-2, ALOS-PALSAR, Sentinel-1和高分三号等可以提供相比过去更高的分辨率。更高的分辨率使得海杂波背景更加地不均匀，因此对舰船目标检测方法也提出了更高的要求。SAR采用微波相干成像，因此单通道SAR图像的本质是复值的，这是其和普通光学红外遥感图像的显著不同之处和核心区别之一。然而，在完全发展的相干斑的假设条件下，以恒虚警率(Constant False Alarm Rate, CFAR)为代表的传统单通道SAR图像舰船目标检测方法主要研究幅度信息^[9,12-20]，认为相位服从均匀分布，不包含任何有效信息。同时，深度学习在传统图像目标检测识别领域取得了突飞猛进的发展，但是绝大多数研究者采用深度学习方法对SAR图像处理时也仅考虑幅度信息，将其作为普通的光学图像进行对待。

考虑到传统观点主要基于过去低分辨率SAR图像，其必须满足以下条件^[1]：(1)均匀媒质的一个分辨单元内有大量的散射单元；(2)斜距距离远大于雷达波长；(3)以雷达波长尺度衡量媒质表面非常“粗糙”；但这些并非完全适用于如今的高分辨率SAR图像^[21-24]，也就是说完全发展的相干斑的部分假设条件在高分辨率情形下并非严格成立。这主要是因为更高的分辨率通常意味着更小的分辨单元，同时一个分辨单元内包含的散射子也更少。对于舰船目标而言，其中一些散射子的强散射可能会对整个后向散射起到主导性的作用。单通道复值SAR图像既包括了幅度信息，也包括了不同于幅度的描述目标状态的相位信息，因此亟需对单通道SAR图像的相位或者复值信息进行研究，这将能在复值层面更完整地描述舰船目标和海杂波的区别，也有助于更深刻地理解SAR图像不同于普通光学遥感图像的复值本质，对提升单通道SAR图像舰船目标检测性能具有重要意义。

此外，近十年来，复值信号处理理论和方法得到了快速发展，Ollila等人^[25-27]、Eriksson等人^[28,29]和Novey等人^[30-32]的先驱工作进一步揭示了复值信号处理和传统实信号处理的区别，Schreier和Scharf二人^[33]对复值数据统计信号处理方法进一步进行了总结。在此基础上，El-Darymli等人^[23,24]对单通道复值SAR图像目标分类问题进行了研究，Wu等人^[34,35]对单通道复值SAR图像的地物分类问题进行了研究，Soccorsi等人^[36,37]采用高斯-马尔科夫随机场模型对单视复值(Single Look Complex, SLC)数据进行了建模研究。这些初步结果证明了复值信息对单通道SAR图像解译具有重要作用。

本文面向舰船目标检测，将单通道复值SAR图像统计建模方法分为幅度统计建模、相位统计建模和复值统计建模3个部分，在对单通道SAR图像幅度统计建模方法进行简要介绍的基础上，首先重点对相位统计建模方法进行研究，引入了相位循环统计量并介绍了适用于相位统计建模的Von Mises分布，并且论证了所提邻域相位方向差特征对区分舰船目标和海杂波的可行性。然后本文对基于复广义高斯分布的单通道SAR图像复值统计建模方法进行研究，详细介绍了其概念和生成过程，并且重点对其参数估计方法进行了研究，比较了基于最大似然法和提出的基于高阶矩的形状参数估计方法。同时本文还介绍了作者研究小组在基于复值统计信息的单通道SAR图像舰船目标检测的部分最新研究成果。最后本文展望了进一步研究方向。

2. 单通道SAR图像幅度统计建模方法

对于传统CFAR舰船目标检测算法而言，对海杂波幅度概率密度函数进行精确的描述至关重要。概率密度函数的精度会影响到检测结果的可靠性和稳健性。但是描述不同状态下的海杂波比较复杂，这是由于海面SAR成像通常受到许多因素的影响，比如海洋表面特性(海况、风速、风向等)、成像几何(入射角、方位角等)、雷达参数(分辨率、波长、极化方式等)以及外部干扰(射频干扰等)等^[2,38,39]。由于已有大量相关文献对幅度统计建模进行阐述^[8,9]，因此本节对此仅进行简要介绍。

当满足完全发展的相干斑条件时，传统观点认为可以采用中心极限定理对海杂波进行描述。因此，海杂波实部和虚部通常被认为服从高斯分布，单视幅度服从Rayleigh分布，多视幅度服从Nakagami分布(Gamma均方根分布)^[1]。当分辨率较低时，这些分布可以对SAR图像海杂波进行较好的拟合。但是当SAR图像分辨率提高时，海杂波呈现较强的拖尾特性^[40]，和Rayleigh分布或Nakagami分布存在明显不同。

为了描述海杂波的拖尾特性，K分布^[40]、对数正态^[41]和Weibull^[42]模型等被提出来。虽然对数正态和Weibull分布模型是经验的，但可以较好地描述Rayleigh分布无法描述的高分辨率SAR图像海杂波。此外，K分布是一种重要的分布模型，其将海洋回波认为是具有不相关事件的两个分量的调制结果。K分布可以通过假设分辨单元内散射单元数目服从负二项式分布得到，也可以利用服从Rayleigh分布的幅度和服从Gamma分布的纹理表征乘积模型得到^[1]。K分布是基于物理散射过程推导的，可以退化为Rayleigh分布，因而具有显著优势。但是K分布依然无法对一些高海况引起的极不均匀海洋杂波进行较好地描述。因此，许多学者尝试对K分布进行了改进。改进方法主要可以分为两类^[7]。其中一种观点认为相干斑不是完全发展的，在高分辨率条件下中心极限定理不再适用。从这点出发，提出了一些相关的分布模型。例如，当在均匀海杂波中存在一个强反射时，可以采用Nakagami-Rice分布^[43,44]。此外，重拖尾Rayleigh分布(Heavy-Tailed Rayleigh, HTR)^[45]以及最近提出的广义K分布(Generalized-K, GK)^[46]同样属于这一类型。其中，重拖尾Rayleigh分布认为海杂波实部和虚部服从零均值的对称 $\alpha$ 稳态分布^[47]，其可以描述比传统Rayleigh分布拖尾更长或者更像脉冲的分布；广义K分布是K分布的扩展，可以被用于描述单通道或者多通道SAR图像等多种不同情景^[11,46,48]。另一种改进K分布的思路认为依然可以采用完全发展的相干斑假设，但是需要加上一些额外的分量(例如热噪声或Rayleigh残留)进行校正。K+Noise(K+N)^[49], K+Rayleigh(K+R)^[50], Pareto+Noise(P+N)^[51], KA^[52]和KK^[53,54]等分布均是根据这种原理来改进传统的K分布，因此也能较好地描述海杂波中的尖峰或者热噪声等。在这些模型中，因为需要估计的参数太多，KA和KK分布在实际中较难应用。K和K+N分布都可以认为是K+R分布的特例。同时，相比于P+N分布，K+R分布可以更好地描述海杂波^[55,56]。此外，G⁰分布和广义Gamma分布也是在SAR图像应用中使用较多的模型，其精度相对K分布具有一定的优势^[57]。

在检测舰船目标的实际应用中，除了需要考虑统计分布模型的精度，还需要考虑更多因素^[7,8]：

第1点就是要考虑模型的复杂度以及在实际应用中的可行性。在前面提到的一些模型中，例如HTR, K+R以及P+N等，虽然理想情况下拟合精度更高，但是其参数估计方法复杂不稳定；同时即使准确求出所有参数，其阈值求解方法同样严重依赖于复杂的数值模型。这就使得它们虽然可用于某些SAR图像处理应用，例如滤波^[58]，但它们并不是舰船目标CFAR分布模型的理想选择。

第2点就是要考虑模型的适应性。这里的适应性不仅指对不同海况的描述能力，同时指对不同类型的SAR图像，例如对单视/多视图像和幅度/强度图像等的适应能力。K分布是最为广泛使用的分布模型之一，其中最主要的原因就是K分布不仅可以拟合单视幅度图像，还可以拟合多视幅度图像，且幅度K分布和强度K分布的参数之间存在明确的关系。

第3点就是要辩证看待基于物理模型和经验模型的分布。通常基于物理模型的分布可以更好地解释SAR图像中的一些现象。但是，由于海杂波复杂的电磁散射机理，物理模型也很难描述各种不同状态下的海杂波分布。另一方面，许多物理模型本身就基于经验模型，并没有很好地解释其中的原因。相反地，虽然部分经验模型没有严格的物理模型支撑，但却能够较好地拟合海杂波。这说明基于物理模型和经验模型的分布之间并没有严格的区分。因此在实际舰船目标检测应用中，只需模型能够较好地拟合真实杂波背景即可，而不需过多关注其建立的依据。

3. 单通道SAR图像相位统计建模方法

回波信号的相位一方面来源于电磁波的传播路程，一方面来源于多个散射点随机散射的综合结果。SAR图像最终记录的相位信息是幅角的主值，其在极化或干涉应用中具有重要的作用。但在文献中很少见到关于直接利用单通道SAR图像中的相位信息进行应用的报道。很少一部分学者在这方面进行了研究。早期，Rihaczek和Hershkowiz^[59,60]在这方面出版了两部著作，但是没有引起足够的重视。MIT林肯实验室的Jao等人^[61]提出了一种相干空间滤波器，但是只适用于检测直线目标。最近，德国DLR的Datcu教授等人^[62]重新提出了一系列对米级高分辨率SAR图像相位信息进行描述的构想。加拿大北部雷达实验室的El-Darymli等人^[23,24]对此进行了一些初步的研究。本节结合相位数据的特点，将引入循环数据统计量的方法对相位进行统计描述，然后介绍适用于相位建模的von Mises分布以及可以描述相位信息的邻域相位方向差特征。

3.1 相位循环统计量

相位是一种典型的循环数据(circular statistics)，其分布是一种循环分布(相位 $- \pi$ 和 $\pi$ 在圆上首尾相接在同一个点)。对这种数据的统计描述有别于传统的非循环数据。例如考虑两个角度，1°和359°，按传统方法计算会认为其均值直接为180°。但是稍加思考就会发现，将其均值认为是0°会更恰当。针对循环数据需要使用专门的统计量才能对其进行有效描述。在解译这些统计量时要注意将所有相位当作矢量，然后首尾相连，才能更好地理解^[63]。

设有一组相位观测量 ${x_1}$ , ${x_2}$ , ···, ${x_n}$ , ···, ${x_N}$ ，其中 ${x_n}$ 表示第 $n$ 个相位观测值， $N$ 是其总数，令

$\left. \begin{aligned} &{C_p} = \sum\limits_{n = 1}^N {\cos \left( {p{x_n}} \right)} ,\quad {{\bar C}_p} = \frac{{{C_p}}}{N}\\ &{S_p} = \sum\limits_{n = 1}^N {\sin \left( {p{x_n}} \right)} ,\quad {{\bar S}_p} = \frac{{{S_p}}}{N} \end{aligned}\right\}$

(1)

${R_p} = \sqrt {{C^2_p} + {S^2_p}} ,\quad {\bar R_p} = \frac{{{R_p}}}{N}$

(2)

${T_p} = \left\{ \begin{aligned} & {{\rm{tan}}{^{ - 1}}\left( {\frac{{{{\bar S}_p}}}{{{{\bar C}_p}}}} \right),\qquad\quad \;\;\;\;{{\bar S}_p} > 0, {{\bar C}_p} > 0} \\ & {{\rm{tan}}{^{ - 1}}\left( {\frac{{{{\bar S}_p}}}{{{{\bar C}_p}}}} \right) + {\rm{\pi }},\;\;\;\;\;\;\;\;\; {{\bar C}_p} < 0} \\ & {{\rm{tan}}{^{ - 1}}\left( {\frac{{{{\bar S}_p}}}{{{{\bar C}_p}}}} \right) + 2{\rm{\pi }},\;\;\;\;\;\;\; {{\bar S}_p} < 0,{{\bar C}_p} > 0} \end{aligned} \right.$

(3)

则根据Fisher^[64]以及Mardia和Jupp^[65]的论述，对应传统非循环数据的部分统计量可以给出以下一些相应定义：

(1) 循环均值(mean direction)

$\mu = {T_1}$

(4)

$\mu$ 是所有相位的平均值，代表了平均方向。

(2) 循环方差(circular variance)

$V = 1 - {\bar R_1}$

(5)

$V$ 在0～1之间变化，表示所有相位关于相位平均值的偏差。

(3) 循环标准差(circular standard deviation)

$v = \sqrt { - \ln {{\bar R}_1}}$

(6)

(4) 组合长度(resultant length)

$R = {R_1} = \sqrt {{C^2_1} + {S^2_1}}$

(7)

(5) 平均组合长度(mean resultant length)

$\bar R = {\bar R_1} = \frac{{{R_1}}}{N}$

(8)

${\bar R_1}$ 在0～1之间变化，若其接近1，则表示相位之间的变化较小。

(6) 循环离散度(circular dispersion)

$\delta = \frac{{1 - {T_2}}}{{2{{\bar R}_1}^2}}$

(9)

其通常被用于计算区间置信度。

(7) 循环偏度(circular skewness)

$S = \frac{{{{\bar R}_2}\sin ({T_2} - 2{T_1})}}{{{{(1 - {{\bar R}_1})}^{{3 / 2}}}}}$

(10)

(8) 循环峰度(circular kurtosis)

$K = \frac{{{{\bar R}_2}\cos ({T_2} - 2{T_1}) - {{\bar R}_1}^4}}{{{{(1 - {{\bar R}_1})}^2}}}$

(11)

图1给出了TerraSAR-X SM模式图像中2个舰船目标和2个海杂波切片的幅度图像示例。其分辨率约为3 m，极化方式为VV。图2是图1中图像中心黄色矩形区域对应的相位图像。由图2可以看出，舰船目标和海杂波的相位图像存在一定的差别，这主要体现在相位图像的纹理上。舰船目标相位图像的纹理存在一些更小的细节，其强散射点方位向和距离向上存在一些明显的条纹模式。但是这些差别由于条纹的相互重叠和随机性通常是难以描述的。图3是舰船目标和海杂波的相位直方图。由图3可以看出，舰船目标和海杂波的相位都近似服从均匀分布，相位平均值和方向存在很大的随机性，而相位长度都近似为0。表1进一步列出了舰船目标和海杂波的部分循环统计量结果。可以看出，舰船目标和海杂波的方差、标准差、偏度以及峰度等其它循环统计量之间不存在明显的差异。这是直接利用相位图像的一些统计信息进行目标识别的一个难点。

图 1 TerraSAR-X舰船目标和海杂波幅度图像示例

Figure 1. TerraSAR-X ship target and sea clutter amplitude images

下载: 全尺寸图片幻灯片

表 1 TerraSAR-X舰船目标和海杂波相位图像循环统计量结果

Table 1. Circular statistical results of TerraSAR-X ship target and sea clutter

循环统计量	均值	方差	标准差	平均长度	偏度	峰度
(a) 舰船目标1	2.21	0.95	2.46	0.05	0	–0.01
(b) 海杂波1	–2.72	0.96	2.49	0.04	–0.01	0.01
(c) 舰船目标2	–0.34	0.98	2.84	0.02	0.01	–0.01
(d) 海杂波2	–2.87	0.98	2.83	0.02	0.01	0

下载: 导出CSV

| 显示表格

图 2 TerraSAR-X舰船目标和海杂波相位图像

Figure 2. TerraSAR-X ship target and sea clutter phase images

下载: 全尺寸图片幻灯片

图 3 TerraSAR-X舰船目标和海杂波相位直方图(以玫瑰图形式展示，红线表示平均方向)

Figure 3. TerraSAR-X ship target and sea clutter phase histograms (Presented in rose charts. Red line indicates mean direction)

下载: 全尺寸图片幻灯片

3.2 Von Mises分布

与一元正态分布相对应，Von Mises分布^[66]通常被称为圆上的正态分布(循环正态分布)，是一种描述方向数据的重要模型。它是N维球面上Von Mises-Fisher分布的一种特例。作为正态分布在圆上的模拟，其可以描述SAR图像中相位的分布。循环随机变量的概率密度函数表达式为

${p_{vM}}\left( {x|\mu ,\kappa } \right) = \frac{{{{\rm e}^{\kappa \cos \left( {x - \mu } \right)}}}}{{2\pi {I_0}\left( \kappa \right)}}$

(12)

其中， ${I_0}\left( \kappa \right)$ 是修正的0阶Bessel函数， ${I_0}\left( \kappa \right) = \dfrac{1}{\pi }\displaystyle\int_0^{\pi} {{{\rm e}^{\kappa \cos \left( x \right)}}} {\rm d}x$ 。参数 $\mu$ 和 $\dfrac{1}{\kappa }$ 是正态分布中均值和方差的模拟量。参数 $\mu$ 是循环平均值，注意其跟传统的均值意义不同，它是位置的度量。 $\kappa$ 表征相位聚集的程度(所以 ${1}/{\kappa }$ 是正态分布方差的模拟量)。如果 $\kappa$ 为0，分布将变为均匀分布；如果 $\kappa$ 增大，分布将趋于以 $\mu$ 为均值 ${1}/{\kappa }$ 为方差的正态分布。图4展示了不同 $\kappa$ 值下的von Mises概率分布情况。

图 4 不同参数下von Mises分布示例

Figure 4. PDFs for von Mises distribution at different parameters

下载: 全尺寸图片幻灯片

Von Mises分布参数估计可以采用最大似然估计方法。其中参数 $\mu$ 是循环均值，其最大似然估计结果定义为^[64]

$\hat \mu = {\rm{Arg}}\left\{ {\frac{1}{N}\sum\limits_{n = 1}^N {\exp \left( {{\rm j}{x_n}} \right)} } \right\}$

(13)

其中 $N$ 是相位观测值的总数， ${x_n}$ 表示一个特定的相位观测值。

为了估计参数 $\kappa$ ，首先需要估计平均组合长度的平方统计量 ${\bar R^2}$ ^[64]

${\bar R^2} = {\bar{ x}}{\bar{ x}}^* = {\left( {\frac{1}{N}\sum\limits_{n = 1}^N {\cos {x_n}} } \right)^2} + {\left( {\frac{1}{N}\sum\limits_{n = 1}^N {\sin {x_n}} } \right)^2}$

(14)

${\bar R^2}$ 的期望值为^[64]

$\left\langle {{{\bar R}^2}} \right\rangle = \frac{1}{N} + \frac{{N - 1}}{N}\frac{{{I_1}{{\left( \kappa \right)}^2}}}{{{I_0}{{\left( \kappa \right)}^2}}}$

(15)

且式(16)统计量为 $\dfrac{{{I_1}{{\left( \kappa \right)}^2}}}{{{I_0}{{\left( \kappa \right)}^2}}}$ 的无偏估计

$\bar R_e^2 = \frac{N}{{N - 1}}\left( {{{\bar R}^2} - \frac{1}{N}} \right)$

(16)

因此，当 $N$ 足够大时，求解下述等式可以得到参数 $\kappa$ 的最大似然估计结果^[64]

$\bar R = \frac{{{I_1}\left( \kappa \right)}}{{{I_0}\left( \kappa \right)}}$

(17)

其中， ${I_m}\left( \kappa \right)$ 为定义如下的Bessel函数

${I_m}\left( \kappa \right) = \frac{1}{\pi }\int_0^{\pi} {\exp \left( {k\cos x} \right)} \cos \left( {mx} \right){\rm d}x$

(18)

3.3 邻域相位方向差

单通道SAR图像中的相位图像通常近似服从均匀分布，直接对其进行建模看起来是没有意义的。相位图像必须经过一定的特殊处理才能更好地提取其中的信息。考虑极化或干涉SAR图像应用中相位信息的使用，必须具有两幅以上同一地区的图像才能有效提取一些信息，即其中必须有1幅图像作为参考图像使用，从而形成1幅主图像和1幅从图像。

考虑使用相位图像的邻域平均相位图像来生成一个参考图像，再通过处理原始相位图像和邻域平均相位图像的相对关系来获取相位信息，本文将这种特征称为邻域相位方向差(Neighborhood Phase Direction Difference, NPDD)。具体步骤为：

(1) 由复值SAR图像 $f\left( {u,v} \right)$ 得到其对应的相位图像 $\theta \left( {x,y} \right)$

$\theta \left( {x,y} \right) = {\rm{Arg}}\left( {f\left( {u,v} \right)} \right)$

(19)

(2) 计算相位图像在不同大小邻域 ${W_r}$ 内的循环平均值图像 ${\hat \theta _r}\left( {x,y} \right)$

${\hat \theta _r}\left( {x,y} \right) = {\rm{Arg}}\left( {\frac{1}{N}\sum\limits_{{W_r}} {\exp \left( {{\rm j}\theta \left( {x,y} \right)} \right)} } \right)$

(20)

其中 $r$ 表示邻域半径，规定若 $r > 1$ , ${W_r}$ 则表示不包括中心点的邻域；若 $r = 1$ , ${W_r}$ 表示中心点本身； $N$ 是邻域内相位观测值的数量。

(3) 最终由不同平均相位图像之差得到相对相位图像

${\rm{NPD}}{{\rm{D}}_{{{st}}}}\left( {x,y} \right) = {\rm{Arg}}\left( {\frac{{\exp \left( {{\rm j}{{\hat \theta }_s}\left( {x,y} \right)} \right)}}{{\exp \left( {{\rm j}{{\hat \theta }_t}\left( {x,y} \right)} \right)}}} \right)$

(21)

其中 ${\hat \theta _s}\left( {x,y} \right)$ 和 ${\hat \theta _t}\left( {x,y} \right)$ 分别表示邻域半径为 $s$ 和 $t$ 的平均相位图像， ${{\rm NPDD}_{st}}\left( {x,y} \right)$ 表示 ${\hat \theta _s}\left( {x,y} \right)$ 和 ${\hat \theta _t}\left( {x,y} \right)$ 的差值图像。注意求解复值商的幅角，相当于求解分子分母之间的相位夹角，因此这里将 ${\rm{NPD}}{{\rm{D}}_{{{st}}}} \left( {x,y} \right)$ 称为 ${\hat \theta _s}\left( {x,y} \right)$ 和 ${\hat \theta _t}\left( {x,y} \right)$ 的差值图像。对于平滑过渡的同质表面，如海洋表面的局部区域，邻域相位方向差理论上为0。而对于相位发生剧烈跳变之处，邻域相位方向差是不为0的。邻域相位方向差可以在某种程度上反映物体表面几何结构的变化情况，若表面几何变化较平缓，则邻域相位方向差较小，否则邻域相位方向差将较大。这就是邻域相位方向差工作的基本原理。在计算邻域相位方向差的过程中，由于 $s$ 和 $t$ 代表邻域范围的大小，因此当SAR图像分辨率较高时，二者取值可稍大，分辨率较低时，二者取值应该比较小，确保比较的邻域范围要足够小，以减少范围增大、散射点增多带来的随机性。一般实际中取 $s = 3$ , $t = 1$ 即可。

El-Darymli等人^[67-69]提出了一种名为后向散射相对相位图像(Backscatter Relative Phase Image, BRPI)的特征来对相位进行描述，其通过原始SAR图像相位图像和经过核卷积之后的相位图像之差得到。但其没有对该过程的具体含义进行很好地解释。El-Darymli等人提到使用不同形式的核(通常中间为0，周围为1)来进行卷积运算，其实质则是挑选不同的邻域相位进行计算。本节的步骤(2)明确了其卷积操作实质就是求取这些相位的循环平均值。根据循环平均值的概念可知，这实际是在求取这些邻域相位的平均方向。由此可知，本节的步骤(3)的意义是求取不同邻域相位的平均方向之间的夹角。最终通过组合不同的邻域，可以求得多个不同的夹角，以此来希望获得能反应目标信息的一些相位信息。后向散射相对相位图像是邻域相位方向差图像当 $t = 1$ 时的特例，邻域相位方向差图像则是包括了所有邻域平均相位图像之差的一般形式。本文提出的邻域相位方向差已经包括了所有后向散射相对相位图像的有效情形。

以取 $s = 3$ , $t = 1$ 为例，获得3.1节中舰船目标和海杂波对应的邻域相位方向差图像 ${\rm{NPD}}{{\rm{D}}_{31}}$ 。如图5所示， ${\rm{NPD}}{{\rm{D}}_{31}}$ 图像主体以0为主，说明 ${\hat \theta _3}$ 和 ${\hat \theta _1}$ 的方向角对大部分背景而言是一致的。而舰船目标强散射点引起的相位图像中的条纹得到增强，在 ${\rm{NPD}}{{\rm{D}}_{31}}$ 图像中变得更加明显。图6是对应的邻域相位方向差图像的直方图。由直方图结果可以看出，舰船目标和海杂波的邻域相位方向差不再服从均匀分布，相位平均值均指向0附近，而相位长度存在一定区别。

图 5 TerraSAR-X舰船目标和海杂波邻域相位方向差图像

Figure 5. TerraSAR-X ship target and sea clutter NPDD images

下载: 全尺寸图片幻灯片

图 6 TerraSAR-X舰船目标和海杂波邻域相位方向差直方图(以玫瑰图形式展示，红线表示平均方向)

Figure 6. TerraSAR-X ship target and sea clutter NPDD histograms (Presented in rose charts. Red line indicates mean direction)

下载: 全尺寸图片幻灯片

表2详细列出了邻域相位方向差的循环统计量结果。对比原始相位图像循环统计量结果可以看出，邻域相位方向差的循环统计量呈现出一定的规律，这也可以直接从邻域相位方向差的图像直观反映得到。具体而言，相比于海杂波，舰船目标的循环方差、标准差更大，同时长度和峰度更小，这说明舰船目标的邻域相位方向差是更分散的，这和邻域相位方向差的工作原理一致。非0邻域相位方向差主要出现在二面角存在的地方，验证了邻域相位方向差的基本工作原理。因为即使邻域范围很小，舰船目标和海面交界处一般会发生相位的剧烈跳变，相反海杂波的相位变化缓慢，在邻域平均的过程中相位变化得到了抵消从而相位差大多数接近0。此外，舰船目标和海杂波的均值和偏度都近似为0，这说明二者平均方向是一致的且关于平均方向都是对称的。

表 2 TerraSAR-X舰船目标和海杂波的邻域相位方向差图像循环统计量结果

Table 2. Circular statistical results of TerraSAR-X ship target and sea clutter NPDD images

循环统计量	均值	方差	标准差	平均长度	偏度	峰度
(a) 舰船目标1	0	0.27	0.79	0.73	0	0.46
(b) 海杂波1	0	0.21	0.69	0.79	0	0.53
(c) 舰船目标2	0	0.37	0.96	0.63	0	0.35
(d) 海杂波2	0	0.21	0.69	0.79	0.01	0.53

下载: 导出CSV

| 显示表格

4. 单通道SAR图像复值统计建模方法

本节主要采用复广义高斯分布对单通道复值SAR图像进行建模研究。复广义高斯分布是一类直接描述复值数据的重要分布，其形状参数是描述非高斯性的一种重要复值信息。本节将首先对广义高斯分布进行介绍，然后对其复值形式即复广义高斯分布进行介绍，并介绍其经典的参数估计方法^[30,31]，最后介绍文献[70]中提出的形状参数快速估计方法。

4.1 广义高斯分布简介

广义高斯分布是一组对称分布族，包括了从超高斯(sup-Gaussian)分布到亚高斯(sub-Gaussian)分布等很多种情形，在极限情况下还包括了所有连续均匀分布和实线的截断^[30,31]。这使得广义高斯分布具有非常广泛的应用。

广义高斯分布通过改变高斯分布的衰减速率而得到，其概率密度函数如下^[30,31]

$\begin{split} {p_{{\rm{CGD}}}}\left( {x|\mu ,\sigma ,\beta } \right) =\,& \frac{\beta }{{2\sigma \varGamma \left( {{1 / \beta }} \right)}}\exp \left[ { - {{\left( {\frac{{\left| {x - \mu } \right|}}{\sigma }} \right)}^\beta }} \right],\\ & x \in {\mathbb{R}}\\[-10pt] \end{split}$

(22)

其中 $\mu$ 为位置参数或者说是均值， $\sigma$ 为标准差， $\beta$ 为形状参数， $\varGamma \left( \cdot \right)$ 为伽马函数。拉普拉斯分布是广义高斯分布当 $\beta = 1$ 时的特例，高斯分布是广义高斯分布当 $\beta = 2$ 时的特例。当 $\beta < 2$ 时，广义高斯分布为超高斯分布；当 $\beta > 2$ 时，广义高斯分布为亚高斯分布；而当 $\beta \to 0$ 和 $\beta \to \infty$ 时，广义高斯分布的概率密度函数分别趋近于脉冲函数和均匀分布。图7展示了0均值下不同 $\beta$ 和 $\sigma$ 下的广义高斯分布概率。从中可以明显看出，形状参数 $\beta$ 越小，概率分布函数顶峰越尖锐；标准差 $\sigma$ 越小，概率分布函数形状越集中于均值 $\mu$ 。

图 7 不同参数下的广义高斯分布示例

Figure 7. PDFs for CGGD at different parameters

下载: 全尺寸图片幻灯片

4.2 复广义高斯分布及其生成

广义高斯分布的复值形式，即复广义高斯分布，更适合于用来描述SAR图像。这是因为SAR聚焦后的原始数据是复值形式的。一个复值随机变量可以直接表示为 $Z = I + {\rm j}Q$ ，其中 $I$ 和 $Q$ 分别代表其实部和虚部；也可以表示为双实变量的矢量形式 ${{{z}}_b} = {\left[ {I,Q} \right]^{\rm{T}}}$ ，同时假设 ${\rm{E}}\left[ {{{{z}}_b}} \right] = {\bf{0}}$ ，则其增广向量形式可以表示为 ${{z}} = {\left[ {z,z^*} \right]^{\rm{T}}}$ , ${{z}} = \left[ {\begin{array}{*{20}{c}} 1& {\rm j} \\ 1& { - {\rm j}} \end{array}} \right]{{{z}}_b}$ ，因此双变量协方差矩阵为^[30,31]

${{{C}}_b} = {\rm{E}}\left[ {{{{z}}_b}{{z}}_b^{\rm{T}}} \right] = \left[ {\begin{array}{*{20}{c}} {\sigma _I^2}&\rho \\ \rho &{\sigma _Q^2} \end{array}} \right]$

(23)

增广协方差矩阵为

$\begin{split} {{C}} \,& = {\rm{E}}\left[ {{{z}}{{{z}}^{\rm{H}}}} \right] \\ & = \left[ {\begin{array}{*{20}{c}} {\sigma _I^2 + \sigma _Q^2}&{\left( {\sigma _I^2 - \sigma _Q^2} \right) + {\rm j}2\rho } \\ {\left( {\sigma _I^2 - \sigma _Q^2} \right) - {\rm j}2\rho }&{\sigma _I^2 + \sigma _Q^2} \end{array}} \right] \end{split}$

(24)

采用增广协方差矩阵形式，复广义高斯分布概率密度函数如下^[30,31]

$\begin{split} {p_{\rm CGGD}}\left( {{{Z}}|\beta ,{{C}}} \right) =\,& \frac{\beta }{{\pi \sqrt {\left| {{C}} \right|} }}\frac{{\varGamma \left( {{2 / \beta }} \right)}}{{\varGamma {{\left( {{1 / \beta }} \right)}^2}}}\; \\ &\cdot \exp \left\{\!\! { - {{\left[ {\frac{{\varGamma \left( {{2 / \beta }} \right)}}{{2\varGamma \left( {{1 / \beta }} \right)}}\left( {{{{Z}}^{\rm{H}}}{{{C}}^{ - 1}}{{Z}}} \right)} \right]}^\beta }} \! \right\} \end{split}$

(25)

其中 $\beta$ 代表形状参数， ${{C}}$ 代表增广协方差矩阵， ${\rm{H}}$ 代表共轭转置。令 $c\left( \beta \right) = \dfrac{{\varGamma \left( {{2 / \beta }} \right)}}{{\varGamma \left( {{1 / \beta }} \right)}}$ ，可将式(25)转换为形式上更类似于广义高斯分布的形式

$\begin{split} {p_{\rm CGGD}}\left( {{{Z}}|\beta ,{{C}}} \right) =\,& \frac{{\beta c\left( \beta \right)}}{{\pi \sqrt {\left| {{C}} \right|} \varGamma \left( {{1 / \beta }} \right)}}\; \\ & \cdot \exp \left\{ { - {{\left[ {\dfrac{1}{2}c\left( \beta \right)\left( {{{{Z}}^{\rm{H}}}{{{C}}^{ - 1}}{{Z}}} \right)} \right]}^\beta }} \right\} \end{split}$

(26)

不同于广义高斯分布，当 $\beta = 1$ 时，复广义高斯分布为复高斯分布；当 $\beta < 1$ 时，复广义高斯分布描述超高斯分布；当 $\beta > 1$ 时，复广义高斯分布描述亚高斯分布。注意当 ${{C}}$ 为对角阵时，即 $\left( {\sigma _I^2 - \sigma _Q^2} \right) \pm {\rm j}2\rho = 0$ 时，该复值随机变量是圆性的。

通常自然地物的实部和虚部分布是相同的，都接近一个高斯分布；而人造目标的实部和虚部分布可能存在一定的差异，且其实部和虚部都更接近于广义高斯分布的尖峰形态。由于形状参数可以反映概率分布函数顶峰的尖锐程度，表征数据的非高斯性，因此通常用形状参数可以区分自然地物和人造目标在形态上的差异。图8展示了舰船目标和海杂波的实部和虚部图像，而图9展示了舰船目标和海杂波直方图分布的异同。

图 8 舰船目标和海杂波的实部和虚部图像

Figure 8. Real and imaginary parts of ship target and sea clutter images

下载: 全尺寸图片幻灯片

图 9 舰船目标和海杂波实部和虚部的直方图

Figure 9. Real and imaginary histograms of ship target and sea clutter images

下载: 全尺寸图片幻灯片

复广义高斯分布样本的生成类似于其概率密度函数的推导过程。首先产生两个标准化的随机双变量，然后将其代入增广形式，最后经过一个变换产生希望的协方差。给定形状参数 $\beta$ ，增广协方差矩阵 ${{C}}$ ，以下过程生成 $N$ 个独立的复变量^[30,31]：

(1) 产生 $n = 1,2,\;···\;,N$ 个复值样本

$\qquad\ z\left( n \right) = {\rm{gamrnd}} {\left( {{1 / {\beta ,1}}} \right)^{{1 / {2\beta }}}}\exp \left( {{\rm j}2\pi {\rm{rand}} } \right)$

(27)

(2) 标准化复方差

$w = {z / {{\eta _c}\left( c \right)}},\;{\eta _c}\left( c \right){\rm{ = }}\frac{{\varGamma \left( {{2 / \beta }} \right)}}{{\varGamma \left( {{1 / \beta }} \right)}}$

(28)

(3) 形成增广向量

${{w}}{\rm{ = }}{\left[ {w,\;{\rm{conj}}\left( w \right)} \right]^{\rm{T}}}$

(29)

(4) 计算 ${{C}}$ 的矩阵平方根得到变换矩阵

${{T}} = {\rm{sqrtm}} \left( {{C}} \right)$

(30)

(5) 进行变换

${{v}} = {{Tw}}$

(31)

其中gamrnd, sqrtm, conj以及rand均为MATLAB里面提供的函数。

4.3 复广义高斯分布参数估计方法

4.3.1 基于最大似然估计法的复广义高斯分布参数估计

Novey等人^[30,31]基于最大似然法提出了一种复广义高斯分布的参数估计方法。其是目前文献中使用最为广泛的复广义高斯分布估计方法。据谷歌学术统计，其被引超过100余次，是当前主流的复广义高斯分布参数估计方法。但是，在实际应用中Novey的方法计算效率非常低，且在一些情况下估计精度较低。其具体步骤为^[30,31]

(1) 首先使用矩估计方法估计初始形状参数 $\beta$

$\begin{split} & \,\beta = \arg {\min _{\beta} }\;\left| {{{f}} \left( \beta \right)} \right|\; \\ & \, {{f}} \left( \beta \right) = \left( {\frac{{{\rm{E}}\left[ {{x^4}} \right]}}{{{\rm{E}}{{\left[ {{x^2}} \right]}^2}}}} \right) \!+\! \left( {\frac{{{\rm{E}}\left[ {{y^4}} \right]}}{{{\rm{E}}{{\left[ {{y^2}} \right]}^2}}}} \right) \!-\! \left( {\frac{{3\varGamma \left( {{1 / \beta }} \right)\varGamma \left( {{3 / \beta }} \right)}}{{\varGamma {{\left( {{2 / \beta }} \right)}^2}}}} \right) \\ \end{split}$

(32)

同时利用样本协方差估计初始增广矩阵 ${{C}}$ 。

(2) 然后使用Newton-Raphson方法更新 $\beta$

${\beta ^{n + 1}} = {\beta ^n} - \left( {{{{{\partial L} / {\partial \beta }}} \Bigr/ {{{{\partial ^2}L} / {{\partial ^2}\beta }}}}} \right)$

(33)

其中 $L$ 表示概率密度函数的对数，同时如果 $\beta < 1$ 则使用定点方法更新 ${{C}}$

${{C}} = \left( {{{2\beta {{\left( {\frac{1}{2}c\left( \beta \right)} \right)}^\beta }} \Biggr/ N}} \right){\sum {\left( {{{{Z}}^{\rm{H}}}{{{C}}^{ - 1}}{{Z}}} \right)} ^{\beta - 1}}{{ZZ}}{^{\rm{H}}}$

(34)

4.3.2 基于高阶矩的复广义高斯分布形状参数快速估计

形状参数是复广义高斯分布的一个重要参数，根据Novey的方法对其进行估计计算效率很低，因此有必要提出新的估计方法。文献[70]基于高阶矩和形状参数的关系提出了一种快速的复广义高斯分布形状参数的估计方法。其推理过程如下^[70]：

根据复广义高斯分布和对称Kotz型分布的关系可知非0矩可以表示为^[29,71]

${M_{n,n}} = {\rm{E}}\left[ {{U^n}} \right] = {\alpha ^{{n / \beta }}}\frac{{\varGamma \left( {{{\left( {n + 1} \right)} / \beta }} \right)}}{{\varGamma \left( {{1 / \beta }} \right)}}$

(35)

其中 ${\rm{E}}\left[ U \right] = {\rm{E}}\left[ {{{\left| Z \right|}^2}} \right]$ 。因此

$\frac{{{\rm{E}}\left[ {{U^2}} \right]}}{{{{\left( {{\rm{E}}\left[ U \right]} \right)}^2}}} = \frac{{\varGamma \left( {{1 / \beta }} \right)\varGamma \left( {{3 / \beta }} \right)}}{{{{\left( {\varGamma \left( {{2 / \beta }} \right)} \right)}^2}}}$

(36)

根据式(36)即可确定高阶矩和形状参数的关系。

同时注意到在文献中应用最为广泛的一种复信号峰度高阶矩(Complex Signal Kurtosis, CSK)，其可以表示为^[29,72,73]

${\rm{CSK}} = {\mu _{2,2}} - 2 - {\left| {{\mu _{2,0}}} \right|^2}$

(37)

其中 ${\mu _{l,m}}$ 为标准矩^[29]。因此

${\rm{CSK}} + {\left| {{\mu _{2,0}}} \right|^2} + 2 = \frac{{{\rm{E}}\left[ {{U^2}} \right]}}{{{{\left( {{\rm{E}}\left[ U \right]} \right)}^2}}}$

(38)

考虑到复值SAR图像均值通常为0，同时假设复随机变量 $Z$ 是圆性的并且存在有限的4阶矩，则有 ${\mu _{2,0}} = 0$ ，且

${\rm{CSK}} = \frac{{\varGamma \left( {{1 / \beta }} \right)\varGamma \left( {{3 / \beta }} \right)}}{{{{\left( {\varGamma \left( {{2 / \beta }} \right)} \right)}^2}}} - 2$

(39)

由此结果可知，复信号峰度可以由一个关于 $\beta$ 的显示表达式表示。图10展示了复信号峰度和形状参数的关系。注意由于负的复信号峰度值无法用对数形式表示，复信号峰度被分为正和负的两部分，也就是说纵轴表示的是复信号峰度的绝对值。由图可以看出，复信号峰度随形状参数增大而单调递减。当 $\beta = 1$ 时， ${\rm{CSK}} = 0$ ，此时样本服从复高斯分布；当 $\beta < 1$ 时， ${\rm{CSK}} > 0$ ，此时样本服从复超高斯分布；当 $\beta > 1$ 时， ${\rm{CSK}} < 0$ ，此时样本服从复亚高斯分布。虽然 $\beta$ 很难通过复信号峰度的显式表达式表示，但是可以通过数值计算方法通过复信号峰度求解 $\beta$ 。文献[70]采用查找表(LookUp Tables, LUTs)的方法来解决这个问题。图11展示了具体流程图，总结起来其可以分为3步^[70]：

图 10 复信号峰度和形状参数的关系

Figure 10. Relationship between CSK and the shape parameter

下载: 全尺寸图片幻灯片

图 11 基于复信号峰度的复广义高斯分布形状参数估计方法流程

Figure 11. Flowchart of shape parameter estimation of CGGD based on CSK

下载: 全尺寸图片幻灯片

(1) 生成查找表。根据其关系生成2种不同精度的查找表，一个低精度查找表和一个高精度查找表。 $\beta$ 采用0～10之间的数，因为绝大多数的SAR图像分布范围都在该区间以内。 $\beta$ 的间隔使用两种不同的精度，分别为1e-1和1e-3。然后根据关系可以得到低精度查找表(大约1 kb)和高精度查找表(大约100 kb)。这样应该有助于满足不同使用者的需求。

(2) 估计样本的复信号峰度。样本的复信号峰度可以通过式(37)得到。这里将估计得到的复信号峰度值记为 $\widehat {{\rm{CSK}}}$ 。最后 $\widehat {{\rm{CSK}}}$ 将被用于查找表。

(3) 使用查找表得到形状参数。在查找表中搜索 $\widehat {{\rm{CSK}}}$ 的最近邻值。然后就可以知道最近邻复信号峰度值在查找表中的索引，根据这个索引就可以得到对应的形状参数。本文将估计得到的形状参数值记为 $\hat \beta$ 。本文在查找过程中将使用高精度查找表。

注意查找表只要生成过1次，即可反复使用。因此在实际操作过程中，估计形状参数只需步骤2和步骤3。该方法中的查找表和相应的MATLAB代码可以通过网站^[74]自由获取。

使用4.2节中生成数据的方法来仿真复广义高斯分布样本。在样本仿真生成的过程中， ${{C}}$ 每回都是随机产生的，而 $\beta$ 从0.1变化到4，其中间隔为0.1。每一个特定大小的 $\beta$ 将会被测试100次，因此对于每个 $\beta$ 本节将生成100组不同的复广义高斯分布样本。最后Novey的方法以及文献[70]中提出的基于复信号峰度的方法将被用来估计这些样本的形状参数。同时设置了不同的样本数量( ${\rm{size}}$ )，分别为500, 1500, 5000, 50000。

估计精度采用均方误差(Mean Squared Error, MSE)表示，结果如图12所示。具体地，当样本数量较小时( ${\rm{size}} \le 500$ )，对于较小的 $\beta$ ( $\beta \le 2.5$ 大致)，提出的方法和Novey的方法效果近似。此时，Novey的方法和提出的方法都倾向于得到比真值更大的估计结果，但是总体的偏差是比较小的。当样本的数量不断增加，提出方法的均方误差显著下降并且具有比Novey方法更好的性能。当样本数量 ${\rm{size}} > 1500$ ，对于较大 $\beta$ 值的估计性能显著提升。当样本数量增加到50000时，提出方法取得了非常好的估计结果，均方误差接近于0。这同时证明了提出方法相合性更好，也就是说只要样本数量足够多，提出方法就能取得较好的估计结果。对于Novey的方法而言，当样本数量增大时，其性能同样得到了提升。但是，Novey方法的均方误差是波动的，近似于一个锯齿状但有许多零点。有趣的是当样本数量增加时，这种现象会越加明显。这些0点大部分来自于Novey方法设置的初始值，其余一些来源于最大似然估计方法。当样本数量增大时，Novey的方法倾向于把形状参数估计为临近设置的初始形状参数，因此其结果会呈现出一个锯齿状。

图 12 Novey的方法和基于复信号峰度的方法的均方误差结果对比

Figure 12. MSE results of Novey’s method and our method

下载: 全尺寸图片幻灯片

图13展示了耗时对比。由图可以看出，当样本数量为500时，Novey的方法耗时约8.3e-2 s，提出的方法耗时约1.5e-3 s，效率提升了55倍；当样本数量增加到50000时，Novey的方法耗时约6.3 s，而提出的方法耗时约4.5e-3 s，效率提升了1400倍。随着样本的增多，提出方法在时间消耗方面的优势变得更加明显。这是因为Novey的方法由于增多的样本耗费了更多的时间，而提出的方法仅在计算复信号峰度时增加了时间并且计算复信号峰度这个过程始终是非常简洁高效的。

图 13 Novey的方法和基于复信号峰度的方法的单次估计平均时间消耗对比

Figure 13. Average time consumption comparison for a single test of Novey’s method and our method

下载: 全尺寸图片幻灯片

总结而言，基于大量仿真实验数据的结果表明基于复信号峰度的方法在估计精度方面比Novey的方法具有一定的优势；在耗时方面，基于复信号峰度的方法比Novey的方法则具有显著的优势。

5. 基于复值统计信息的单通道SAR舰船目标检测部分研究结果

作者研究小组在基于复值统计信息的单通道SAR图像舰船目标检测和分割等方面进行了初步研究^{[21,22,74,75]}，结果表明采用单通道复值SAR图像进行舰船目标检测可以更好地利用复值信息，可以有效缓解传统基于幅度图像进行检测产生的一些问题。部分最新研究成果和结论如下。

5.1 基于复信号峰度的复杂环境背景舰船目标检测

针对复杂环境背景下的舰船目标检测问题，文献[21]根据复信号峰度对海杂波和舰船目标敏感的优点，设计了一种基于复信号峰度的舰船目标检测方法。其检测基本原理是海杂波的复信号峰度在不同分辨率、波段、极化、入射角成像参数条件下近似恒定，维持在0附近一个较窄的区间内^[21]，如图14(a)所示，这和其幅度或强度在不同成像条件下的巨大变化有明显区别，如图14(b)所示。而舰船目标呈现出和海杂波完全不同的复信号峰度特性。舰船目标的复信号峰度在很大一个区间内剧烈变化并且远大于0，如图14(c)所示(此处舰船目标样本来源于经OpenSARShipFilter工具^[76]过滤后的OpenSARShip数据集^[77,78])。因此通过设置一个全局阈值即可以有效检测舰船目标，这种性质非常有助于海洋背景中舰船目标的快速筛选，例如可以较好地克服宽幅SAR图像入射角变化带来的辐射不均衡问题。

图 14 不同成像条件下海杂波和舰船目标的复信号峰度对比

Figure 14. CSK plots of sea clutter of typical sea clutter and ship targets from different acquisitions

下载: 全尺寸图片幻灯片

此外，欧盟JRC的Santamaria等人^[79]通过SUMO软件对地中海地区超过1万幅的Sentinel-1图像进行批量舰船目标检测时，发现由于受到射频干扰的影响，结果中存在大量虚警，不得不舍弃其中不少数据或者特定地区的数据。Sentinel-1中的射频干扰如图15所示。文献[22]对舰船目标检测中的射频干扰鉴别问题进行了研究。通过分析舰船目标和射频干扰产生机理的不同，文献[22]发现了射频干扰复值信息的两个特点，即射频干扰呈现非常弱的非圆性并且主要呈现亚高斯性，而舰船目标呈现强非圆性和超高斯性；其虽然和舰船目标具有近似的强度，但是其非圆性和非高斯性却截然不同。因此通过非圆性和非高斯性等复值特征的判断可以有效区分射频干扰和舰船目标等人造目标^[22,70]。文献[22]中提出的基于非圆性和非高斯性方法的总体精度为97.94%，而基于幅度图像纹理特征的SVM方法精度仅为88.91%，说明了复值信息的有效性。

图 15 Sentinel-1图像中的射频干扰现象(绿圆表示舰船目标)

Figure 15. Ships affected by RFIs in Sentinel-1 images (The green circles represent ships)

下载: 全尺寸图片幻灯片

文献[21]中提出的基于复信号峰度的方法可以较好地结合非圆性和非高斯性，通过对多种典型复杂环境场景进行实验，结果表明基于复信号峰度的检测方法可以有效减少射频干扰引起的虚警，同时减轻方位模糊造成的影响。此外，提出的方法能够有效检测密集目标状态下的舰船目标，要优于传统基于幅度图像的CFAR检测方法。更多细节可参考文献[21,22]。

5.2 结合幅度信息和复信号峰度的鲁棒舰船目标迭代分割

针对单通道复值SAR图像舰船目标分割问题，文献[21,75]中提出了一种结合幅度信息和复信号峰度的舰船目标复值图像鲁棒迭代分割方法。

其分割原理为：设有一个同时包含舰船目标和海杂波的舰船切片，首先假设海杂波的实部和虚部都服从高斯分布，同时只考虑海杂波中的实部或虚部及其对应的实信号峰度。不妨以正态分布为例，如图16(a)所示，可以看出，对于高斯分布而言，较小的点集中在0附近而较大的点是分散的。由于实信号峰度值的大小更多地取决于离群值而不是峰值的^[80]，因此，当幅度阈值 $t$ 很小时(如0.5倍幅度均值)，包含点的分布的离群倾向弱于高斯分布，此时实信号峰度小于0；当阈值 $t$ 不断地增大时，包含的离群点越来越多，这个分布的离群倾向越来越强，此时实信号峰度会不断增大，如图16(b)所示；最后，当阈值 $t$ 增大到理想的分割阈值 $t^*$ 时，此时这个分布是一个包括了所有点的高斯分布，因此实信号峰度等于0。总结起来，当阈值 $t$ 由一个很小的值变化到最终的阈值 $t^*$ 的过程中，实信号峰度的值将由–1.2 (其实质为均匀分布的实信号峰度，注意高斯分布的顶部局部区域可近似为一个均匀分布)变化到0(即高斯分布的实信号峰度)。由于舰船目标的像素值通常远大于海杂波或者说相较于海杂波离群倾向更明显，因此迭代过程在包含舰船目标像素之前就会终止；否则，实信号峰度的值将会大于0。

图 16 高斯分布迭代分割过程的示例

Figure 16. Illustration of the iteration process for a Gaussian distribution

下载: 全尺寸图片幻灯片

提出的迭代分割方法实际上可以归纳为一种从混合信号(高斯信号和超高斯信号)中提取高斯信号的方法。在提取高斯信号的过程，同样得到了超高斯信号，整个过程完成了舰船目标和海杂波的分离。因此，提出的方法也可以单独作为一种从SAR图像自然背景中分割人造目标的方法。不同于常规的SAR图像阈值分割方法，该方法能够有效利用SAR图像中海杂波的复值信息，从而确定合理的分割阈值，直接对复值SAR图像进行分割。更多细节可参考文献[21,75]。

图17展示了对部分Sentinel-1 SM和高分三号UFS模式舰船目标复值切片进行迭代分割的结果，并和经典的Otsu图像分割方法进行了对比。Otsu图像分割方法源于对普通光学图像的应用，在SAR图像分割中其应用对象通常是16位幅度图像拉伸映射后的8位图像。由结果可以看出，迭代分割方法的结果更加稳定和可靠；而Otsu方法容易过分割，特别是对高分三号UFS模式这种高分辨率数据中的一些杂波动态范围比较大的切片。总的来说，Otsu方法的效果非常依赖于拉伸映射后的8位图像质量，而提出的迭代分割方法可以有效利用SAR图像目标和背景复值统计特性的区别，结果更加地鲁棒。

图 17 舰船目标切片复信号峰度迭代分割和Otsu分割对比结果

Figure 17. Comparison of Otsu and CSK iteration segmentation results

下载: 全尺寸图片幻灯片

6. 结束语

本文面向舰船目标检测，较为全面系统地综述了单通道复值SAR图像统计建模方法。不同于传统惯例着重于单通道SAR图像幅度统计建模方法的研究，本文重点对单通道SAR图像相位和复值统计建模方法进行了研究，同时有针对性地介绍了作者研究小组的一些最新研究成果。随着SAR图像分辨率的提升以及对舰船目标检测性能要求的提高，本文研究将对理解单通道SAR图像的复值本质特性及其应用具有重要意义，同时对促进相关研究者的工作具有积极作用。但是关于单通道SAR图像复值信息挖掘及其应用需要进一步深入研究的还有很多，主要包括以下3个方面：

(1) 超高分辨率单通道SAR图像复值信号处理。目前，国内外已有多部机载SAR系统具有亚米级分辨率成像能力。例如，中科院毫米波机载系统可以获得0.15×0.15 m的分辨率^[34,35]，国外MEMPHIS系统^[81]可以获得0.06 m的超高分辨率，同时星载SAR系统的分辨率也在不断提高。随着SAR分辨率的提高，完全发展的相干斑的假设条件变得更加地不适用。超高分辨率SAR不仅可以获取目标更丰富的幅度信息，其蕴含的复值信息如非圆性等也将在目标识别领域发挥更重要的作用。发展超高分辨率单通道SAR图像复值信号处理方法和理论具有重要的研究意义，对充分释放超高分辨率SAR的应用潜能具有重要作用，是一个值得研究的发展方向。

(2) 单通道复值信息和多模式SAR图像处理。单通道复值信号处理侧重于对单通道SAR图像本身复值信息的提取，而随着未来SAR技术的进一步发展，方位向多通道的宽幅高分SAR、MIMO-SAR、双/多基SAR、多极化SAR、干涉SAR、TomoSAR、视频SAR等多模式SAR^[3,4]都不再局限于单通道SAR图像。单通道SAR图像复值信号处理可被认为是多模式SAR图像处理的重要基础组成部分，因此通过单通道SAR图像复值信息提取有助于促进多模式SAR图像处理，同时多模式SAR图像处理也有助于促进单通道SAR图像复值信息处理方法在多种应用场景中的突破，具有非常广阔的发展空间。

(3) 单通道复值SAR图像解译和深度学习。以卷积神经网络为代表的深度学习在光学图像目标检测识别中取得了令人瞩目的成绩。但是SAR图像成像机理和传统光学图像存在显著差异，图像特征差异也非常明显，直接套用光学领域的深度学习方法在SAR图像解译领域性能受到极大限制。作为SAR图像区别于光学图像的一种本质特征和根本区别，SAR图像的复值信息将有助于深度学习实现更有效的特征表达和数据驱动，是一个非常重要的突破口。将复值SAR图像解译和深度学习进行有机融合是一个非常有前景的研究方向，有望成为突破SAR视觉智能解译瓶颈的重要途径。

图 1 卷积神经网络结构

Figure 1. The structure of CNN

下载: 全尺寸图片幻灯片

图 2 Dropout示意图

Figure 2. The schematic diagram of Dropout

下载: 全尺寸图片幻灯片

图 3 残差块结构图

Figure 3. The structure of residual block

下载: 全尺寸图片幻灯片

图 4 Inception结构图

Figure 4. The structure of inception

下载: 全尺寸图片幻灯片

图 5 降噪自编码器层结构图

Figure 5. The structure of residual DAE

下载: 全尺寸图片幻灯片

图 6 网络结构

Figure 6. The structure of network

下载: 全尺寸图片幻灯片

图 7 网络功能框图

Figure 7. The block diagram of network function

下载: 全尺寸图片幻灯片

图 8 HRRP数据示意图

Figure 8. Schematic diagram of HRRP

下载: 全尺寸图片幻灯片

图 9 卷积神经网络结构图

Figure 9. The structure of CNN

下载: 全尺寸图片幻灯片

表 1 飞机参数(m)

Table 1. Parameters of planes (m)

机型	机长	机高	机宽
安26	23.80	9.83	29.20
奖状	14.40	4.57	15.90
雅克42	36.38	9.83	34.88

下载: 导出CSV

表 2 训练集加入瑞利噪声结果

Table 2. Recognition results based on Rayleigh noise training set

	信噪比(dB)	–5	–3	0	3	5	10	15	20	测试集未加入噪声
图6改进CNN	平均识别率	91.14	93.98	96.14	97.37	97.71	98.58	99.24	99.53	99.81
	安26识别率	88.83	92.97	96.47	97.80	98.50	98.93	99.13	99.40	99.90
	奖状识别率	92.13	95.43	97.93	99.20	99.53	99.63	99.77	99.90	99.97
	雅克42识别率	92.47	93.53	94.03	95.10	95.10	97.17	98.83	99.30	99.57
图6改进CNN (训练集中未加噪声)	平均识别率	51.51	57.89	63.41	65.14	65.72	66.40	66.58	66.68	99.81
	安26识别率	0	0	0	0	0	0	0	0.03	99.47
	奖状识别率	99.80	99.83	99.77	99.87	100.00	100.00	100.00	100.00	99.97
	雅克42识别率	54.77	73.83	90.47	95.57	97.17	99.20	99.73	100.00	100.00
图9的CNN	平均识别率	85.30	90.43	94.90	96.66	97.30	98.19	98.36	98.44	99.63
	安26识别率	72.13	79.27	86.53	90.57	92.10	94.73	95.20	95.47	99.97
	奖状识别率	83.90	92.17	98.23	99.50	99.83	99.90	99.93	99.93	99.67
	雅克42识别率	99.87	99.87	99.93	99.90	99.97	99.93	99.93	99.93	99.27
图9的CNN (训练集中未加噪声)	平均识别率	33.86	34.50	38.12	51.82	61.32	66.62	66.70	68.20	99.97
	安26识别率	0	0	0	0	0	0	0.10	4.60	99.90
	奖状识别率	1.57	3.50	14.37	55.47	83.97	99.87	100.00	100.00	100.00
	雅克42识别率	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00

下载: 导出CSV

表 3 训练集加入高斯白噪声结果

Table 3. Recognition results based on White Gaussian noise training set

	信噪比(dB)	–5	–3	0	3	5	10	15	20	测试集未加入噪声
图6改进CNN	平均识别率	85.99	91.61	96.04	98.07	99.00	99.68	99.88	99.94	99.94
	安26识别率	86.63	90.20	94.00	96.30	98.07	99.30	99.70	99.83	99.87
	奖状识别率	83.27	92.77	97.00	99.23	99.60	99.90	100.00	100.00	100.00
	雅克42识别率	88.07	91.87	97.13	98.67	99.33	99.83	99.93	100.00	99.97
图6改进CNN (训练集中未加噪声)	平均识别率	55.93	58.89	63.91	69.10	70.94	75.49	82.73	93.01	99.90
	安26识别率	12.00	13.47	15.17	20.37	21.33	29.63	49.43	79.40	99.83
	奖状识别率	82.20	84.57	90.57	95.57	97.63	99.63	99.97	100.00	100.00
	雅克42识别率	73.60	78.63	86.00	91.37	93.87	97.20	98.80	99.63	99.87
图9的CNN	平均识别率	80.61	87.23	93.74	96.59	97.59	98.87	99.08	99.21	99.19
	安26识别率	71.37	78.50	86.77	91.73	93.90	96.87	97.63	97.80	97.73
	奖状识别率	80.67	89.53	97.07	99.20	99.43	99.87	99.77	99.93	99.93
	雅克42识别率	89.80	93.67	97.40	98.83	99.43	99.87	99.83	99.90	99.90
图9的CNN (训练集中未加噪声)	平均识别率	39.92	42.74	47.09	53.90	59.90	73.09	86.69	97.47	99.98
	安26识别率	0.57	0.83	1.77	3.37	5.30	22.20	60.27	92.47	99.97
	奖状识别率	21.73	28.90	40.33	58.83	74.63	97.10	99.83	99.97	100.0
	雅克42识别率	97.47	98.50	99.17	99.50	99.77	99.97	99.97	99.97	99.97

下载: 导出CSV

表 4 不同噪声类型识别结果

Table 4. Recognition results based on different noise types

	信噪比(dB)	–5	–3	0	3	5	10	15	20	测试集未加入噪声
方案1	平均识别率	76.73	82.21	88.31	92.12	93.41	96.23	97.97	99.06	99.79
	安26识别率	76.47	80.53	83.90	85.90	86.30	90.90	94.53	97.40	99.57
	奖状识别率	73.37	78.43	87.10	92.67	95.37	98.23	99.47	99.80	99.80
	雅克42识别率	80.37	87.67	93.93	97.80	98.57	99.57	99.90	99.97	100.00
方案2	平均识别率	81.68	87.18	92.53	94.93	95.52	96.70	96.94	97.08	97.18
	安26识别率	63.63	72.97	82.20	87.23	88.50	91.03	91.57	91.87	92.13
	奖状识别率	81.60	88.77	95.50	97.77	98.30	99.27	99.50	99.60	99.63
	雅克42识别率	99.80	99.80	99.90	99.80	99.77	99.80	99.77	99.77	99.77

下载: 导出CSV

表 5 删除结构加入噪声结果

Table 5. Recognition results based on deleted structure

	信噪比(dB)	–5	–3	0	3	5	10	15	20	测试集未加入噪声
图6改进CNN	平均识别率	79.98	86.83	93.11	96.18	97.11	98.59	98.72	98.99	99.03
	安26识别率	73.50	82.03	89.60	94.33	94.97	97.10	97.30	97.70	97.80
	奖状识别率	81.90	89.50	96.63	98.77	99.50	99.77	99.90	99.90	99.87
	雅克42识别率	84.53	88.97	93.10	95.43	96.87	98.90	98.97	99.37	99.43
删除第1个残差块	平均识别率	79.59	84.82	90.76	93.66	94.36	95.68	95.88	96.09	96.04
	安26识别率	68.97	74.70	80.67	84.67	85.97	88.23	88.63	89.00	89.07
	奖状识别率	82.50	88.43	96.80	98.83	99.50	99.73	99.83	99.83	99.83
	雅克42识别率	87.30	91.33	94.80	97.47	97.60	99.07	99.17	99.43	99.23
删除第2个残差块	平均识别率	79.06	85.53	91.51	94.44	95.41	96.78	96.88	97.11	97.12
	安26识别率	65.93	73.77	92.80	88.03	89.63	92.30	92.63	92.90	93.07
	奖状识别率	94.60	91.00	96.27	98.47	99.20	99.50	99.50	99.63	99.57
	雅克42识别率	86.63	91.27	95.57	96.93	97.40	98.53	98.50	98.80	98.73
删除第1, 2个残差块	平均识别率	74.51	82.38	89.67	92.31	93.21	93.93	94.28	94.54	94.44
	安26识别率	58.33	69.17	78.00	81.87	83.13	84.60	85.57	86.27	85.93
	奖状识别率	68.50	80.40	92.97	96.60	97.93	98.33	98.53	98.57	98.57
	雅克42识别率	96.70	97.57	98.03	98.47	98.57	98.87	98.73	98.80	98.83
删除inception模块	平均识别率	78.08	84.42	90.10	92.81	93.71	95.14	95.26	95.48	95.58
	安26识别率	73.73	79.13	85.17	88.37	89.13	90.97	91.03	91.30	91.30
	奖状识别率	82.67	88.97	95.53	97.53	98.43	99.20	99.23	99.27	99.27
	雅克42识别率	80.83	85.17	89.60	92.53	93.57	95.27	95.50	95.87	96.17
删除降噪自编码器层	平均识别率	78.59	83.60	90.00	93.62	94.83	96.32	96.99	97.40	97.33
	安26识别率	71.97	76.90	83.63	88.80	90.37	92.30	93.33	93.73	93.53
	奖状识别率	77.73	83.07	91.20	94.67	96.33	97.70	98.63	99.07	99.17
	雅克42识别率	86.07	90.83	95.17	97.40	97.80	98.97	99.00	99.40	99.30

下载: 导出CSV

参考文献(17)

[1]	王俊, 郑彤, 雷鹏, 等. 深度学习在雷达中的研究综述[J]. 雷达学报, 2018, 7(4): 395–411. doi: 10.12000/JR18040 WANG Jun, ZHENG Tong, LEI Peng, et al. Study on deep learning in radar[J]. Journal of Radars, 2018, 7(4): 395–411. doi: 10.12000/JR18040
[2]	LUNDÉN J and KOIVUNEN V. Deep learning for HRRP-based target recognition in multistatic radar systems[C]. 2016 IEEE Radar Conference (RadarConf), Philadelphia, USA, 2016: 1–6. doi: 10.1109/radar.2016.7485271.
[3]	盖晴晴, 韩玉兵, 南华, 等. 基于深度卷积神经网络的极化雷达目标识别[J]. 电波科学学报, 2018, 33(5): 575–582. doi: 10.13443/j.cjors.2017112101 GAI Qingqing, HAN Yubing, NAN Hua, et al. Polarimetric radar target recognition based on depth convolution neural network[J]. Chinese Journal of Radio Science, 2018, 33(5): 575–582. doi: 10.13443/j.cjors.2017112101
[4]	JITHESH V, SAGAYARAJ M J, and SRINIVASA K G. LSTM recurrent neural networks for high resolution range profile based radar target classification[C]. The 3rd International Conference on Computational Intelligence & Communication Technology, Ghaziabad, India, 2017: 1–6. doi: 10.1109/ciact.2017.7977298.
[5]	徐彬, 陈渤, 刘宏伟, 等. 基于注意循环神经网络模型的雷达高分辨率距离像目标识别[J]. 电子与信息学报, 2016, 38(12): 2988–2995. doi: 10.11999/JEIT161034 XU Bin, CHEN Bo, LIU Hongwei, et al. Attention-based recurrent neural network model for radar high-resolution range profile target recognition[J]. Journal of Electronics &Information Technology, 2016, 38(12): 2988–2995. doi: 10.11999/JEIT161034
[6]	徐彬, 陈渤, 刘家麒, 等. 采用双向LSTM模型的雷达HRRP目标识别[J]. 西安电子科技大学学报, 2019, 46(2): 29–34. XU Bin, CHEN Bo, LIU Jiaqi, et al. Radar HRRP target recognition by the bidirectional LSTM model[J]. Journal of Xidian University, 2019, 46(2): 29–34.
[7]	刘家麒, 陈渤, 介茜. 基于注意力机制和双向GRU模型的雷达HRRP目标识别[J]. 雷达学报, 2019, 8(5): 589–597. doi: 10.12000/JR19014 LIU Jiaqi, CHEN Bo, and JIE Xi. Radar high-resolution range profile target recognition based on attention mechanism and bidirectional gated recurrent[J]. Journal of Radars, 2019, 8(5): 589–597. doi: 10.12000/JR19014
[8]	赵飞翔, 刘永祥, 霍凯. 基于栈式降噪稀疏自动编码器的雷达目标识别方法[J]. 雷达学报, 2017, 6(2): 149–156. doi: 10.12000/JR16151 ZHAO Feixiang, LIU Yongxiang, and HUO Kai. Radar target recognition based on stacked denoising sparse autoencoder[J]. Journal of Radars, 2017, 6(2): 149–156. doi: 10.12000/JR16151
[9]	张欢. 基于射频隐身的机载雷达系统软件实现及HRRP目标识别研究[D]. [硕士论文], 南京航空航天大学, 2016. ZHANG Huan. RF stealth based airborne radar system simulation and HRRP target recognition research[D]. [Master dissertation], Nanjing University of Aeronautics and Astronautics, 2016.
[10]	GOODFELLOW I, BENGIO Y, and COURVILLE A, 赵申剑, 黎彧君, 符天凡, 等译. 深度学习[M]. 北京: 人民邮电出版社, 2017: 199–226. GOODFELLOW I, BENGIO Y, and COURVILLE A, ZHAO Shenjian, LI Yujun, FU Tianfan, et al. translation. Deep Learning[M]. Beijing: Posts and Telecommunications Press.
[11]	HINTON G E, SRIVASTAVA N, KRIZHEVSKY A, et al. Improving neural networks by preventing co-adaptation of feature detectors[J]. arXiv:1207.0580, 2012.
[12]	NIELSEN M A. Neural Networks and Deep Learning[M]. San Francisco, USA: Determination Press, 2015.
[13]	HE Kaiming, ZHANG Xiangyu, REN Shaoqing, et al. Deep residual learning for image recognition[C]. 2016 IEEE Conference on Computer Vision and Pattern Recognition, Las Vegas, USA, 2016: 770–778.
[14]	SZEGEDY C, LIU Wei, JIA Yangqing, et al. Going deeper with convolutions[C]. 2015 IEEE Conference on Computer Vision and Pattern Recognition, Boston, USA, 2015: 1–9.
[15]	XING Chen, MA Li, and YANG Xiaoquan. Stacked denoise autoencoder based feature extraction and classification for hyperspectral images[J]. Journal of Sensors, 2016, 2016: 3632943.
[16]	DU Lan, LIU Hongwei, WANG Penghui, et al. Noise robust radar HRRP target recognition based on multitask factor analysis with small training data size[J]. IEEE Transactions on Signal Processing, 2012, 60(7): 3546–3559. doi: 10.1109/tsp.2012.2191965
[17]	RICHARDS M A, 邢孟道, 王彤, 李真芳, 等译. 雷达信号处理基础[M]. 2版. 北京: 电子工业出版社, 2017: 64–67. RICHARDS M A, XING Mengdao, WANG Tong, LI Zhenfang, et al. translation. Fundamentals of Radar Signal Processing[M]. 2nd ed. Beijing: Publishing House of Electronics Industry, 2017: 64–67.

施引文献

期刊类型引用(6)

1.	黄绮雯，张慧华，韩尚宇，纪晓磊. 弧形负泊松比蜂窝结构抗冲击性能的有限元分析. 南昌航空大学学报(自然科学版). 2024(02): 76-88 . 百度学术
2.	王波，丁逸凡，王坤，陆远，张青林. 基于ZYNQ的通用智能超表面控制系统设计. 激光杂志. 2024(10): 26-32 . 百度学术
3.	孔亚盟，王国玉，冯德军，王俊杰. 相位调制表面周期非均匀间歇调制方法. 兵工学报. 2023(04): 1209-1216 . 百度学术
4.	齐家瑞，姚斌，王均委，杨忠钰，穆帅虎，向泰. 可控超材料微波反应腔的加热效率及均匀性研究. 云南大学学报(自然科学版). 2023(04): 896-903 . 百度学术
5.	于博，庄书磊，王正心，王曼诗，郭兰军，李鑫煜，郭文瑞，苏文明，龚诚，刘伟伟. 基于纳米印刷技术的双螺旋太赫兹可调超表面. 物理学报. 2022(11): 345-352 . 百度学术
6.	张昭，张磊，郭江川，姚欣欣，谭治军. 声子晶体增材制造的研究进展. 电焊机. 2021(08): 11-22+176 . 百度学术