干涉合成孔径声呐(interferometric synthetic aperture sonar，InSAS)^[1-2]是一种高分辨率的三维成像声呐，它由在方位向上平行安装的两个合成孔径声呐^[3-5](synthetic aperture sonar，SAS)组成，具有成像分辨率与工作频率和成像距离无关的优点^[6-7]。InSAS和干涉合成孔径雷达(interferometric synthetic aperture radar，InSAR)的原理是一样的，它们的工作流程主要包括合成孔径成像、复图像配准^[8]、干涉图生成、相位滤波、相位解缠和高程估计。InSAR中的一些信号处理的方法也被应用到InSAS中，但InSAS和InSAR也存在一些差别，主要体现在: (1)InSAR大多采用重复轨道单天线模式，即在不同的时间、不同轨道上对同一地区进行成像；而InSAS均采用单轨道双天线模式，即在声呐平台上以一定基线长度为间隔安装两条接收阵，一次经过目标区域即可获得两幅声呐图像，图像对之间的偏差相比于InSAR图像要小，因此，InSAS复图像配准要比InSAR复图像配准容易。(2)InSAR波长一般为厘米级，而基线长度通常是波长的数千倍，因此，干涉条纹非常密集。而InSAS的基线长度一般是波长的十几倍，干涉条纹比较稀疏，即干涉图相位变化相对比较缓慢。

在InSAS信号处理中，复图像的精确配准是获取高精度的数字高程模型(digital elevation model，DEM)的保证，配准质量的好坏直接影响到后续的干涉图生成和数字高程重建。复图像配准^[9-10]的基本步骤包括选取同名点对、几何变换模型构建^[11-12]和输入图像重采样。当两幅图像间的同名点对确定以后，就可以利用同名点对的坐标和偏移量进行几何变换。在InSAR几何变换模型构建中应用最广泛的是整体多项式模型，其原理是利用已知同名点的坐标和偏移量，通过最小二乘法进行曲面拟合，求得多项式的系数，该模型简单易于求解。在InSAR系统中，发射阵到目标的距离(斜距)一般是几百千米，成像区域的范围相比斜距要小得多，因此，用整体多项式曲面拟合就能达到较好的效果。但在InSAS系统中，斜距一般为几十米到几百米，成像区域的范围和斜距基本属于同一个量级，所以利用整体多项式曲面拟合很难达到较好的效果。在InSAR图像配准方法中，研究比较多的是基于特征的配准^[13-15]，但是InSAS图像不同于InSAR图像，有时水下地形特征没有InSAR那么明显，所以有时基于特征的配准在InSAS中并不适用。

针对前面提到的问题，本文提出一种基于分段曲面拟合的InSAS复图像配准方法。在InSAS复图像配准中还未见分段配准的思想。首先，将整个图像沿着距离向分成几段，相邻分段之间有部分重叠区域；其次，对每个分段，利用该段内的同名点的坐标和偏移量单独进行多项式曲面拟合，得到该分段的偏移量曲面；然后，将每个分段曲面拟合的结果拼接成整体的偏移量曲面；最后，利用偏移量对辅图像进行插值重采样，即可得到配准后的图像。利用仿真数据和真实的湖上试验数据将本文方法和传统的整体曲面拟合方法进行对比，实验结果表明，本文方法配准精度更高，更适合于大场景InSAS复图像配准。

1 分段曲面拟合方法

本文主要讨论InSAS复图像配准中的几何变换模型构建环节。由于InSAS采用的是单轨道双接收阵工作模式，上下两条接收阵平行安装，且声呐平台沿着方位向直线运动，因此，复图像配准过程中主、辅图像方位向偏移量可忽略不计，只考虑斜距方向偏移量。

1.1   多项式几何变换模型

本文选取同名点是根据相干系数最大原则，具体操作过程是: 在主图像上选取一定大小的窗口(如$ M\times N $)，在辅图像的搜索区域内(如$ k\times l $，$ k>M $且$ l>N $)进行滑动，每滑动一次，计算一次相干系数的值。当在整个搜索区域内滑动完以后，计算出相干系数的最大值，最大值点就是要找的同名点。对整个图像中的同名点进行筛选，保留相干系数值大于某阈值的同名点作为有效同名点。

设有效同名点的个数为$ n $，其坐标为$ ({x}_{1}, {y}_{1}) $，$ ({x}_{2}, {y}_{2}) $…$ ({x}_{n}, {y}_{n}) $，斜距的偏移量为d(x₁，y₁)，$ d({x}_{2}, {y}_{2})\cdots d({x}_{n}, {y}_{n}) $，则多项式几何变换模型为:

式中，$ {a}_{ij} $为待求的多项式系数；$ d(x, y) $为点$ (x, y) $处的偏移量值；$ N $为多项式阶数，一般二阶多项式能满足要求。令$ N=2 $，式(1)可以写成:

式(2)是$ n $个方程6个未知数的方程组，至少需要6个同名点的数据才能解出未知系数。一般情况下$ n $远大于6。令

则式(2)可以写为:

利用最小二乘法可以求得式(4)的解为:

将式(5)代入式(1)可以得到多项式曲面:

将整个图像中所有点的坐标代入式(6)可以求出每个点的斜距偏移量，这些偏移量构成一个二维曲面。每个同名点配准误差$ d{(x, y)}_{\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}} $为拟合的偏移量$ d(x, y) $与其真实的偏移量$ d({x}_{i}, {y}_{i}) $的差，即$ d{(x, y)}_{\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}}=d(x, y)-d({x}_{i}, {y}_{i}) $。

1.2   分段曲面拟合模型

1) 将整个图像沿距离向分割成几部分，如图‍1所示。图 1中，实线为分段线，两个相邻的分段之间设置一定的重叠区域，如虚线所示。本文采取的距离向分段策略是均匀分为5段，相邻分段之间的重叠区域设置为本段像素数的20%。

图  1  分段示意图

Figure  1.  Sketch Map of Partition of the Whole Image

下载: 全尺寸图片 幻灯片

2) 搜索位于每一段(包含重叠区域)内的有效同名点，利用这些同名点的坐标和偏移量进行多项式曲面拟合，得到每段的偏移量曲面。

3) 将每段拟合后的偏移量曲面拼接成一个整体的偏移量曲面。拼接时，关键是要处理好重叠区域像素点的偏移量值。对非重叠区域内的像素点，直接将该分段的拟合偏移量的结果作为最终的整体偏移量曲面的结果，而对于重叠区域中的像素点，如图 1中的点A，采用加权线性组合的方式计算其最终的偏移量值，计算公式为:

式中，$ f\left(A\right) $是待求的点在整体曲面中的偏移量值；$ {f}_{1}\left(A\right) $、$ {f}_{2}\left(A\right) $分别是点$ A $在第1分段和第2分段拟合的曲面中的偏移量值；$ {w}_{1} $和$ {w}_{2} $是归一化权重，$ {w}_{1}={r}_{2}/{l}_{1} $，$ {w}_{2}={r}_{1}/{l}_{1} $，$ {l}_{1}={r}_{1}+{r}_{2} $，$ {l}_{1} $是重叠区域的宽度(像素数)，$ {r}_{1} $是点$ A $到其所在的重叠区域左边界的距离(像素数)，$ {r}_{2} $是点$ A $到其所在的重叠区域右边界的距离(像素数)。

2 实验结果与分析

2.1   仿真实验

仿真的场景是平地上沿距离向并排放置的两个圆锥，其幅度图如图 2所示。用传统的整体曲面拟合法和本文的分段曲面拟合法进行配准，配准后的干涉图分别如图 3(a)和图 3(b)所示。通过比较可以发现，图 3(b)比图 3(a)的条纹更清晰，质量更好，特别是左侧部分的干涉条纹区别更明显。

图  2  仿真场景幅度图

Figure  2.  Amplitude Image of Scene Simulation

下载: 全尺寸图片 幻灯片

图  3  配准后干涉图

Figure  3.  Interferograms After Registration

下载: 全尺寸图片 幻灯片

图 4给出了两种方法的同名点的配准误差统计直方图。从图 4中可以看出，分段曲面拟合的配准误差主要集中在-0.2~0.2像素，而整体曲面拟合的配准误差主要集中在-0.4~0.4像素，由此可见，分段曲面拟合的配准误差明显小于整体曲面拟合的配准误差。

图  4  同名点的配准误差直方图

Figure  4.  Histogram of Registration Errors of Control Points

下载: 全尺寸图片 幻灯片

为了更直观地比较两种方法的配准效果，统计图 3的干涉图中间部分100×100像素的矩形区域内的残余点数目。在干涉相位图中，当4个相邻的像素点构成的环路中相位差的和不为0时，将左上角的点记为残余点，残余点示意图如图 5所示。

图  5  残余点示意图

Figure  5.  Sketch Map of Residue Points

下载: 全尺寸图片 幻灯片

残余点的具体计算公式为:

式中，$ \Delta 1 $为$ \phi (i, j) $与$ \phi (i, j+1) $之间的相位差的缠绕值；$ \Delta 2 $、$ \Delta 3 $、$ \Delta 4 $分别是如图 5所示对应点之间的相位差的缠绕值；$ W $表示缠绕运算(模$ 2\mathrm{\pi } $运算)；$ \phi $为干涉图中的缠绕相位值；若$ q\ne 0 $，将左上角的点$ (i, j) $记为残余点。配准质量越好，干涉图就越清晰，残余点就越少；反之，若配准不精确，干涉图中残余点就越多。

图 6(a)和图 6(b)分别是用整体曲面拟合法和本文分段曲面拟合法进行配准后的局部残余点分布图，残余点数目分别为1 743和777。

图  6  局部残余点比较

Figure  6.  Comparison of Local Residue Points

下载: 全尺寸图片 幻灯片

可以看出，图 6(b)中的残余点比图 6(a)的残余点少很多，由此可见，分段曲面拟合的配准结果要比整体曲面拟合的结果好。为了定量比较两种拟合方法的配准效果，选择均方根误差(root mean squared error，RMSE)、相关系数、残余点数目作为评价标准，结果见表 1(RMSE为无量纲)。从表 1可以看出，本文分段曲面拟合法的配准质量明显高于整体曲面拟合法。

表  1  仿真数据配准结果比较

Table  1.  Comparison of Simulation Data Registration Results

拟合方法	RMSE	相关系数	残余点数
整体曲面拟合	0.305 6	0.939 7	78 746
分段曲面拟合	0.068 7	0.967 1	28 173

下载: 导出CSV 

| 显示表格

2.2   真实湖上试验数据处理结果

选择大场景的真实湖上试验数据进行处理，该地形有沟壑、田埂，地形变化比较明显。采用的湖上试验数据参数为: 声速为1 500 m/‍s，载频为150 kHz，带宽为40 kHz，脉宽为0.02 s，脉冲重复间隔为0.2 s，距离采样频率为80 kHz，阵元大小为0.08 m。湖上试验数据幅度图如图 7所示。图 8给出了两种方法的同名点的配准误差直方图，可以看出，本文方法的配准误差更小。图 9(a)和图 9(b)分别是整体曲面拟合配准后和分段曲面拟合配准后的干涉图。

图  7  湖上试验数据幅度图

Figure  7.  Amplitude Image of Lake Trial Data

下载: 全尺寸图片 幻灯片

图  8  同名点的配准误差直方图

Figure  8.  Histogram of Registration Errors of Control Points

下载: 全尺寸图片 幻灯片

图  9  配准后干涉图

Figure  9.  Interferograms After Registration

下载: 全尺寸图片 幻灯片

由图‍9可以发现，在局部地形变化明显的地方，分段曲面拟合的配准结果更好，条纹更清晰，如图 9中的白色矩形框内的部分。

为了更清晰地比较两种拟合方法的配准效果，将图 9中两个白色矩形框内的部分进行局部放大。图 9(a)和图 9(b)中较大的矩形框局部放大结果分别如图 10(a)和图 10(b)所示，可以看出，图 10(b)的条纹比图 10(a)更清晰，由相位反映的地形变化更明显。

图  10  图 9中较大矩形框放大图

Figure  10.  Enlarged Image of Bigger Rectangle in Fig.9

下载: 全尺寸图片 幻灯片

图 11(a)和图 11(b)分别是图 9(a)和图 9(b)中左下角的矩形框的局部放大图，可以看出，图 11(b)的左侧第一个条纹比图 11(a)的左侧第一个条纹清晰很多。由此可见，在局部地形变化明显的地方，实用分段曲面拟合法进行配准比实用整体曲面拟合法进行配准更有优势。

图  11  图 9中左下角矩形框放大图

Figure  11.  Enlarged Image of Left-Down Rectangle in Fig.9

下载: 全尺寸图片 幻灯片

对两种方法定量的评价如表 2所示。由表 2可以看出，无论是均方根误差，还是相关系数和残余点数目，分段曲面拟合法的配准结果都优于整体曲面拟合法的配准结果。

表  2  真实数据配准结果比较

Table  2.  Comparison of Real Data Registration Results

拟合方法	RMSE	相关系数	残余点数
整体曲面拟合	0.198 4	0.580 2	81 956
分段曲面拟合	0.119 1	0.582 4	69 455

下载: 导出CSV 

| 显示表格

3 结语

本文针对大场景InSAS图像用传统的整体多项式曲面拟合法配准精度低的问题，提出了一种基于分段曲面拟合的InSAS复图像配准方法。将大场景图像沿距离向分成多段，每段单独进行多项式曲面拟合，再把各段拟合的结果拼接成整体的偏移量曲面。将本文方法与传统整体曲面拟合方法进行比较，并以均方根误差、相关系数、残余点数目作为评价标准，通过仿真实验和真实的湖上试验数据处理实验，验证了本文的分段曲面拟合配准法比整体曲面拟合配准法的配准精度更高。后续将主要在配准效率上进行并行化方法研究，以满足实时处理的要求。