面向InSAR稀疏控制点测图的同名点提取方法

姜丽敏 陈曙暄 向茂生

①(中国科学院电子学研究所微波成像技术国家级重点实验室 北京 100190)

②(中国科学院研究生院 北京 100190)

③(北京航天自动控制研究所 北京 100854)

1 引言

InSAR数据处理中，干涉参数(基线长度、基线角和干涉相位偏置)直接影响最终获取数学高程模型(DEM)的精度[1]。国内外关于这些参数偏差的标定技术，大部分采用航带内和航带间重叠分景的外定标方案，该方案存在两个问题：(1)大区域地形测绘时，需要在每景干涉像对中实地布放足量且分布合理的地面控制点(Ground Control Points,GCPs)。工作量大、作业效率低；而且野外采集GCPs受地形条件限制，存在某些测区如荒山、沼泽等难以实现GCPs布放；(2)忽略了相邻景的邻接关系，会造成相邻景重叠区域高程不一致。因此，各景单独外定标方法难以实现经济高效的大面积地形测绘。为推动我国InSAR技术的实用化，开展InSAR稀疏控制点测图方法研究具有较强的现实意义和理论价值。

所谓InSAR稀疏控制点测图，就是充分利用相邻影像之间的HPs(即具有同一地理位置的特征)并联合少量GCPs，对测绘区域内所有景数据的干涉参数进行整体标定。因此，HPs的自动提取显得尤为重要。目前，国内外在SAR影像自动配准方面关于同名点提取的研究虽有一些成果[2,3]，但这些方法无需考虑干涉技术的特殊问题，不适用于InSAR稀疏控制点测图，前期工作也证明了这一点。

SAR的侧视成像模式和相干成像机理使得自动提取HPs较为困难。首先，相邻影像重叠域的成像视角差异，造成重叠区存在灰度不一致、旋转变形、尺度缩放以及仿射变形；其次，同一分辨单元内各散射点随机散射信号相互叠加产生的相干斑噪声，降低了图像的质量，掩盖了图像的细节结构，很大程度上影响SAR影像的特征提取。这就要求相邻影像中提取的特征，必须对旋转、尺度缩放、仿射变形、视角变化等由于SAR的特殊成像机制造成的图像变化因素，保持一定的不变性和适应性；最后，基于InSAR的地形测绘是以SAR复数据提取的相位信息为信息源获得地形高程，所以最终匹配的特征所在位置处的干涉相位必须能够正确反映其对应的高程信息。

受上述条件约束，基于光学影像的特征提取算法难以直接应用于SAR图像。本文围绕不变特征检测，结合近年来计算机视觉领域关于局部不变量描述符的相关研究[4]，利用SIFT[5]和SURF[6]的不变特性，针对本文特殊应用对干涉相位质量的约束，提出了以质量图为导引的特征筛选方法。

2 SIFT和SURF简述

二者均可分4步：尺度空间[7]极值检测、特征点精确定位、主方向确定和特征向量生成。

2.1 极值检测

SIFT利用差分高斯(Difference of Gauss,DoG)近似拉普拉斯高斯(Laplace of Gauss,LoG)作为极值点检测的基础，而SURF借助于积分图像[8]利用箱式滤波器近似LoG作为Fast-Hessian极值点检测的基础。

(1)DoG极值检测 差分高斯尺度由图像高斯尺度空间表述中相邻尺度层函数相减得到，即

其中G(x,y,σ)是标准差为σ的高斯函数，I(x,y)为图像，k是常量，*为卷积运算。DoG极值检测就是检测DoG函数在尺度空间上的极值点。

(2)Fast-Hessian极值检测SURF是基于Hessian矩阵的行列式实现极值检测：

2.2 特征点定位

两者均采用文献[9]的定位方法，前者基于DoG函数，后者基于Hessian行列式。

2.3 特征描述符

两者都是对特征点邻域强度分布的描述，包括特征点主方向确定和特征向量构成。

(1)SIFT描述符

第1步 主方向确定。利用兴趣点邻域像素的梯度方向分布为其指定方向参数，像素点(x,y,σ)的梯度幅值和方向为

第2步 特征向量生成。由兴趣点邻域的梯度方向直方图，建立一个128维的特征向量。

(2)SURF描述符

第1步 主方向确定。特征点的圆形邻域以π/3为单元，将单元格内所有像点在水平和垂直方向的一阶Harr小波响应dxi,dyi(i=1,…,m)(m为单元格内像素点总个数)分别相单元格的局部方向向量。邻域内各单元Hy/Hx的最大值即为特征点的主方向。

第2步 基于Harr小波响应的特征描述符生成。利用积分图像计算Harr小波响应，进而构成64维特征向量。

3 同名点提取算法

3.1 算法的约束条件

(1)成像机理约束 同航带相邻影像属于同一飞行航迹，但机载平台受其适应性限制，同航带一致处理会使运动误差沿航向积累，造成运动误差的空变效应增强，最终导致成像质量下降，两幅影像相干性降低，不利于后续干涉处理，故干涉数据处理采用重叠分景方式，如图1(a)。各分景成像时对平台的运动补偿相互独立，具有各自的理想航迹线，因此同一飞行航迹的相邻影像重叠区间存在尺度缩放、小角度旋转、仿射变形等。

同航向相邻航带属于两次飞行航迹，成像几何见图1(b)。重叠区位于航带Ⅰ远距端，航带Ⅱ近距端。SAR侧视成像几何和斜距投影使得重叠区在航带Ⅱ的地距采样间隔大于航带Ⅰ，例如重叠度50%，近、远距和中心视角分别为30°,60°和45°，此时重叠测区在两航带的地距分辨单元满足2关系，进而造成相邻影像重叠区的几何畸变差异较大(图1(c))。图1(c)中，SAR影像的近距压缩效应使得重叠测区等地距间隔摆放的6个圆点在两幅影像的几何畸变不一致；当两航带的重叠测绘范围非常小时，该局部几何畸变将更加显著。SAR影像的回波强度是分辨单元内无数个独立散射体散射信号的相参叠加，而散射体的后向散射系数除与自身参数(表面粗糙度、几何形状、介电常数)有关外，还与SAR系统参数(波长、极化和成像视角)密切相关[10]，视角变化会造成相邻影像重叠区存在较大的视角灰度差异；此外，很难保证两次航迹完全平行，非平行航迹会使重叠域存在旋转变形。

图1 相邻景成像几何

图2 SAR影像中的叠掩、阴影、多路径效应

综上所述，两种情况都要求算法对仿射、灰度、旋转和尺度变化等具有一定的不变性和适应性。因此，本文采用具有不变特性的SIFT和SURF描述兴趣点的邻域纹理信息。

(2)干涉定标约束 InSAR稀疏控制点测图利用HPs具有相等高程作为重叠数据块的参考基准，因此必须减少HPs所在位置处的干涉相位噪声，使其能够正确反映地物的真实高程。相位噪声产生的4种因素[11]：(1)系统噪声；(2)两通道图像的失配；(3)地形形变；(4)几何关系模糊。前两种噪声在成像和图像配准时得以良好抑制；地形形变产生的干涉相位在双天线单航过干涉方式下可以忽略。几何关系模糊引入的相位噪声是本文考虑的重点，主要有：叠掩、阴影和多路径反射(图2)。

图2中点a,b的回波延迟相等，两点回波相叠加产生叠掩，叠掩处的干涉相位不能反映任一点的高程。地形坡度大于入射余角时入射波无法照射到的黑色区域在SAR图像上呈现为阴影，阴影区没有回波信号，信噪比非常低。ac边垂直于地表面，其散射形式复杂多样，会产生多路径效应[12]，典型的反射路径有：ac边直接反射—线①,ac和地表面的等效二面角反射器产生二次反射—线②,ac和地表面相互作用产生三次反射—线③，后两者的回波延迟路径较前者长。上述3种现象多见于地形起伏剧烈和强散射特性目标邻域，例如建筑物、高压塔等人工目标，因此 HPs要远离这类目标群，减少 HPs处相位噪声引入的高程误差，避免误差传递和累积。

几何关系模糊引入的上述诸多干涉相位噪声在干涉图中主要体现为干涉复图像对的相关性较低或不均匀，为此本文采用质量图导引方法避免其影响。

3.2质量图导引的特征筛选方法

3.1节的第(2)个约束条件要求对正确匹配点对的相位数据进行评价。文献[13,14]从不同角度分析论证了干涉相位标准差σφ与相关系数γ，视数N的关系，证明σφ的Cramer-Rao极限为式(4)。N确定时，σφ为γ的函数，图3为N=1,2,4,8,10,12,16,32时σφ随γ的变化，显然σφ随γ的增大而减小。基于此结论，本文以质量图为导引剔除位于相干性低的水域、阴影和相位噪声大的人工目标邻域的匹配点对，保证HPs的相位质量。

图3 视数N不同时σφ随γ的变化关系

复干涉像对u1,u2的相干质量图由式(5)计算，式中上标*为复共轭运算，E取均值。尽管σφ随γ呈现图3变化趋势，但单个像点所包含信息非常有限，相关系数大并不意味该像点的相位噪声小：首先，若强散射目标后向散射信号的旁瓣峰值大于其邻域弱散射目标的主瓣峰值，则强目标淹没弱目标，进而使邻域后向散射较弱的像素点对应的干涉相位偏离真实值；其次，后向散射方式相对复杂的人工物体尤其是硬目标，散射系数的动态范围达+40 dB[10]，成像几何差异造成其自身及邻域像点的γ在[0.10,0.99]间波动。因此，为避免邻域相关性差的点影响该点干涉相位的正确解缠，利用相关性评估相位噪声必须在像点的局部邻域内进行统计分析，需要确定3类参数——像点自身相关系数门限，邻域尺寸W和邻域相关性判定。

(1)邻域尺寸W确定根据成像几何和地物高度计算某点产生叠掩和多路径的最大范围W/2。图4中H为平台高度；地面点A1和树上某散射点A2到达平台的斜距相等，记为R;θ1,θ2分别为A1和A2的视角；x1,x2分别为A1和A2相对星下点的水平地面距离。由SAR图像的近距压缩效应(3.1节)知，相同高度地物在近距和远距对应的W1和W2满足W1＜W2。因此，根据图像远距的最高地物计算x1,x2进而确定邻域尺寸W=2(x2-x1)。

其中Q=1为一发双收，Q=2为自发自收，λ为波长，h,R1分别为像点的高程和其到主天线的斜距，B,α分别为基线长度和基线角，θ为视角。σφ0确定后，由式(4)得：

(3)邻域相关性 为得到判定邻域相关性的定量结论，以中国科学院电子学研究所在山西某地获取的X波段InSAR数据为例，分析3种典型邻域的相关性——不包含人工物体、与人工物体接边和包含人工物体，进而得到本文对邻域相关性的判定准则。图5中A,B,C区为尺寸相等且分别对应上述3种典型邻域的质量图以及相应的相关系数归一化直方图。

图4 邻域尺寸W确定

图5 平地和居民区的质量图及相关系数直方图

首先，图5(a)中A区，B区和C区各自质量图的均值和方差分别为 0.9741和 0.0290,0.9727和0.0318,0.9312和0.1211。数值表明这3种邻域其质量图的均值和方差均无显著变化，也就是说，根据兴趣点邻域质量图的均值和方差难以确定兴趣点是否受强散射特性目标影响。

其次，该3种典型邻域各自相关系数归一化直方图(图5(b)，图5(c)，图5(d))在γ∈[0.80,0.95]内的均值和方差分别为0.0795和0.1002,0.0869和0.1074,0.1475和0.1180。比较此3组数据可见，随着邻域逼近人工物体，尽管相关系数归一化直方图的均值和方差有所增大，但增长幅度非常有限且没有一个明显的门限。由此可见，根据邻域相关系数归一化直方图的均值和方差也难以确定兴趣点是否受强散射特性目标影响。

最后，就邻域像素点比例在相关系数γ＞0.95和γ∈[0.80,0.95]区间内的差值而言，A区，B区和C区却相差悬殊，分别为75.8109%,73.4821%和47.3833%。显然，包含建筑区的C区在这两个区间的像素点比例之差小于50%，而不包含建筑区的A区和与建筑区接边的B区均明显大于50%。这是因为裸地的后向散射特性较一致，不受邻近强点目标影响时InSAR两通道接收的SAR回波信号具有较强的相关性，而人工物体的后向散射往往具有很强的方向性，使得人工物体在两通道的回波信号去相关，进而使人工物体及其邻域的相关性波动非常大。

通过对大量平地和居民区的相关性分析，邻域质量图中γ＞0.95和γ∈[0.80,0.95]的像素点比例之差大于等于50%时，相关系数归一化直方图更接近理想情况，而邻域质量图的均值和方差却没有明显变化。因此，本文将邻域质量图中γ＞0.95和γ∈[0.80,0.95]的像点比例之差大于等于50%作为：邻域相关性判定准则。

3.3 算法流程

本文提出的两套 HPs提取算法(简称“SIFT-HPs”,“SURF-HPs”)，其基本流程为：

(1)利用SIFT或SURF检测多尺度特征并建立特征描述符；

(2)利用 Euclidean相似性度量距离并按 NN/SN准则[5]进行双向匹配：首先在第2幅影像的特征集合中寻找与第 1幅影像特征集合相匹配的点集，然后对中已匹配的点在A中寻找与其匹配的点，从而建初始匹配对集合；

(3)利用 RANSAC鲁棒算法[15]估计两幅影像的单应矩阵，初始匹配对中满足的内点称为“候选HPs”。具体步骤为：

(a)从第(2)步得到的初始匹配对集合中随机抽取4组；

(b)计算当前抽样所确定的单应矩阵和一致点集

(c)若当前的一致点集大于前一次的一致点集，则保留当前的一致点集和单应矩阵；

(d)由自适应算法终止抽样过程，获得最大一致集；

(e)由最终得到的最大一致集计算它所对应的单应矩阵，满足的兴趣点对则为“候选HPs”。

需要指出的是，阈值ε的取值越小，说明两个兴趣点的匹配精度越高。InSAR稀疏控制测图主要利用HPs具有相等高程这一信息，考虑到1 m范围内(米级分辨率的 SAR 图像上大约占3～4个像素点)自然地形的高程起伏通常非常小，因此文中ε=3 ；

(4)若“候选HPs”的相干系数γ≥，则判定W×W邻域内γ＞0 .95和γ∈[0.80,0.95]的像点比例是否相差50%。

4 实验及性能分析

为验证本文方法的有效性，采用中国科学院电子学研究所在陕西某地获取的干涉数据进行实验，实验区包含平地、道路、居民区以及玉米地等典型地形。

4.1 同航带实验

单景影像尺寸为3072×13384，航向重叠36%。重叠域远距端的最高物体建筑物约10 m，由式(6)计算W=10结合0.25 m地距采样间隔得质量图导引的邻域尺寸41 × 4 1，为保证HPs均匀分布采用分块处理。表1为从影像尺度空间表达的不同起始层开始时，SURF-HPs和SIFT-HPs分别在各阶段对应的匹配对数量，数值表明：

表1 同航带提取的HPs

图6 同航带第1分块相邻景HPs提取

(1)从第2尺度层开始检测有以下优点：(a)相当于对SAR图像的相干斑噪声进行了滤波处理，减小其对检测性能的影响，进而增加特征点的匹配稳定性；(b)本实验使用内存2 GHz的Intel(R)Core 2 双核 CPU，采用 Release优化版对程序进行优化，SIFT-HPs和 SURF-HPs计算速度分别提高8 s和97 s，而相应的 HPs分别减少1和0，说明从第 2尺度层开始既能够提高检测速度又不牺牲检测性能。

(2)如2.1节所述，SURF是基于梯度的梯度实现尺度空间极值检测，检测的兴趣点在SAR影像上梯度信息丰富的区域(如建筑物、地形接边处等)比较密集；加之基于Harr小波响应的特征描述符使得其对邻域梯度信息丰富的兴趣点的描述比较完备，造成其获取的匹配点对多集中在人工物体附近，因此SURF-HPs提取的HPs数量低于SIFT-HPs。

篇幅所限，仅将SURF-HPs在两分块的关键步骤示于图6和图7。图6(a)中圆圈代表第1分块相邻影像分别检测的2857和2644个兴趣点，连线指示初始匹配对，其中有一显著误配对(箭头所示)。图6(b)为满足的内点和质量图导引的HPs(非椭圆)。位于远距端的第2分块(图7)包含大量含玉米秆的田地，数据获取时当地3 m/s的风速吹动玉米叶子引起两通道回波去相关，致使该区域相关性较低，相关系数归一化直方图峰值为2.5970%且在γ∈[0.80,0.95]内的波动较大(均值和方差见表2)；另一方面，相同散射类型的地物在远距端的回波信号较之近距弱且信噪比差，使得图6(b)中A区和图7中B,C,D区所示裸地的相关性依次降低，各区相关系数归一化直方图在γ∈[0.80,0.95]的均值和方差以及质量图中γ＞0 .95和γ∈[0.80,0.95]的像点比例见表2，数据显示D区相干性稍逊于玉米地。上述两个因素造成SIFT-HPs和SURFHPs在该分块提取的 HPs非常有限，分别为2和0(见表1)。

图7中第11,13,15,17号点邻域质量图和相关系数归一化直方图的统计值见表2。11号点位于C区上方且邻近玉米地，理论上其邻域相关性应介于二者之间，实验数据与理论分析相吻合。位于B区右边但周围场景较复杂的13号点，其邻域相干性远低于B区。尽管15、17号点较11号点远离玉米地，但它们上方有一排东北走向的树且下方D区的相干性较C区差，因此这两点的邻域相干性均比11号点差。相比17号点，15号点左边有田梗，因此其邻域相关性比17号差。这4个点邻域以及上述4个区域的质量图均值和方差间的最大差异分别为0.0324和 0.0343都非常小，再次说明难以利用质量图的均值和方差判定邻域相关性。相反，它们的质量图中γ＞0 .95和γ∈[0.80,0.95]的像点比例之差的差异非常悬殊，最大和最小值分别为72.0430%(A区)和18.0250%(15号)且与自身邻域相关性相吻合，表明本文判定邻域相关性准则的合理性。

图7 同航带第2分块实验

表2 图6，图7部分区域及图7个别点对在前一景的相关性分析

表3 图7部分匹配点对的相关系数

图7中第7,10,20和21号匹配对在相邻干涉像对中的相关系数见表3。建筑区散射方向图的方向性非常尖锐[10]，其回波随姿态角变化的起伏要比典型地面回波随姿态角变化起伏快且显著，造成第20个匹配对在两幅质量图中的相关系数差异悬殊。尽管位于建筑物附近的第7个点对在两幅质量图中的相关系数均大于 0.95，但建筑物突出的角反射器效应通常会产生多路径，使得其周围回波在两通道的相关性波动较大，这在其邻域质量图图7(b)中显而易见，质量图中γ＞0 .95和γ∈[0.80,0.95]的像点比例分别为 4 6.7579%和 43.9619%，从而造成相关系数归一化直方图(图7(b))明显偏离理想形状，但邻域质量图的均值0.9173和方差0.0876却相当好，同样验证了难以利用质量图的均值和方差判定邻域相关性。第10号点对恰好位于相对X波段近似光滑的沥青马路上，散射回波类似镜面反射且对视角非常敏感，致使它在相邻干涉像对的相关系数仅为0.360327和0.543882。

4.2 相邻航带实验

影像尺寸分别为3072×13384和2048×7900，旁向重叠50%，航向重叠70%，分辨率分别为0.5 m和1 m。不同分辨率使得两幅影像重叠域间同一地物存在较大的尺度缩放和明显的仿射变形(见图8(a))；重叠域在两幅影像的成像视角范围分别为[4 5o,60o]和[3 0o,45o]，地物后向散射系数随视角增大而衰减[10]造成图8(a)左边影像的整体视觉灰度小于右边。

左边影像最远距的高层物体高压塔约20 m，根据本文方法得质量图导引的邻域尺寸为48×48。对比SIFT-HPs和SURF-HPs各阶段的匹配对数目(表4)知，两影像存在较大的几何畸变时 SIFT-HPs性能更佳。图8(a)为SIFT-HPs中满足的内点，方框为4对HPs。31,33号点在左右两干涉像对的相关系数分别为(0.967697,0.968168)和(0.985176,0.952455)，然而前者邻域质量图几乎不受低相关区道路的影响(图8(b))，而后者邻域质量图包含道路(图8(c))使得质量图中γ＞0 .95和γ∈[0.80,0.95]的像素点比例之差33.0430%远低于前者72.2181%，因此后者邻域相关系数归一化直方图图8(c)偏离理想形状。尽管13号点位于两块玉米地间的裸地上，但玉米秆的去相关和远距回波信噪比低等因素使其质量图中γ＞0 .95和γ∈[0.80,0.95]间的像素点比例之差为 - 27.0673%，进而使邻域相关系数归一化直方图呈现图8(d)。

图8 相邻航带实验

表4 相邻航带提取的HPs

5 结论

本文着重研究 InSAR稀疏控制点测图中的同名点提取，借助于SIFT和SURF特征不变描述符，提出两套质量图导引的方法，该方法能有效剔除相干性低、相位噪声大等难以用于干涉数据处理的特征点对。同航带相邻景和同航向相邻航带相邻景的大量实验验证了本文方法的有效性和可行性，从而降低了InSAR大区域测图的人工作业量。

[1]Madsen S,Zebker H,and Martin J.Topographic mapping using radar interferometry：processing techniques[J].IEEE Transactions on Geoscience and Remote Sensing,1993,31(1)：246-256.

[2]Suri S,Schwind P,Uhl J,et al..Modifications in the SIFT operator for effective SAR image matching[J].International Journal of Image and Data Fusion,2010,1(3)：243-256.

[3]戳涌光,王耀强.基于点特征的多源遥感影像高精度自动配准方法[J].遥感科学与技术,2010,25(3)：404-409.Zhai Yong-guang and Wang Yao-qiang.An automatic and high-precision registration method based on point features for multi-source remote sensing image[J].Remote Sensing Technology and Application,2010,25(3)：404-409.

[4]Mikolajczyk K and Schmid C.A performance evaluation of local descriptors[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2005,27(10)：1615-1630.

[5]Lowe D G.Distinctive image features from scale-invariant keypoints[J].International Journal of Computer Vision,2004,60(2)：91-110.

[6]Bay H,Ess A,Tuytelaars T,et al..Speeded-up robust features (SURF)[J].Computer Vision and Image Understanding,2008,110(3)：346-359.

[7]Lindeberg T.Scale-Space Theory in Computer Vision[M].Netherlands：Kluwer Academic Publishers,1994：31-162.

[8]Simard P Y,Bottou L,Haffner P,et al..Boxlets：a fast convolution algorithm for signal processing and neural networks[C].Proceedings of the 1998 Conference on Advances in Neural Information Processing System II,MIT Press Cambridge,MA,USA,1999：571-577.

[9]Brown M and Lowe D G.Invariant features from interest point groups[C].In British Machine Vision Conference,Cardiff,Wales,2002：656-665.

[10]Ulaby F T,Moore R K,and Fung A K.Microwave Remote Sensing Volume II Radar Remote Sensing and Surface Scattering and Emission Theory[M].United States of America：Addison-Wesley Publishing Company,1982：231-362.

[11]Rodriguez E and Martin J M.Theory and design of interferometric synthetic aperture radars[J].IEE Proceedings-F,1992,139(2)：147-159.

[12]Henderson F M and Lewis A J.Manual of Remote Sensing Volume 2[M].America：American Society for Photogrammetry and Remote Sensing,1998：705-732.

[13]Li F K and Goldstein R M.Studies of multibaseline spaceborne interferometric synthetic aperture radars[J].IEEE Transactions on Geoscience and Remote Sensing,1990,28(1)：88-97.

[14]Lee J S,Hoppel K W,Mango S A,et al..Intensity and phase statistics of multilook polarimetric and interferometric SAR imagery[J].IEEE Transactions on Geoscience and Remote Sensing,1994,32(5)：1017-1028.

[15]Fischler M and Bolles R C.Random sample consensus：a paradigm for model fitting with applications to image analysis and automated cartography[J].Communications of the ACM,1981,24(6)：381-395.