灾后复杂地形区域的测量测绘模型设计*

林元茂，李 建，韩 立

(1.重庆工程职业技术学院，重庆 江津402260；2.成都理工大学 生态资源与景观研究所，四川 成都 610059；3.成都理工大学 国家环境保护水土污染协同控制与联合修复重点实验室，四川 成都 610059)

我国是世界上自然灾害种类最多的国家，自然灾害对地理状况会产生很大的影响，以2008年5月12日汶川地震为例，地震共造成69 227人死亡，374 643人受伤，17 923人失踪，严重破坏地区超过10万km2，其中，极重灾区共10个县(市)，较重灾区共41个县(市)，一般灾区共186个县(市)[1]。地震造成地震重灾区和极重灾区许多地方的地形地貌出现极大的变化，大量的建(构)筑物都遭受严重的破坏，同时也会给灾后地形的测量测绘工作造成极大的困难。因此，可靠的灾后复杂地形区域的测量测绘方法是目前急需解决的问题。许多学者针对这一情况进行了研究，王旸[2]提出基于重建三维地形的虚拟重建方法，可以获得相关的地形复杂图像，但该方法未对地形图像进行干涉解缠，大大降低了测量测绘的精确度。乐林等[3]提出一种估计地形复杂度的LiDAR点云多尺滤波方法，该方法需要对大量云数据进行多尺度滤波处理，复杂度高，测量效率低。采用基于单一线性模型的测量测绘方法测量灾区复杂地形过程中，不能滤波处理生成的DEM和外部DEM范围内的系统粗差，不能过滤掉无价值高程点，测量效果差。

为了解决上述方法中存在的问题，本文设计一种用于灾后复杂地形区域的测量测绘模型。通过遥感技术进行定标，进而获得灾区复杂地形的精准影像，然后运用线性回归方法得到高程点，最后通过实验验证本文方法的准确性、有效性。

1 灾后复杂地形区域的测量测绘模型设计

1.1 精化灾区地形DEM

灾后复杂地形区域遥感InSAR影像受到解缠相位残余、雷达阴影等造成的相位粗差干扰，会对最终灾区地形DEM(离散单元法)精度造成影响[4]。因此，需要精化灾区地形DEM，获取高精度的灾区地形高程图。在干涉相位误差的基础上，对灾区地形DEM进行精化，精化过程在外部DEM辅助下完成，对灾区遥感InSAR影像进行干涉定标，建立干涉影像和外部DEM的精确关系，将定标后精度高的外部DEM从地理坐标变换到雷达距离-多普勒坐标系下，实现干涉影像与外部DEM的精准对应，在解缠干涉影像图中通过多模型的线性回归逐个像素精确的去除相位趋势，最后将无价值高程点滤除，获取精确地灾后复杂地形区域高程点[5]。

本文模型主要是将初始斜矩、相位偏置、基线长度、基线倾角等参数通过现有地面定标点坐标求得的过程，现有地面定标点通常由布设角反射器所得[6]。通过干涉定标可以将灾区地形机载InSAR图像参数对DEM精度的影响消除，将基线量测精度提高，由此达到提高最终灾区地形DEM精度的目的。由现有的国内外研究结果可知，基线估计精度通常都很低，尤其是机载InSAR，因为其基线过短，细微的改变都将影响最终的DEM精度，因此基线定标的单位需要精确到毫米。

1.1.1 定标方法和流程

设定一三维坐标系，距离向方向用X轴表示，方位向方向用Y轴表示，垂直于X轴和Y轴组成的平面用Z轴表示。两个天线的位置坐标用A0和A1表示，A0和A1的间距就是基线长度，第一个天线中心在坐标系内位置坐标用A0(0,0,H)表示，第一个天线中心距离坐标原点的高度为H，由不同坐标原点决定H的大小，第二个天线中心在坐标系中的坐标用A1(x,y,h)表示，用o1表示A1点在XY平面的投影。

两个天线至同一地面控制点的斜矩用R1，R2表示。将待定标参数确定，基线长度用At表示，基线倾角用Aa表示，初始斜矩用P0表示，干涉相位用φ表示。则依据相高转换公式，可得：

F(At，Aα，P0，φ)=0。

(1)

求导上述公式，得到方程为：

W=DΔx-l。

(2)

式中：W是公式(1)的误差，公式(2)内每个控制点处对其待定标参数的导数用系数矩阵D表示。待定标参数的改正数用Δx表示，公式(2)的初值用l表示。依据最小二乘准则，可得：

WTW=min。

(3)

可以求得改正数：

引理2.1[34] 设(X;≤,→,→→,1)为伪BCI-代数，则下列结论成立：对任意x,y,z ∈ X,

Δx=(DTD)-1Dl。

(4)

将以上公式通过迭代处理至改正数小于目标阈值时，迭代结束。

定位时需要确定初始基线倾角、基线长度等初始定标参数，参数可以通过全站仪测量雷达天线罩各角点获取[7]；通过影像配准、干涉处理、相位解缠等处理干涉影像；通过雷达影像内的量测获取定标点的像点坐标；利用定标点的像点坐标，在解缠相位图中得到各定标点的解缠相位；通过最小二乘法对公式(4)进行迭代求解，获取新的系统参数改正数；将改正数与系统参数初始值相加，用处理过的系统参数实施相高转换；最后用检查点检测高程精度。

为了确保干涉影像同外部DEM间的高精度对应，需要对干涉定标后的高精度外部DEM进行变换处理，将DEM从地理坐标变换到雷达距离-多普勒坐标系中。

1.1.2 线性回归分析

灾区的地形趋势复杂多变，造成灾区复杂地形干扰影像中存在解缠干涉问题，使得灾区地形DEM的精度大大降低，因此将多模型的线性回归策略与地形趋势相拟合，通过多模型的线性回归方法，分析解缠干涉影像图中的相位趋势[8]。因为灾后地形复杂的原因，需要选取多个回归模型，将解缠干涉影像图中的干涉相位依据区域增长算法分成多个子影像块，但子影像块不超过10个，按照地形变化分割干涉相位。依据以下步骤，对所有影像块实施回归分析。

采用式(5)将雷达距离-多普勒坐标系中的外部DEM转换成干涉相位：

(5)

式中：外部DEM转化得到的干涉解缠相位用ηz表示；干涉解缠系数用λ表示；基线长度和斜矩方向上传感器同目标间距离用B和R表示；入射角同基线以及水平方向间的夹角为θ和α；外部DEM高程值为Q。

(6)

1.2 获取有价值高程点

解缠干涉图在相位趋势除掉后可以转化成有效灾区地形高程图，通过公式(7)计算得到高程差函数模型为：

(7)

InSAR获取的高程值和外部DEM的距离范围可以去除解缠干涉图的系统粗差[13]。用μ表示误差的均值，δ表示标准偏差，依据Chebyshev理论，在区间μ±4δ时，不论误差服从呈何种分布，误差的概率最小为94%。将高程差高于4倍外部DEM误差的无价值高程点过滤掉，求解InSAR与对应外部DEM得到高程点的绝对高程差[14]，挑选出绝对高程差高于4倍的外部DEM误差的高程点，重构剩下的高程点得到精确地灾后复杂地形区域DEM。

将无价值高程点滤除后，因为InSAR影像在一些区域相干性过低，导致在影像中形成连续无价值的高程点，这种无价值高程点被称为“空洞”。超过20个像素的连续区域称为影像空洞，其保留了无价值的高程点，利用双线性插值方法在较小的连续区域中得到高程点，很好地解决“空洞”问题[15]。

依据式(7)获取的有效高程图以及有价值高程点，可实现灾后复杂地形区域的准确测量测绘。

2 测绘模型性能分析

为了测试灾后复杂地形区域测量测绘模型的测绘精度、测绘效率等性能，进行了如下实验检验。

2.1 测绘对象

实验选取某灾后横断山脉测图区域进行验证，此区域属高原高山地形，起伏较大，地貌复杂，容易发生地震、山体滑坡、雷击、冰雹、塌方、暴雪等自然灾害。地势高、高差大是此测图区域的地貌特点。测区内最高海拔7 556 m，最低海拔540 m，海拔平均为2 000～6 000 m。该横断山脉整体地貌景观为东西并列，南北纵贯。在该区域中选取一地势较为平缓的试验区，该区域最高海拔高程为3 312 m，最低海拔高程为3 094 m，高差为218 m。

2.2 灾区地形DEM的精化及结果

将1：100 000的外部DEM引入到实验中，外部DEM将实验地区的地理坐标变换成雷达距离-多普勒坐标时，需要依据高程值设置的原始地理编码和轨道参数进行分析。在此基础上，采用本文模型对实验区域InSAR幅度影像实施模拟，得到精化灾区地形DEM后的SAR幅度影像模拟效果图如图1所示。本文模型模拟的灾区地形InSAR幅度影像同真实InSAR纹理基本一致，本文模型能够高精度模拟灾后复杂地形区域的InSAR影像，为该地区的准确测量提供准确的分析依据。

图1 SAR幅度影像模拟效果

2.3 灾区遥感InSAR影像的干涉定标结果

为了获取与干涉图在几何关系上相对应的外部DEM，可以通过模拟InSAR幅度影像与真实InSAR幅度影像一一对应精化地理编码表，由外部DEM转化获取解缠的干涉相位。采用解缠前后的干涉相位图如图2所示。

图2 解缠前后的干涉相位图

经过干涉处理后比例为1：100 000的外部DEM，高程精度优于26 m，设置4倍外部DEM精度约100为InSAR高程点的滤波阈值。把图2b解缠后的干涉相位图分成3个子影像块，可根据区域增长算法来分割。在3个子影像块内依次选取相干性大于0.5的像素点分析线性回归，进而清除相位趋势。将InSAR生成的灾区地形DEM与其对应的外部DEM相减后，过滤掉绝对高程差高于100 m的高程点。实验共有15 999 879个高程点，理论上应该有4 798 347个高程点被清除，由于影像内有“空洞”，所以只有4 701 324个高程点被滤除，剩下的点在InSAR生成的高程结果中被保留，从最终的图3中1:50 000的DEM结果中可以看出，精化后的灾区地形DEM比之前的DEM精确很多，说明采用本文研究的结果可提高灾后复杂地形区域绘图的精确度，增强区域测量测绘效果。

图3 1:50000DEM结果

2.4 不同模型的测量效率对比分析

把实验灾区地形InSAR影像分成不相同的数据块，将单一线性模型作为对比模型，以模型的测量效率为指标，进行对比实验，并得到两种不同模型的测量效率对比结果如图4所示。采用单一线性对比模型对实验灾区地形进行测绘，在0～45 s时间段内，测量的面积难以超过600 km2，当45～55 s时，测量的面积大幅度增加，最大测量面积约为1 000 km2，而后的5s测量面积依然不超过600 km2；采用本文研究结果模型对实验灾区地形进行测绘，仅在25～30 s时间段内，测量面积较小，不超过600 km2，而其余时间段中，测量面积始终保持在不低于1 000 km2的范围内。对比两种结果充分说明，在较短的时间内，本文模型即可达到较大的测量面积，且始终保持在1 000 km2左右的测量范围，采用该模型的测量效率更高。

图4 本文模型测量效率

2.5 不同模型的测量精度对比分析

测量点的高程是测试模型测量精度的一个重要指标，高程能够直接影响测量精度的高低。为此，将单一线性作为对比模型，进行测量精度的实验测试。在实验区内选取10个测量点，测量点在InSAR图像中均匀分布，本文模型和单一线性模型检测10个测量点得高程结果以及实际高程结果对比结果，分别用表1与表2表示。分析表1得出，对复杂灾区地形进行测绘时，在第9个测量点时，本文模型测得的高程为3 184.11 cm，实际高程为3 168.54 cm，计算高程误差为15.57 cm，是本次实验的最大误差；在第6个测量点时，本文模型测得的高程为3 167.21 cm，实际高程为3 164.52 cm，计算高程误差为2.69 cm，是本次实验的最小误差。10个测试点的误差区间为±15 cm。

分析表2得出，采用单一线性模型对复杂灾区地形进行测绘时，在第1个测量点时，单一线性模型测得的高程为3 265.96 cm，实际高程为3 284.64 cm，计算高程误差为-18.68 cm，是本次实验的最大误差；在第3个测量点时，单一线性模型测得的高程为3 288.44 cm，实际高程为3 283.41 cm，计算高程误差为5.03 cm，是本次实验的最小误差。10个测试点的误差区间为±19 cm。

表1 测量点高程及测量误差对比

表2 单一线性对比模型高程与测量点实际高程对比表

对比表1和表2的实验结果可以看出，单一线性模型的最大高程误差比本文模型的最大高程误差大3.11 cm，其最小高程误差比本文模型的最小高程误差大2.53 cm，10个测试点的误差区间要比本文模型10个测试点误差区间范围大出±4 cm。由此结果充分说明，本文模型的高程误差范围较小，相对的测量精度较大。

综合上述四组实验结果得出，本文研究得到的灾后复杂地形区域的测量测绘模型得到的SAR幅度影像模拟效果较好，解缠前后的干涉相位图十分清晰，保留的高程点较为准确，且经过实验对比分析，测得所设计模型的测量效率和测量精度均远远高于单一线性模型的测量效率和测量精度，具有更高的应用价值。

3 结论

针对灾后复杂地形区域，采用基于单一线性模型的测量测绘方法测量时，无法有效过滤掉无价值高程点、测量精度低的问题，本文研究得到高精度的灾后复杂地形区域的测量测绘模型，依据高分率雷达影像的外部DEM辅助实现测量测绘，并且将多模型的线性回归策略与地形趋势相拟合，通过多模型的线性回归逐个像素精确的去除相位解缠中多余相位，得到有效高程图，从中采集精确地灾后复杂地形区域高程点。实验验证了本文模型在灾后横断山脉这种地势高、起伏大地貌复杂的地形情况下，模拟的灾区地形InSAR幅度影像同真实InSAR纹理基本一致，为该地区的准确测量提供准确的分析依据；表明本文模型有着可靠的精确度与稳定度，能满足灾后复杂地形区域的测量测绘要求，具有较高的实用性。