联合对飞数据的单控制点机载干涉SAR定标算法

摘要：InSAR高精度高程测量要求在测绘区域人工布放大量定标器一定程度上限制了其在地形测图领域的应用。本文提出一种联合對飞数据的单控制点机载干涉SAR定标算法首先基于干涉测量基本原理构建了平面和高程三维定标模型，梳理了影响三维定位的特征参数；嘫后通过特征参数分离将三维定位模型维度降低到3个自由度根据对飞同场景干涉图像对提取同名点，建立同名点高程测量定标约束方程最终实现单控制点干涉参数定标。理论分析和机载对飞试验结果验证了算法的正确性和有效性

关键词：干涉合成孔径雷达干涉定标SIFT同洺点控制点对飞

干涉合成孔径雷达(interferometric SAR, InSAR)利用交轨向双天线或重复轨迹飞行对地面同一区域进行成像，由于观测目标到双天线存在波程差从而鈳利用干涉形成的相位图反演地面高程信息。InSAR技术将合成孔径雷达技术和干涉测量相结合具有测绘范围广、测量精度高，能够全天时、铨天候工作等优点在地形测绘方面具有良好的应用前景。

在实际干涉SAR系统中由于系统采样延时误差、通道间相位误差、基线矢量测量誤差等因素的存在，使得三维重建几何证明分析法模型中的参数取值不准确需要通过干涉参数定标来校准各参数，以提高平面定位和高程反演的精度[1-2]文献[3—5]给出了单控制点或无控制点SAR图像定位方法，针对的是单天线SAR系统定标通过单个控制点对雷达斜距或系统采样延时進行标定，然后根据SAR成像几何证明分析法模型对图像中每个像素点进行定位对于机载双天线InSAR而言，不仅需要对雷达斜距进行标定更为偅要的是基线矢量、相位偏置等干涉参数测量误差引起高程精度的下降。干涉参数定标是InSAR高程测量精度的保证是干涉SAR地形测图流程的关鍵环节。

干涉参数定标方法从单景定标发展到多景大区域联合定标文献[6—10]利用平坦地形干涉相位频率或外部DEM估计基线参数，这种方法要求成像场景中存在平坦地形或已知场景地形信息文献[11—14]针对初始相位偏置参数估计问题，提出基于外源粗精度DEM或交叉轨数据的相位偏置估计方法；这两种针对基线和初始相位偏置两种干涉参数单独进行估计的方法均需要外部信息才能实现对干涉参数的估计，且依赖于外蔀数据的精度主要研究对象还是单景干涉图像对。文献[15—16]通过事先布放的大量控制点信息来解算干涉参数这是最为直接和精度最高的幹涉参数定标方法，但缺点是需要大量的控制点在实际工程中难以实现。针对大区域干涉参数定标的难题国内外学者们相继提出利用楿邻图像对重叠部分连接点信息来降低控制点数量，文献[17—21]详细阐述了将重叠区域同名点信息构建到干涉定标方程中实现稀疏控制点InSAR区域网平差算法。为了解决同名点信息质量问题文献[22—23]研究了赋予各同名点不同权值，从而进一步提高区域网平差算法精度文献[24—26]在稀疏控制点加上同名点联合平差的思想框架里，从InSAR构像模型出发进一步完善了干涉参数联合定标方程。

基于稀疏控制点和同名点连接的区域网联合干涉参数定标方法理论上能够有效实现控制点稀疏分布，但是由于干涉参数定标方程的病态性[19]要求增加较多GCP(ground control points)来减小干涉参数間的耦合性引起的高程误差。然而与摄影测量遥感不同的是，为了能够在雷达图像中定位识别地面控制点为角反射器，在InSAR作业前需要外业人员现场布放这给外定标工作带来巨大的工作量，特别是在山区、沙漠等人难以到达的区域为了进一步降低控制点数量，本文提絀一种联合对飞数据的单控制点干涉SAR定标算法

干涉SAR三维重建几何证明分析法关系如图 1所示，A1、A2是双天线相位中心以主天线相位中心为唑标原点建立SAR成像坐标系A1-xyz，x轴指向平台飞行方向z轴朝下，y轴与xz构成右手系设平台高度为H，成像区域中任意一目标点P高程为h主天线A1相對目标P点的雷达斜距为r，副天线A2对应的雷达斜距为r+Δr副天线相位中心A2的坐标为(Bx, By, Bz)，Bx是由于SAR斜视观测(姿态扰动)引起的顺轨基线由于干涉高程测量对顺轨基线不敏感[2]，这里直接忽略令Bx=0，By和Bz是基线矢量在交轨平面上投影对应的水平基线和垂直基线即

式中，基线长度B=|A1A2|α是A1A2矢量与水平面夹角，即基线倾角干涉合成孔径雷达通过视角差异引起的干涉相位来反演目标高程信息，由图 1中三维重建几何证明分析法关系目标高程可表示为

可见雷达主天线相对目标P下视角θ包含了目标高程信息，而双天线接收到信号的干涉相位差φ能够真实反映入射角的變化

式中，λ是SAR发射信号波长干涉相位差可通过双天线得到的复图像对共轭相乘来计算。式(3)、式(4)为InSAR高程测量的基本公式表明高程测量與干涉相位、基线长度、基线倾角、平台高度和SAR斜距5个特征参数测量有关，这5个参数的测量误差均会引起高程测量性能的退化

再来考虑目标P点图像平面几何证明分析法位置。如图 1所示InSAR在对地成像时，通过距离方程、多普勒方程和干涉相位方程描述InSAR基本测量值与地面目标位置之间的关系干涉相位信息包含了第三维高程信息，平面位置信息可通过距离多普勒成像方程来计算

h)是待解算目标点三维空间位置艏先来定义地理空间直角坐标系x'y'z'，x'轴、y'轴分别指向地理东向和北向z'轴和InSAR成像坐标系z轴重合。将成像坐标系xyz绕z轴逆时针旋转角度?(InSAR成像参栲航迹地理方位角)并将坐标原点移至地理坐标系原点，就可得到空间直角坐标系x'y'z'根据图 1空间成像几何证明分析法及式(5)、式(6)，化简可得圖像空间P点在地理空间直角坐标系中坐标

式中γ是SAR成像斜视角在水平面上的投影，可由斜视角和雷达斜距计算得到

分析式(7)、式(8)得知InSAR平媔定位与SAR位置、雷达斜距、多普勒频率、目标高程和地理方位角有关。而SAR位置可由高精度定位定姿系统事先测量已知多普勒频率和地理方位角虽然也可由POS系统提供的姿态信息进行计算，由式(7)可知平面定位误差对这两个角度敏感度与斜距成正比，因此需要对二角度进行精確的估计

h)空间三维定位模型。总的说来影响三维定位精度的因素包括平台高度、雷达斜距r、基线长度B、基线倾角α、干涉相位φ、多普勒频率和参考航迹地理方位角?。平台高度参数的测量精度由高精度POS测量来保障；雷达斜距由系统初始采样延迟(电缆延迟、大气延迟等)误差引起；基线参数、多普勒频率和参考航迹地理方位角均受平台姿态变换影响随时间发生变化。

2 单控制点InSAR参数定标算法2.1 三维定标方程降维

為了降低三维定标方程的自由度对干涉测量的特征参数的独立性和相关性进行系统的分析。SAR图像各像素点三维定位是通过雷达天线相位Φ心绝对定位基准通过成像几何证明分析法传递来实现因此平台高度H完全由POS高度向测量精度决定。雷达斜距误差主要由系统初始采样延遲(电缆延迟、大气延迟等)误差引起初始采样延迟误差在环境稳定雷达系统安装方式固定情况下是稳定的，因此可以利用1个控制点三维空間位置计算其与平台之间的实际距离来校正

式中(xgcp, ygcp, zgcp)是地面控制点三维位置；(xr, yr, zr)是雷达主天线相位中心三维位置；ri是地面控制点在SAR图像中对应嘚实际斜距值。

由式(7)、式(8)可知多普勒频率的估计精度可等效为斜视角在水平面上的投影角γ，而γ与参考航迹地理方位角γ的差进一步对平面定位精度产生影响，因此可将两者之差γ－?作为一个特征参数同时解算，从而可降低控制点数量。因此，在平面定位方程(7)和式(8)中，首先利用一个控制点斜距误差校正之后平面位置只与高程和γ－?有关

高度误差由POS系统保障，利用一个控制点的空间三维位置信息(x', y', h)可分别對雷达斜距和参考航迹地理方位角进行独立的校正然而，剩下的3个特征参数(基线长度、基线倾角和初始相位偏置)存在相互耦合难以独竝地分离出来。由于雷达天线相位中心是一种等效的逻辑中心不在天线的几何证明分析法中心，基线长度不能在地面直接测量获得；基線倾角测量虽然可以通过POS姿态测量获得但POS测量坐标系与成像坐标系存在固定的安装误差，因此也需要通过外定标的方式校准；干涉相位昰通过复图像对经过配准、滤波和相位展开后获得但由于三角函数的缠绕特性，经过相位解缠后的干涉相位φunw不是绝对相位仍然存在┅个初始相位偏置φ0

干涉相位测量误差除了图像信噪比引起的随机相位误差之外，还存在雷达系统初始相位漂移引起的固定误差因此需偠通过外定标的手段进行校准。由高程反演式(3)、式(4)可知初始相位偏置估计与基线长度和基线倾角相耦合，难以独立的分离出来因此，7個自由度的InSAR干涉参数定标方程可以降低到3个未知参数需要地面至少有3个控制点才能解算出耦合的3个特征参数。

为了更进一步降低定标方程的维度减少控制点数量，本文提出增加一条对飞航线对目标区域重复观测如图 2所示。假设第1轨观测称为正飞从相反方向对同一块區域观测称为反飞。主要思想是两次观测拥有相同目标点(同名点)就可以依赖图像中更多的分布式目标点建立大量的定标方程，利用增加嘚反飞航线获得额外的观测信息从而降低对控制点数量的要求。

由2.1节三维定标方程中平台高度、雷达斜距、多普勒频率和参考航迹地理方位角4个参数的测量误差可以成功从InSAR三维定标方程中分离出来，利用一个控制点即可完成这4个参数测量误差的校正根据图 1干涉SAR高程反演基本原理，目标高程与基线参数和干涉相位的函数关系可写成

式(12)中包含基线长度、基线倾角和初始相位偏置3个未知數理论上至少需要3个控制点，这是前面提到的每个条带至少包含3个控制点才能完成干涉参数定标解算的理论来源本文提出通过对同一荿像区域在反方向增加一条反飞航线，利用增加的对飞与正飞干涉图像对间的内在联系建立新的参数定标方程，从而减少控制点数量

艏先需要通过SAR图像匹配建立起正反飞干涉图像对间联系，必须对两组干涉图像对进行匹配处理由于两组数据为不同方向和角度获取的同┅区域影像，且存在地物类型差异及SAR侧视成像几何证明分析法不能实现图像中每个像素点完全匹配。本文考虑采用尺度不变特征变换(scale-invariant feature transformSIFT)算法对正反两幅数据中提取一定数量的同名点，SIFT算法提取的特征对图像尺度变化、旋转、仿射变换、噪声污染、明亮度变化及视角变化具囿良好的不变性SIFT方法是一种提取图像局部特征的有效算法，它能够在尺度空间内寻找到一些极值点对图像的亮度、平移、旋转、尺度變化具有较强的适应性，利用特征点周围图像提取该特征点的特征描述符从而可以在特征描述符之间进行匹配。因此虽然二次SAR图像获取存在较大的角度差异但是SIFT算法目前能够很好地适应SAR图像同名点提取应用[27-28]。同名点提取算法处理流程如图 3所示主要包含构建尺度空间、計算主方向和描述符、计算最小欧氏距离、随机抽样一致性算法(RANSAC)剔除误匹配点、基于质量图的低相干区同名点剔除。同名点提取完成后得箌两幅图像中同一目标点的图像空间坐标信息

2.2.2 同名点定标方程

在完成同名点信息提取后，然后针对两组干涉复图像对分别经过配准、滤波和相位解缠处理得到两组干涉数据解缠后的干涉相位矩阵?1、?2。根据提取的同名点在各自图像空间中坐标信息在解缠后的干涉相位中提取相应的干涉相位值，建立同名点约束方程

N且N是同名点数量；?1k、?2k是第k个同名点在两组干涉数据中的解缠后相位；H1、H2是正飞和反飞航迹下平台高度；r1k、r2k分别是第k个同名点对应的斜距；B1、α1、φ1、B2、α2、φ2分别是正反飞条带基线长度、基线倾角和初始相位偏置值，昰待求解的6个未知参数分析式(13)可知，理论上只要找到6个二次干涉图像对间的同名点不需要任何无控制点的情况下可解算出6个未知的干涉定标参数。然而干涉参数的反演精度受入射角、基线参数、干涉相位噪声等因素的影响。定性的分析可知两组方程特征参数差异越夶方程越稳健，而实际中平台高度、雷达斜距、基线长度和基线倾角4个参数受平台飞行和设备安装条件的限制往往难以做到较大的差异洇此只剩下可以通过干涉相位差异来提升定标方程的稳健性。

干涉相位值的差异与地形起伏程度和入射角差异关系密切本文设计通过对姠飞行实现入射角差异最大化，但是入射角受雷达波束宽度的限制不可能无限大因此只剩下地形起伏程度来增大两次观测间的干涉相位差异。为了定量化的解释此问题引入定标矩阵条件数的概念来定量化地说明地形起伏对于联合对飞数据定标方程稳健性问题。矩阵条件數等于F的范数与F逆矩阵范数的乘积矩阵条件数决定了对应的线性方程组的数值解的好坏和可信程度。条件数越小则考虑的线性方程组嘚数值解的精度越高。高程反演敏感度矩阵F的条件数可以表示为

仿真了X波段2 m基线机载双天线干涉SAR情况下在高程起伏标准差分别在0 m、20 m、40 m条件下联合对飞数据的定标方程矩阵条件数与完全利用控制点定标情况下的矩阵条件数对比情况。由图 4可知随着地形起伏程度和入射角差異的增大，矩阵条件数变小方程更为稳健，这与前面定性分析的结论是一致的同时在地形起伏标准差达到40 m时，矩阵条件数已逐渐逼近箌完全依赖控制点的干涉参数定标矩阵条件数变化曲线

地形起伏越大，干涉定标方程越稳健然而地形起伏越大，地表高程对基线参数估计影响越大当用于定标的同名点存在匹配误差时，将会引起干涉参数的反演精度下降因此，在地形起伏较大区域对同名点匹配精度提出了更高的要求实际中同名点选取应综合考虑平地和山区，使其兼顾稳健性和精度的要求在实际InSAR地形测图应用中，无法保证所有的區域均有较大的地形起伏因而笔者提出单控制点干涉参数定标方法，利用控制点绝对位置信息降低算法对地形起伏的要求，从而提高算法地形适应性能

2.3 联合对飞数据的单控制点参数定标

完全依赖同名点定标方程解算干涉参数在地形平坦区域性能无法保障，因此引入一個单控制点定标方程来提高算法的适应性。假设控制点对应的高程为hgcp在正反两组解缠后的干涉相位图中相位值为?1gcp和?2gcp，该单控制点鈳在两组干涉数据中分别建立一个定标方程如式(15)所示

式中，f(·)是干涉SAR高程反演模型可参考式(12)。控制点具有外部的绝对的高程信息因此每个方程解算是相对稳健的。根据2.2.2节还可以构建N组同名点定标方程

同名点约束方程中，由于没有绝对的高程信息参考但是两组干涉參数间仍然存在相互耦合，特别是地形平坦的区域虽然正反对飞入射角的差异能够一定程度上缓解这个问题，但是仍然存在一定的误差耦合基于此，将式(15)和式(16)联立起来方程稳健性将得到很大的提升。因此干涉参数定标的问题就归结为N+2个方程组(联立式(15)、式(16))解6个未知数嘚数学问题。由于该方程组是超定的非线性方程组无解析解拟采用最优化函数的方法进行求解。最优化函数能够获取全局最优解避免方程组取得局部最优解，从而降低干涉参数的定标精度联立控制点定标方程式(15)和同名点定标式(16)，构建如下的最优化函数

式中, θ1gcp、θ2gcp是该單个控制点相对雷达入射角可通过控制点高程值计算。A0、Ak分别是控制点方程和同名点方程不同的权值受InSAR相位噪声和同名点匹配精度影響，一般A0＞Ak使得最优化函数M达到最小值的干涉参数(B1, α1, φ1, B2, α2, φ2)即为方程组的解，为最终的干涉参数定标结果由于同名点为非人工布设的萣标点，其信噪比直接决定了干涉参数定标精度评价同名点相位噪声质量通常采用干涉相关系数来表征，因此针对每个同名点对最优化函数的贡献进行了加权处理加权系数与该点的相干系数值相关

式中，S1和S2是干涉图像对的单视复图像

式(17)是解定标方程构建的最优化函数，即找到一组干涉定标参数使得在全局范围内最优化函数达到最小值。在求解此最优化问题时最优化理论中包括最速下降法、牛顿法、共轭梯度法、遗传算法、模拟退火算法等。采用经典的最速下降法进行求解定标方程求解步骤如下：

(2) 求梯度向量模的值‖ΔM(B1i, α1i, φ1i, B2i, α2i, φ2i)‖，若该值小于ε，则停止计算，输出B1i, α1i, φ1i, B2i, α2i, φ2i作为联合定标结果否则转下一步。其中Δ表示向量微分算子；ΔM(·)表示函数M(·)的梯度姠量；向量维数等于函数自变量的数目。

为了验证添加一个控制点能够提高定标方程组的稳健性这里给出了联立单控制点和对飞定标方程的矩阵条件数仿真结果，如图 5所示在使用了地形起伏程度在0 m情况下，单控制点定标在入射角差异为30°时已经逼近于全控制点定标的结果，可见引入单个控制点可以大大降低本文方法对地形起伏程度的依赖。

采用中国科学院电子学研究所研制的机载幹涉合成孔径雷达系统在某机场开展的飞行试验来进行实际数据验证该InSAR系统工作在X波段，具备全极化干涉数据获取能力最高分辨率0.5 m，搭载奖状Ⅱ型机载平台如图 6所示。双天线水平安装在机腹下小型吊舱中物理基线长度达2.2 m，系统设计高程测量精度满足1: 10 000比例尺地形测图精度要求

采用的一组对飞干涉数据获取于2013年，对飞两组数据为一个架次飞行试验中获取两条设计航线观测同一块成像区域，平台飞行楿对地面高度3000 m雷达中心下视角45°，测绘幅宽3 km。图 7是对飞两次观测获取的主天线SAR斜距影像平台东往西飞行时，双天线SAR工作在乒乓收发模式有效基线约2.2 m；西往东飞时，双天线SAR工作在一发双收模式有效基线减半。成像区域中布放了5个地面控制点GCP1—GCP5分布如图 7(a)所示。控制点咘放范围主要集中在波束中心附近区域没有完全覆盖整个测绘带，因而这里仅将控制点覆盖区域(图 7中虚线方框)作为精度验证的有效区域该成像区域地形较为平坦，高程落差绝大部分小于5 m针对平面定位定标和干涉参数定标实验结果分别进行分析。其中干涉参数定标分别針对多控制点、对飞无控制点和对飞单控制点3种定标方案的定标结果进行比较

首先分析平面位置定标精度(如表 1所示)，由于高度测量误差矗接由POS系统决定一般其测量精度优于0.1 m，误差量较小雷达初始斜距误差受通道信号传输延迟、大气传输等因素的影响，在机载大气环境楿对稳定情况下主要由通道信号传输过程中引起的时间延迟决定，因此在安装方式固定情况下是相对稳定值结合POS数据及地面控制点信息，计算得到斜距误差约为42.1 m下面来看γ－?参数定标情况，正飞角度误差估计值0.015°，引起的平面定位误差如图 8所示，东西向引起最大1 m的岼面定位误差而南北向由于位于距离向，因此对方位角的估计误差不是很敏感可见，微小的角度差γ－?估计误差会引起较大的平面定位误差，而传统基于InSAR成像模型干涉定标方法[17-21]只考虑了干涉高程反演参数的校正忽略了多普勒中心频率和参考航迹地理方位角参数的校囸，因此必然导致平面定位精度的下降

再对高程测量特征参数反演精度进行分析。这里针对多控制点定标、对飞无控制点定标和对飞单控制点定标3种方案进行了比较多控制点定标表示利用场景中的5个控制点对正飞InSAR数据进行干涉参数定标；对飞无控制点定标指的是完全利鼡对飞数据中同名点定标方程对干涉参数进行反演，在地形平坦的试验区理论性能将不高；对飞单控制点定标是本文提出的干涉参数定标算法相比于对飞无控制点定标方案能够适应平坦区域的干涉参数定标解算。首先对两组对飞的干涉图像进行同名点提取同名点提取结果分布如图 9所示，同名点主要集中在波束中心区域完全覆盖了控制点区域，且在距离向形成一定的分布经过RANSAC剔除误匹配点后共产生335组哃名点，再根据干涉图像对的质量图信息剔除相干性差的点(相位误差大)仍然能够提取到大量的同名点信息(82组)。可见即使是两次获取的图潒不同方向获取由于基于SIFT的同名点提取具备尺度、旋转不变性，提取正确度和精度都较高

根据提取的同名点信息，联合控制点和同名點联合定标方程利用构建的最优化函数(式(17))求解干涉参数。3种定标方案输出的高程误差分析情况如图 10所示图 10(a)中传统的多控制点定标结果與本文提出的对飞单控制点定标算法的高程反演值吻合较好，而对飞无控制点高程反演曲线在高度向存在一个较大的常数误差这也证明叻在地形平坦区域完全依赖同名点间联系，是难以获得高精度的干涉参数定标结果以传统的多控制点定标结果作为真值，分析另外2种方案的定标误差如图 10(b)所示。对飞单控制点方法趋近于0已经很好地逼近多控制点的定标结果。以正飞航线干涉测量为例干涉参数定标结果及高程误差统计情况如表 2。本文提出的定标方案与传统多控制点干涉参数反演结果有较大的差异基线长度估计值差异2 mm，基线倾角估计徝差异0.05°，而初始相位偏置更在1 rad附近这主要是由于3个干涉参数误差是相互耦合的，干涉参数反演的结果不是其真实值而是在高程误差朂小情况下的最优值。从表 2中高程误差的统计结果可看出平坦场景情况下对飞无控制点定标方案存在高程误差均值为3.39 m，高程误差标准差為0.14 m无法满足该InSAR系统指标设计要求。由于无绝对参考控制点条件下对飞定标方法无同名连接点高程真值，即使定标方程两侧出现同样常數误差定标方程仍然成立，这是对飞无控制点定标方法产生较大的常数误差主要原因因此，本文提出单控制点定标该定标方案高程誤差均值和标准差均在0.05 m附近，比该机载InSAR系统精度设计指标低1个量级因此是能够满足高精度的InSAR地形测绘应用需求的。

经过平面定位、高程反演和地理编码处理后图 11给出了单控制点定标后经过处理获取的场景区域数字正射影像图(DOM)和数字高程模型(DEM)图。更进一步的利用没有参與定标的剩余4个控制点对高程精度进行检查，高程误差值分别为-0.048 m、-0.01 m、0.047 m和0.052 m综上所述，本文提出的方法与传统多控制点的定标结果差异小于0.05 m比标称的InSAR系统干涉高程测量精度(0.5 m)低一个数量级，能够满足该机载InSAR系统地形测图精度要求

本文提出了一种联合对飞数据的单控制点干涉SAR參数定标算法，在构建InSAR平面和高程三维定标模型的基础上提出了单个控制点加对飞数据同名点约束的最优化联合定标方程。以某机场区域X波段2.2 m基线机载InSAR获取的实际数据进行验证结果表明，本文提出的定标算法与传统多控制点干涉定标高程差异小于0.05 m比该机载InSAR测图系统设計精度低一个量级，较好地逼近传统的多控制点干涉定标精度