利用机载双天线InSAR数据生成高精度DEM的试验研究——以大面积丘陵地区为例

张延冰，郭华东，韩春明

(1.中国科学院遥感与数字地球研究所，北京 100094;2.中国科学院大学，北京 100049)

0 引言

合成孔径雷达干涉测量(interferometric synthetic aperture Radar，InSAR)技术具有全天时、全天候，测量精度和效率高等特点，已成为当前快速获取高精度地面三维信息的重要手段之一[1]。在热带雨林以及地形复杂的高山、极地地区，InSAR已成为获取DEM的重要技术手段。双天线InSAR系统具有时效性和自动化程度高等特点，且不存在时间失相关问题[2-3]，在地形制图方面有着巨大的优势。自1974年Graham利用机载InSAR技术获取了能满足1∶25万比例尺地形图要求的高程数据以来，该技术越来越多地应用到地形测量中，并被证实了具有获取高精度DEM的能力[3-4]。

在科技发达国家，机载InSAR技术已被应用于地形测绘方面，用以提供业务化运行的服务和产品。如InterMap公司的STAR-3i系统已实现商业运行多年，实现了巨大的经济效益。国内机载InSAR系统研究相对较晚[5]，在测绘领域的应用仍处于起步阶段。2004年，中国科学院电子学研究所成功研制了“X波段机载双天线InSAR系统”[6]，并进行了地形测绘飞行试验，在地势平坦地区取得了较好的效果。但在地形复杂区域，由于缺乏精确的干涉参数和影像外方位元素，在数据处理时遇到了干涉质量、影像纠正等难题[7]，所获得的DEM还没有达到较高的精度。

本文在已有应用机载InSAR数据生成DEM的技术流程[8]基础上，引入精确干涉参数定标和区域网平差处理技术，提出利用国产机载双天线SAR数据生成丘陵地区大面积高精度DEM的方法。

1 试验区概况及试验设计

1.1 试验区概况

江油试验区位于四川盆地西北部及涪江上游，属于四川省绵阳市;地貌以低山丘陵为主;海拔高度383～717 m。测区包含居民地、河流、耕地、林地、铁路以及公路等地物类型，东西长59 km，南北宽4.6 km。经纬度范围介于 E104°30′00″～105°07′30″，N31°37′30″～ 31°40′之间。

1.2 机载InSAR系统及飞行参数

试验采用X波段机载InSAR系统。该系统频率为9.6 GHz，HH极化、乒乓模式;基线长度为2.33 m，脉冲重复频率(PRF)为3 125 Hz，脉冲宽度为16 μs，平均功率约220 W;获取的InSAR数据地面分辨率优于0.5 m。飞机上安装的高精度位置姿态测量系统(POS/AV510)能获取传感器的空间坐标、旋转角、速度及加速度等信息，参数如表1所示。

表1 飞行参数Tab.1 Flight parameters

1.3 野外数据获取

为实现高精度区域网平差，便于成果精度检验，需要在试验区精密布设角反射器控制点。根据任务需要及实际地形情况，在试验区两端及中部沿距离向分别布设3列(共28个)控制点。

另外，为实现高精度机载InSAR系统参数定标，在测区中部北侧沿距离向均匀布设9个定标控制点。试验区范围及控制点布设如图1所示。

图1 试验区范围及地面控制点分布Fig.1 Location of study area and ground control points distribution

外业实地踏勘并选取控制点，记录控制点点名和点位略图，用GPS设备组网静态观测控制点。机载雷达飞行获取数据时，在控制点上精确布设铝制3面等腰直角角反射器(边长35 cm)。差分GPS基准站采用测区中央双基站布放，利用高精度双频GPS接收机进行连续同步观测。基站在飞机起飞前0.5 h开机，飞机降落后0.5 h关机。

2 机载双天线InSAR获取DEM原理

图2为机载双天线InSAR干涉测量成像示意图。

图2 机载双天线InSAR干涉成像示意图Fig.2 Imaging geometry of airborne dual－antenna InSAR

式(2)—(4)中:c为光速;τ为绝对时间延迟;ρ为距离向采样率;j为P点方位向影像坐标;λ为电磁波波长;φ为该点的相位。

3 数据处理流程

采用自主开发的机载InSAR地形制图处理系统InSARmap软件处理数据，其流程如图3所示。

图3 机载双天线InSAR数据生成DEM处理流程Fig.3 Data processing flowchart of airborne dual－antenna InSAR DEM generation

3.1 干涉参数定标

干涉参数定标是高精度测图的关键技术，也是测图误差控制的关键步骤。干涉参数定标就是要对系统设备的相位偏移和干涉参数偏差进行严格校正，以提高生成DEM的精度[9]。

本文基于敏感度方程模型进行干涉参数定标。首先通过目标点高程方程对待定标的干涉参数求微分，得到敏感度矩阵，建立高程误差对干涉参数偏差的关系;再利用GPS测量定标控制点高程，用最小二乘法求解各个干涉参数偏差估计值，进而校正相应的干涉参数。试验选用基线长度B、基线倾角α、绝对延迟时间τ和相位偏移量φoffset4个参数建立敏感度方程，即

高程差和干涉参数偏差的关系为[10]

式中:Δ为L×1高程误差向量;F为L×N敏感度矩阵;△X为N×1待定标的干涉参数偏差向量;L为定标控制点数目;N为待定标参数个数，此处N=4。采用迭代方法进行干涉参数定标，用初始的干涉参数和相位生成DEM，通过解算敏感度方程得到各干涉参数偏差值，校正相应干涉参数，并重新计算DEM，直到2次计算的DEM均方差小于0.01为止。

3.2 干涉处理

干涉处理步骤包括复图像配准、干涉纹图生成、干涉纹图滤波以及相位解缠等。本研究采用基于快速傅里叶变换(fast fourier transform，FFT)的复相关法[8]进行复图像精配准，配准图像对干涉处理后生成干涉纹图和相干图，并采用中值滤波和Goldstein滤波相结合的滤波方法对复干涉纹图进行滤波处理，采用最小费用流法对滤波后的干涉纹图进行相位解缠，获得影像上每点的相位数据。通过InSAR-map软件中的干涉处理模块可实现整条带数据的自动批量干涉处理，能节省大量处理时间。

图4 区域网连接点分布示意图Fig.4 Control points distribution in bundle adjustment

3.3 区域网平差

由于POS系统测量的航偏角及雷达系统提供的多普勒中心频率精度不够高，在测区稀少控制点条件下仅依靠POS数据及干涉定标参数难以实现大面积、多航带SAR影像的高精度校正。区域网平差技术可以建立影像之间的约束关系，通过调整和精化部分外方位元素，实现稀少控制点条件下的多航带SAR影像的高精度定位[11]。

本文基于 Range-Doppler-Phase(RDP)模型[11]建立载机位置、像点及地面点之间的关系，确定平差参数系数矩阵，建立高精度三维模型误差方程组。利用地面控制点和连接点数据、定标数据、POS数据和相位信息，按最小二乘原理求解各景影像准确的外方位元素(平差参数)和连接点的地面坐标。确定的平差参数包括解缠相位φ偏移量、多普勒中心频率fDC及偏航角θy。由于高程和平面坐标对平差参数敏感程度不同，将高程和平面平差参数分开迭代求解。首先求解高程平差参数(解缠相位φ偏移量改正)，再将改正的高程平差参数代入平面误差方程，求解平面平差参数(多普勒中心频率fDC偏移量，偏航角θy改正数)。

在试验区4个条带内及条带间影像重叠区域选取6 254个连接点，组成区域网，如图4所示。

在试验区布设的37个角反射器控制点(其中28个测区控制点，9个定标点)中，21个(12个测区控制点，9个定标点)用于参与平差计算，16个作为检查点(图1)。

3.4 相高转换

用平差计算方法获得各景影像准确的外方位元素后，就可进行相高转换，以获得影像上各点的高程。由式(1)—(4)可知，获取像点的高程需已知该点的航高H、基线倾角α、基线长度B、绝对时间延迟τ和相位φ。其中，H，α，B和τ延迟由参数定标获得，相位由干涉处理数据及区域网平差参数共同确定。通过相高转换获得影像上各点的高程。

3.5 地理编码

为便于应用，需将原始斜距几何下的InSAR高度数据投影到某一参考系统下。本研究基于正侧视模型，利用载机飞行数据信息，对高度图像进行地理编码。由POS数据获取方位向上每扫描行数据的载机位置(XS，YS)，坐标转换得到其高斯坐标(xs，ys)。为了能更准确地确定目标点位置，将通过平差计算得到的多普勒中心频率fDC引入斜视角β[1]，即

由式(1)求出侧视角θ。设θy为航偏角(由区域网平差结果得到)，利用几何关系得到该像点在高斯平面上的坐标(xp，yp)[10]，即

通过以上过程便可获得每一像点的地理坐标，再进行重采样，将影像的高程h(i，j)赋到重采样后的像点上，得到地理编码的DEM。

3.6 影像镶嵌

本文将试验区分为4个条带进行数据获取。条带间重叠率≥40%;条带内相邻影像重叠率为20%左右。在生成各景DEM之后，需在统一的坐标系统(WGS-84坐标系高斯-克吕格投影3°分带，1985国家高程基准)镶嵌DEM。具体步骤如下:

1)对于非重叠区域的像素，其高度值由覆盖该像素的DEM直接赋与。

2)对于二度重叠区域的像素，其高度值等于覆盖该像素的2个高度值的平均值。

3)对于重叠度≥3区域的像素，采用加权平均的方法[3]确定，将该点上的所有高度值进行排序，把偏离中值超出阈值(设为σ=±1.2 m)的高度值从序列中剔除，剩余高度值按其干涉相干度值进行定权并进行加权平均，得到该点的高度值。

4 试验结果及精度分析

4.1 试验结果

基于本文提出的高精度DEM的生成技术流程，对试验区的1个架次、4条带及76景SLC影像进行地形制图处理，利用自主开发的机载InSAR地形制图处理系统软件快速实现上述处理过程，生成了覆盖超过500 km2的一整幅高精度DEM图像(图5)。

图5 处理生成的DEMFig.5 Generated DEM from test data

4.2 精度分析

为了评价上述流程所生成的DEM精度，本研究将精确布设在试验区但未参与平差的16个角反射器(图1中◆所示)作为检查点，将其平面坐标及高程的实测值与生成值进行比较，评定所生成DEM的内部精度(图6)。

利用在试验区实测的185个GPS检查点进行精度检核。由于这些检查点在SAR影像上的位置特征不明显，刺点将会产生较大的位置偏差，所以未统计其平面误差，仅检核高程。发现其中有4个GPS点落入到阴影区域中，在后续比较中排除这4个点。GPS检查点多位于开阔区域且其测量精度较高，根据GPS检查点的平面坐标将其与DEM数据叠加，用InSAR DEM高程值HDEM减去检查点高程值HGPS得到外部检核的高程误差值，如图7所示。

图6 X，Y方向的平面误差分布Fig.6 Error distribution of X and Y direction

图7 高程误差平面分布及其统计直方图Fig.7 Distribution and statistical histogram of elevation error between GPSs and InSAR DEM

表2列出了16个角反射器检查点(JC01—JC16)、181个地面 GPS高程检查点(GZ01—GZ181)对DEM的检测精度，并对检查点的平面点位中误差和高程中误差分别进行了统计分析。

表2 DEM精度统计Tab.2 Precision statistics of DEM (m)

可以看出，图6中X方向误差和Y方向误差在0值附近分布，其均值分别为0.15 m和0.01 m;在图7左侧，高程误差在正负区间均有分布，且大部分高程误差集中在±1 m之间，其均值为0.057 m，标准差为±0.547 m;在图7右侧，高程误差频率分布曲线略呈左偏态，而非正态分布。高程误差在正区间0.2处两侧分布的点明显偏多，即InSAR生成DEM略高于实测GPS检查点的高程。这主要是由于X波段微波的表面穿透性较弱，InSAR获取的并不是裸露地形表面的高度[2]。

通过精密布设的角反射器检查点的检核，DEM的平面点位中误差为±1.188 m，高程中误差为±0.508 m;动态GPS测量检查点的高程中误差也达到了±0.550 m。可见，应用机载InSAR数据生成的DEM的精度达到了1∶1万丘陵地区一级高程中误差(表3)的要求。

表3 1∶1 万 DEM 技术指标[12]Tab.3 Technical indicator of 1∶10 000 DEM[12] (m)

5 结论

本文提出了利用国产机载双天线InSAR数据生成大面积高精度DEM的技术流程。基于自主开发的机载InSAR地形测图软件处理了高分辨率X波段机载双天线SAR数据，获得了大面积高精度DEM，并对生成的DEM进行精度分析，得出如下结论:

1)对影响高程精度的基线长度、基线倾角、相位偏置等参数，利用布放合理、精密测量的角反射器作为地面控制点进行外定标处理，这对于保证制图精度非常重要。

2)区域网平差不仅可以减少地面控制点数量，还能提高连接点坐标的精度和区域网整体性，是解决面积大、控制点稀少条件下地形测图的关键技术。

3)对于丘陵地区，试验证明了基于国产机载双天线SAR数据的DEM生成技术可满足1∶1万比例尺的地形图制图精度要求;机载InSAR技术可作为复杂地区地形测图制取的一种技术手段。

志谢:感谢中国科学院电子学研究所提供了机载双天线SAR数据。

[1]孙中昶，郭华东，李新武.机载双天线InSAR数据生成高精度DEM 的误差分析[J].高技术通讯，2012，22(2):171-179.Sun Z C，Guo H D，Li X W.Error analysis of high precision DEM generated from airborne dual-antenna interferometric SAR data[J].Chinese High Technology Letters，2012，22(2):171-179.

[2]Li X P，Baker A B，Hutt T.Accuracy of airborne IFSAR mapping[C]//Proceedings of the American Society of Photogrammetry and Remote Sensing，XXII International Congress，Washington，USA.2002.

[3]Wimmer C，Siegmund R，Schwabisch M，et al.Generation of high precision DEMs of the Wadden Sea with airborne interferometric SAR[J].IEEE Transactions on Geoscience and Remote Sensing，2000，38(5):2234-2245.

[4]Mercer B.National and regional scale DEMs created from airborne InSAR[C]//ProcPIA，2007，36(3):W49A.

[5]黄国满，张继贤，赵 争，等.机载干涉SAR测绘制图应用系统研究[J].测绘学报，2008，37(3):277-279.Huang G M，Zhang J X，Zhao Z，et al.Research on airborne SAR interferometry mapping system[J].Acta Geodaetica et Cartographica Sinica，2008，37(3):277-279.

[6]Xiang M S，Wu Y R，Li S E，et al.Introduction on an experimental airborne InSAR system[C]//Proceedings of the 2005 IEEE International Geoscience and Remote Sensing Symposium，Seoul Korea:IEEE，2005(7):4809-4812.

[7]刘艳华，赵 争，黄国满.机载InSAR数据自动生成DEM技术及其在内蒙古丰镇地区的应用[J].地学前缘，2006，13(3):104-107.Liu Y H，Zhao Z，Huang G M.Automatic DEM generation technology by airborne InSAR data and its application to Fengzhen，Inner Mongolia[J].Earth Science Frontiers，2006，13(3):104-107.

[8]Sun Z C，Guo H D，Li X W，et al.DEM generation and error analysis using the first Chinese airborne dual-antenna interferometric SAR data[J].International Journal of Remote Sensing，2011，32(23):8485-8504.

[9]杨怀宁，郭华东，韩春明.机载InSAR敏感度方程定标限制条件的仿真实验[J].高技术通讯，2010，20(10):1049-1054.Yang H N，Guo H D，Han C M.Imitation experiment on the parameter requirements of the sensitivity equation-based calibration for airborne InSAR[J].Chinese High Technology Letters，2010，20(10):1049-1054.

[10]张 薇，向茂生，吴一戎.基于正侧视模型的机载双天线干涉SAR 外定标方法[J].遥感技术与应用，2008，23(3):346-350.Zhang W，Xiang M S，Wu Y R.Studies on outside calibration method based on the boresight model for dual-antenna airborne interferometric SAR[J].Remote Sensing Technology and Application，2008，23(3):346-350.

[11]马 婧，尤红建，胡东辉.FLeberl模型与干涉测量模型相结合的InSAR影像区域网平差[J].红外与毫米波学报，2012，31(3):271-276.Ma J，You H J，Hu D H.Block adjustment of InSAR images based on the combination of FLeberl and interferometric models[J].Journal of Infrared and Millimeter Waves，2012，31(3):271-276.

[12]国家测绘局测绘标准化研究所.CH/T1008—2001基础地理信息数字产品1∶10 000，1∶50 000 数字高程模型[S].北京:国家测绘局，2001.National Institute of Surveying and Mapping Standardization.CH/T1008—2001digital products of fundamental geographic information 1∶10 000，1∶50 000 digital elevation models[S].Beijing:State Bureau of Surveying and Mapping，2001.