基于矩阵拟合的InSAS干涉相位估计方法

白生祥，张森，张学波，郑春弟

(1.海军陆战学院 教研部，广东 广州510431;2.海军工程大学 电子工程学院，湖北 武汉430033)

干涉合成孔径声呐(InSAS)技术是目前海洋探测领域的研究热点和前沿。作为一种新兴的水下三维高分辨成像技术，它兼备了合成孔径成像分辨率高和干涉测深精度高的优点，具有十分广阔的应用前景［1］。InSAS获取水下三维地形的原理与干涉合成孔径雷达(InSAR)相似，一般都需要经过合成孔径成像、复图像配准、干涉相位滤波、相位解卷、以及DEM重构等步骤［2-4］。在干涉处理过程中，这几个步骤依次级联，前一步骤存在的误差会累积、传递到下一步当中。由于复图像配准时要求精度达到亚像素级，一个较小的配准误差就可能导致错误的相位展开结果，给干涉处理带来了很大挑战。因此，研究对配准误差具有宽容性的干涉相位估计方法具有重要的意义。

近年来，干涉信号处理中提出的许多相位滤波方法［5-6］:如圆周期中值滤波、圆周期均值滤波、自适应滤波等，都要求图像精确配准，否则难以准确地恢复真实相位。文献［7］中提出了一种基于联合像素模型的子空间投影方法，能够在配准误差不超过一个分辨单元的条件下准确地估计干涉相位，为干涉信号处理开辟了一条新的途径。许多文献进行了研究改进［8-12］。这些方法主要集中于对联合像素模型的改进上，属于基于子空间分解一类的方法。而由于存在配准误差，这类方法在特征分解时存在信号子空间向噪声子空间扩散的问题，从而造成干涉相位的估计误差，困扰着估计性能的提高。本文基于矩阵拟合的思想提出一种干涉相位估计方法，所提方法采用矩阵拟合代替了特征分解，因而能够从源头上摆脱子空间扩散的影响，获得更加稳健的干涉相位估计性能。

1 联合像素模型

假定InSAS复图像对已精确配准，对同一像素i(对应于同一地面单元)，2个传感器阵列接收的复数据矢量s(i)可以写为

式中:a(φi)=［1 ejφi］T为阵列导向矢量;s1、s2分别为传感器1、2所对应的复图像数据;x(i)为像素i复幅度矢量(即不同传感器接收的复散射矢量);n(i)为加性零均值高斯白噪声;φi为真实地形的干涉相位，也是要估计的参数;⊙表示Hadamard积;s(i)为观测数据矢量，可以描述为联合零均值复周期高斯随机矢量，其对应的联合协方差矩阵Cs(i)为

式中:Rs(i)为像素i的相关矩阵，定义为

在实际情况中，2幅图像很难完全精确配准，当存在配准误差时，基于单像素模型的干涉相位估计方法将失效。在这种背景下，文献［7］提出联合像素模型，利用当前待估计像素及其邻近像素具有相同统计特征这一特性，通过综合利用邻近像素对的相干信息，从统计估计的角度，使干涉相位尽可能得到恢复。典型的联合像素模型构造中像素窗为(3×3，3×3)，此时，联合观测数据矢量s(i)表示为

对应的联合协方差矩阵Cs(i)为

假设联合观测数据矢量s(i)的相邻像素所处的局部区域为平地，则有

2 矩阵拟合相位估计

对子空间拟合的思想进行拓展，使得拟合的空间不仅是信号子空间或噪声子空间，而是一个与接收数据有关的矩阵。为了便于推导，首先从联合协方差矩阵的拟合说起，可以构造一个如下的拟合关系［13-14］:

式(7)给出的拟合问题是一个二维可分离的最小二乘问题，固定φi就可以求出的最小二乘解:

式中:†表示矩阵伪逆。对的求解涉及角度搜索问题和矩阵求逆运算，运算量巨大。为简化计算，选取一个对角矩阵作为T的解取代的最小二乘解，即

式中:τ1、τ2、…、τM是M个复变量。将式(10)代入式(7)，则联合矩阵协方差矩阵拟合可分解为式(11)所示的M个独立的列向量拟合问题:

结合式(6)分析可知，上式方括号内的项为与φi无关的常数，因而运算量可大大降低。将式(12)代入式(11)可得到关于相位估值的解，即

将式(14)代入式(7)，仿照式(9)～(13)的求解过程，得到i的解为

式中

将式(17)代入式(16)并化简后得

经过上面的一系列推导就完成了干涉相位的估计问题，但显然直接采用式(18)进行搜索计算量仍然十分巨大，这里提出一种快速计算方法，通过给出i的闭式解形式避免干涉相位的搜索问题。

为了便于推导，令

容易证明矩阵A为M×M维的Hermitian矩阵，以2N表示偶数M，将A表示为2×2维的子矩阵组成的分块矩阵形式，则式(19)可写为

式中:Amn为A的子矩阵。令

易知B同样为2×2维Hermitian矩阵。令，则

由式(23)可见，当μ+φi=2kπ(k为整数)即φi=2kπ-μ时，代价函数取得最大值。又因为φi为干涉相位，其取值范围为-π≤φi＜π，所以其估计结果i可表示为

3 矩阵拟合方法处理步骤。

采用矩阵拟合方法估计干涉相位可分以下4个步骤。

1)图像粗配准。

利用传统的配准算法(如相关法)对SAS图像进行粗配准。只要配准误差不超过一个分辨单元即可，因此图像配准的复杂程度可大大降低。

2)估计联合协方差矩阵。

构造联合观测矢量s(i)，计算相应的联合协方差矩阵s(i)。在满足各态历经性的条件下，可用样本协方差矩阵按下式进行估计:

式中:2L+1为获得的独立同分布的样本数。

3)估计加权相关系数矩阵。

实际中，相关系数矩阵也是通过样本平均得到，令sm(i)和sn(i)分别表示联合观测矢量s(i)中第m和n个样本，则相关系数可用下式估算:

4)矩阵拟合。

利用式(18)构造参数为φi的拟合关系式，使代价函数达到最大值时，通过式(24)求出所对应的干涉相位i即为最优估计。

对合成孔径声呐图像中的每个分辨单元按照上面的4个步骤逐一操作，就能够精确的恢复真实相位。

4 处理性能验证及分析

下面分别采用InSAS仿真数据和实测数据验证本文方法的有效性。

4.1 InSAS 仿真数据

采用文献［15］中提出的仿真方法生成回波信号，经合成孔径成像处理后，对图像进行相干处理得到干涉相位图。图1为2幅合成孔径图像精确配准时不同方法的处理结果，图1(a)为直接共轭相乘得到的干涉相位图，图1(b)为采用联合子空间投影方法估计的干涉相位图，图1(c)为采用本文方法估计的干涉相位图，图1(d)为2方法处理结果与真实相位的误差直方图。从图中可以看出，在精确配准条件下，2方法都能够正确估计出干涉相位，估计性能相差不大。

图2为2幅图像存在1个像素配准误差时不同方法处理得到的干涉相位。由图可见，此时联合子空间投影和本文方法仍可以较好的恢复干涉相位，两者相比，本文方法估计出的条纹更清晰，估计误差更小，性能更优。图3为不同配准误差条件下，2种估计方法所得结果的均方根误差变化情况。由图可见，当配准误差小于0.5时，2种方法均方根误差差别不大;当配准误差大于0.5时，本文方法呈现出更好的适应性和稳健性，能够保持更好的相位估计性能。

图1 精确配准时不同相位估计方法处理结果比较Fig.1 Comparison of interferogram for accurate coregistration

图2 存在1个像素配准误差时不同方法的处理结果比较Fig.2 Comparison of interferogram for the coregistration error of 1 pixel

图3 均方根误差随配准误差的变化Fig.3 Standard deviation versus coregistration error

4.2 InSAS 实测数据

在国家863计划(2007AA091101)的支持下，海军工程大学在原有基础上成功研制了干涉合成孔径声呐海试样机［16］，下面采用该样机2010年在试验过程中采集的实测数据，对提出的方法进行验证，样机系统主要参数设置见表1。

表1 InSAS样机系统参数Table 1 Parameters of the InSAS system

为了更直观的对2种方法在不同配准误差下的性能进行对比，选择配准误差范围为［0，1］个分辨单元的2幅干涉合成图像进行干涉处理，图像配准误差总体上沿着距离向逐渐降低。从图4可以看出，在干涉合成孔径图像配准精度较高的右侧区域2种方法处理效果差别不大，都能恢复清晰的干涉条纹。但在配准误差较大的左侧区域，联合子空间投影方法由于受到子空间扩散的影响出现大量离散点，相比而言，矩阵拟合方法估计性能更加稳健，具有更好的配准误差适应能力和更优的估计性能。

图4 InSAS实测数据处理结果比较Fig.4 Comparison of interferogram for the data of trial InSAS system

5 结束语

本文提出一种基于矩阵拟合的InSAS干涉相位估计方法，该方法利用了包含在联合协方差矩阵中的相干信息，并采用矩阵拟合的思想进行干涉相位估计，从根本上避免了困扰子空间方法的信号子空间扩散的难题。利用InSAS仿真数据和实测数据对矩阵拟合方法的干涉相位估计性能进行了验证，处理结果表明了所提方法的有效性。

［1］HAYES M P，GOUGH P T.Synthetic aperture sonar:a review of current status［J］.IEEE Journal of Oceanic Engineering，2009，34(3):207-224.

［2］孙大军，田坦.合成孔径声呐技术研究(综述)［J］.哈尔滨工程大学学报，2000，21(1):51-56.SUN Dajun，TIAN Tan.The study of synthetic aperture sonar(SAS)technique(Review)［J］.Journal of Harbin Engineering University，2000，21(1):51-56.

［3］张学波，唐劲松，钟何平.合成孔径声呐多接收阵数据融合CS成像算法［J］.哈尔滨工程大学学报，2013，34(2):240-244.ZHANG Xuebo，TANG Jinsong，ZHONG Heping.Chirp scaling imaging algorithm for synthetic aperture sonar based on data fusion of multi-receiver array［J］.Journal of Harbin Engineering University，2013，34(2):240-244.

［4］HUANG Haifeng，WANG Qingsong.A method of filtering and unwrapping SAR interferometric phase based on nonlinear phase model［J］.Progress in Electromagnetics Research，2014，144(3):67-78.

［5］于晓歆，杨红磊，彭军还，等.一种改进的Goldstein In-SAR干涉图滤波算法［J］.武汉大学学报:信息科学版，2011，36(9):1051-1054.YU Xiaoxin，YANG Honglei，PENG Junhuan，et al.A modified Goldstein algorithm for InSAR interferogram filtering［J］.Geomatics and Information Science of Wuhan University，2011，36(9):1051-1054.

［6］ZHANG Wenge，ZHANG Qin，ZHAO Chaoying，et al.Noise reduction for InSAR phase images using BM3D［J］.Chinese Journal of Electronics，2014，23(2):329-333.

［7］LI Zhenfang，BAO Zheng，LI Hai，et al.Image autocoregistration and InSAR interferogram estimation using joint subspace projection［J］.IEEE Transactions on Geoscience and Remote Sensing，2006，44(2):288-297.

［8］张晓玲，陈钦，韦顺军.基于MUSIC算法的Pol-InSAR相位估计方法［J］.电子科技大学学报，2011，40(5):652-657.ZHANG Xiaoling，CHEN Qin，WEI Shunjun.Estimation method for Pol-InSAR multi-interferometric phase based on the MUSIC method［J］.Journal of University of Electronic Science and Technology of China，2011，40(5):652-657.

［9］李海，吴仁彪，廖桂生.改进的联合子空间投影的InSAR干涉相位估计方法［J］.系统工程与电子技术，2013，35(1):45-52.LI Hai，WU Renbiao，LIAO Guisheng.Improved joint subspace projection method for InSAR interferometric phase estimation［J］.Systems Engineering and Electronics，2013，35(1):45-52.

［10］ZHANG Sen，XU Yanyi，TANG Jinsong.A robust estimation method of interferometric phase based on weighted subspace fitting［C］//The 2nd International Congress on Image and Signal Processing(CISP’09).Tianjin，China，2009:1-5.

［11］ZHANG Sen，TANG Jinsong，CHEN Ming，et al.Image autocoregistration and interferogram estimation using extended COMET-EXIP method［J］.IEEE Transactions on Geoscience and Remote Sensing，2010，48(12):4204-4218.

［12］LIAO Guisheng，LI Hai.Estimation method for InSAR interferometric phase based on generalized correlation steering vector［J］.IEEE Transactions on Aerospace and Electronic Systems，2010，46(3):1389-1403.

［13］王永良，陈辉，彭应宁，等.空间谱估计理论与算法［M］.北京:清华大学出版社，2005:152-156.WANG Yongliang，CHEN Hui，PENG Yingning，et al.Methods and algorithms of spatial spectrum estimation［M］.Beijing:Tsinghua University Press，2005:152-156.

［14］LIU N，ZHANG L R，LIU X，et al.Multibaseline InSAR height estimation through joint covariance matrix fitting［J］.IET Radar，Sonar and Navigation，2009，3(5):474-483.

［15］CHEN Ming，ZHANG Sen，TANG Jinsong.An interferometric synthetic aperture sonar raw signal simulation based on points-scatterer model［C］//International Joint Conference on Computational Sciences and Optimization.Sanya，China，2009:367-369.

［16］张森，唐劲松，陈鸣，等.干涉合成孔径声呐海试样机研制与试验［J］.声学技术，2012，31(2):167-173.ZHANG Sen，TANG Jinsong，CHEN Ming，et al.Development and sea trial of interferometric synthetic aperture sonar［J］.Technical Acoustics，2012，31(2):167-173.