粗DEM辅助的广域PS-InSAR数据预处理方法及应用

沈 龙, 耿 浩, 马 仪, 程志万, 黄双得, 朱 宇, 刘 靖, 张钰松, 索志勇

(1. 电力遥感技术联合实验室(云南电网有限责任公司电力科学研究院)， 云南昆明 650217；2. 北京空间飞行器总体设计部， 北京 100094;3. 西安电子科技大学雷达信号处理国家重点实验室， 陕西西安 710071)

0 引言

星载合成孔径雷达(Synthetic Aperture Radar, SAR)作为一种全天时、全天候的微波遥感成像技术，为高精度、低成本地实现大范围观测以及地表形变监测提供了可能。干涉合成孔径雷达(Interferometric SAR, InSAR)利用具有相干性的两幅SAR图像可获得高精度的数字高程图(Digital Elevation Model, DEM)。永久散射体(Permanent Scatter，PS)干涉技术作为一种地表形变探测方法,克服了传统干涉/差分干涉合成孔径雷达技术中时间去相关、空间去相关和大气延迟等问题。

针对PS-InSAR的研究目前主要存在两个问题：第一是现有算法的研究基本都集中在对小观测区域的研究，观测场景数据大范围覆盖特性考虑不够充分；第二是算法假设数据误差都得到了有效校正。例如，肖行诠等人基于18景TerraSAR-X条带模式高分辨率影像数据，采用时序InSAR技术获取了成都地区尖山和九江2座500 kV变电站及周边电力铁塔的毫米级形变监测结果，观测区域约为几平方千米。朱珺等人采用16景Sentinel-1A数据结合PS-InSAR研究了风电场所在区域在2015年5月到12月期间产生的地表形变。随着形变测量手段的多样化和场景的复杂化，给干涉测量技术提出了新要求。2014年张学东等人利用2008—2010年邻轨ENVISAT ASAR数据进行轨道拼接，获取京沪公路沿线沉降，使跨轨道、多幅影像的大范围PS-InSAR监测成为现实。2016年Sousa等人利用2014—2016年C波段Sentinel-1数据对大坝的变形分量进行提取，显示了C波段传感器在大坝监测的特殊情况下的潜力。2019年 Fryksten等人利用永久散射体技术和结合获取了瑞典乌普萨拉市城区地表沉降，并结合精密水准测量结果验证了PSI方法的有效性。2020年张玲等人利用RADARSAT-2和TerraSAR-X两种不同波段，不同分辨率卫星数据开展唐山市城区主要活动断裂两侧微小差异性形变探测研究，并分析了不同波段卫星对同一地表形变监测的差异性和一致性。

传统PS-InSAR处理算法的处理性能依赖于精密科学轨道数据，但在实际应用中，由于精密科学轨道获取的时效性较差，而采用实时轨道会导致差分干涉相位中存在由于地形起伏导致的残余相位误差；并且超大场景处理时，会存在许多不同于小场景处理的特点与难题。因而，本文针对超大观测场景的PS-InSAR处理，提出了超大观测场景数据预处理的方法，可有效校正由于实时轨道精度较差引起的处理误差，并针对超大场景提出了能够适应配准误差空变性的序列图像配准方法，并结合超大场景处理时多景图像参考基准不一致的问题，利用一种类“区域增长法”的思路，完成超大观测场景的数据拼接问题。

本文第一节分析了轨道误差对观测场景地形相位去除(即差分干涉相位)的影响，并进一步分析了超大观测场景配准时存在的配准误差空变问题;第二节根据第一节的分析提出了基于辅助DEM的系统误差校正方法，并进一步提出了广域图像层级配准方法；第三节给出了改进的PS-InSAR处理流程；第四节利用覆盖云南全境的多景Sentinel-1A数据进行了超大观测场景PS-InSAR的处理实验，验证了本文方法的有效性;第五节总结全文。

1 广域差分干涉处理误差特性分析

1.1 轨道误差对差分干涉相位的影响分析

PS-InSAR处理基于差分干涉SAR处理，而差分干涉处理必须消除由DEM引起的地形相位。然而在实际系统中，由于轨道误差(尤其是只有实时轨道参数的时候)的存在导致SAR图像或DEM的定位产生较大误差，因而会严重影响差分干涉相位的精度，进而影响时间序列处理的稳健性。假定目标真实的坐标点为(,)，其对应的目标高度为(,)，由于轨道误差导致的方位向偏移量为和距离向偏移量为，其对应的目标高度为(+,+)，由于位置偏差引入的高程误差为=(,)-(+,+)，则由高程误差引入的差分干涉相位误差为

(1)

式中，为波长，为雷达斜距，为雷达下视角，为有效垂直基线。由式(1)可知：当高程误差达到1个模糊高度时，引入的差分干涉相位残余误差可达2π。因而需要校正由于轨道误差定位偏差导致的高程误差引入的差分干涉相位误差。

1.2 广域观测图像配准误差空变性分析

图1给出了某一方位时刻InSAR的场景观测示意图。假设轨道完全平行，在这种情况下，整个场景的基线长度和斜视角是固定的，方位向和距离向偏移量不发生耦合。我们认为地球表面是平坦的，卫星运动是平稳的。

图1 雷达干涉广域观测几何对应的配准偏移量空变示意图

根据雷达观测几何，点的主辅图像斜距差：

Δ=-=sin(-)

(2)

和分别为两天线到点的斜距。′点的主辅图像斜距差为

Δ='-'=sin(+Δ-)

(3)

式中，'和'分别为两天线到点的斜距，Δ为′点相对于点的下视角变化，则′的斜距差之差为

Δ=Δ-Δ=cos(-)Δ

(4)

由上式可知，如果将点的数据进行了配准，则′点的数据由于雷达下视角不同，其斜距差会增大，Δ变化较大时，配准偏移量甚至会达到多个距离单元，如下式所示：

Δ=Δ

(5)

式中为距离向采样单元。因而采用整体配准方法对广域图像进行配准已不再适合，必须提出新的方法解决由于广域观测面临的配准偏移量空变问题。

2 DEM辅助的PS-InSAR广域形变探测数据预处理方法

通过上一小节分析可知：必须对存在PS-InSAR使用的数据进行预处理，校正由于轨道误差导致定位错误而引入的DEM误差，并且克服由于广域观测时雷达几何差异导致的配准误差空变。本小节详细讨论对广域形变探测时的数据预处理方法。

2.1 基于辅助DEM的系统误差校正方法

高精度的差分干涉处理结果依赖于星载精密科学轨道数据，然而对于地表形变、滑坡等应急地质灾害，所使用的数据往往是卫星应急拍摄的，因而只具有低精度的实时轨道数据，利用实时轨道数据进行差分干涉测量或者PS-InSAR处理，轨道误差和系统定时误差(也即斜距误差)会严重影响处理结果的精度，因而需要对轨道误差和系统定时误差进行校正。本文提出一种基于低精度辅助DEM数据的系统误差校正方法，其处理流程如下：

第一步：观测场景对应DEM的获取

根据卫星星历参数中观测场景的四角经纬度利用定位方程获取地面观测场景对应区域的DEM，通常情况下，选择的DEM数据范围略大于四角经纬度确定的范围，以确保能够获得足够的信息来估计系统误差。定位方程采用距离多普勒定位方法：

(6)

第二步：基于DEM信息的目标定位

1)在待定位的SAR图像中选择二维格网，然后利用迭代定位方法获取格网像素的高程和经纬度，输入参数包括外部地表数字高程库、卫星状态矢量和格网像素在SAR图像中的坐标等；

2)对格网像素的经纬度进行插值，得到非格网像素的经纬度；

3)根据非格网像素的经纬度插值DEM，取出相应的高程。

第三步：基于入射角模型的模拟图像获取

根据Ulaby模型，由于地形信息引起的不同入射角对应的归一化后向散射系数：

=1-sin

(7)

第四步：基于模拟图像的系统误差计算与校正

模拟SAR图像基于DEM和实时轨道卫星星历参数获得，而实测SAR图像是由真实的数据录取几何进行成像得到，因而将模拟SAR图像和实测SAR图像利用传统方法进行精确配准，获得二维偏移量，然后根据二维偏移量即可获得实时轨道系统误差校正。假定利用配准获得的方位向偏移量为Δ和距离向偏移量为Δ，则方位向延时误差Δ和距离向延时误差Δ为

(8)

式中，为脉冲重复频率，为雷达距离向采样频率。

第五步：实时轨道系统误差修正

(9)

式中，和分别为实时轨道星历参数记录的方位起始时间和斜距时间。为清楚地表明实时轨道数据对差分干涉SAR的影响，对轨道误差校正前后的差分干涉相位图对比如图2所示。理想情况下，干涉相位在去除地形相位后的剩余相位只有噪声，存在误差时会有残余相位存在。从图2(a)和图2(b)的对比结果可知：经过系统误差校正后，差分干涉相位图中的相位条纹变化减少，说明经过系统误差校正后对地形相位的获取更加精确，也说明了系统误差校正方法的有效性。

(a) 系统误差校正前的差分干涉相位

2.2 基于辅助DEM多级图像精确配准方法

对于小区域的观测场景，由于雷达观测几何差异引起的主辅图像配准偏移量通常可以用统一的映射模型来表示，但是随着观测场景的增大，尤其是距离向观测带的增大，使得主辅图像之间存在较大的几何畸变，此时再用统一的映射模型对主辅图像进行配准会引起较大的失配误差。本小节提出一种层级图像配准方法以解决广域观测图像的精确配准，为后续PS-InSAR处理提供高相干的PS-InSAR时间序列图像。方法的具体处理步骤如下：

第一步：图像整体配准

为避免传统配准方法大范围搜索带来运算复杂度增加的问题，本文采用基于雷达观测几何的图像配准方法先对雷达图像进行整体移位，使得InSAR系统的主辅图像能够局部配准。

第二步：像素级图像配准

由于距离向配准偏移量在近距端和远距端是有差异的，所以需要根据粗配准要求(即配准误差达到像素级)确定距离向分块大小，采用的粗配准方法可利用传统的图像配准方法。距离向分块原则如下：在距离向的配准误差偏移量不超过一个分辨单元。

第三步：亚像素级图像配准

广域图像配准时必须采用分块处理，但是分块后的数据存在拟合参数不统一的问题，会导致配准后的图像相干性存在“马赛克”现象，进而会引起处理流程中后续步骤的误差传递，最终影响PS点识别数量及形变速度估计精度损失。本文在传统配准方法的基础上，提出采用基于稀疏控制点三次样条内插的方法获取广域图像的配准偏移量，直接对控制点的配准偏移量进行内插，然后根据配准偏移量直接计算辅图像在主图像中的映射关系，避免整体二次映射带来的配准精度下降。具体处理步骤如下：1) 将粗配准的图像分块；2) 在粗配准的图像上等间隔选择配准控制点；3) 利用最大频谱法对控制点的距离和方位偏移量进行估计；4) 利用图像频谱特性(或相干性)对图像的二维偏移量进行质量评估，剔除低于某一阈值的控制点；5) 对高质量配准控制点的二维偏移量进行三次样条内插，获得广域观测辅图像中每一像素在主图像中的对应位置；6) 对辅图像进行插值，获得配准后的广域干涉图像对。

图3(a)～(c)给出了层级配准不同阶段对应的干涉相位图，随着配准精度的提升，干涉相位图中干涉相位条纹越清晰，干涉相位质量也越高；图3(d)～(f)给出了层级配准时不同阶段对应的相干系数直方图，相干系数均值分别为0.647 2，0.713 1，0.866 4，从直方图分布及相干系数均值可知，层级图像配准方法有效提升了干涉图像对的相干系数，降低了相位噪声，有利于后续的时间序列散射点稳定提取和形变速度估计；图3(g)给出了距离向配准偏移量的空变特性，可以看到配准偏移量在距离向观测带宽内存在近似线性的特征，如果不考虑该空变特性，会导致干涉相位/差分干涉相位生成性能下降，甚至出现无法生成的现象。

图3 基于层级图像配准方法的结果对比

3 基于数据预处理的广域PS-InSAR形变探测流程

根据前述分析，本节给出完整的粗DEM辅助数据预处理的广域PS-InSAR形变探测流程。

3.1 基于数据预处理的PS-InSAR广域形变探测流程图

改进的全流程广域形变探测处理如图4所示，图4中所示灰色部分为广域形变探测中的改进部分。从图4可知：针对传统PS-InSAR处理，本文在广域探测时对时序处理数据进行了预处理，降低或消除了轨道系统误差和空变配准误差对时序差分干涉相位的影响，仿真结果也初步验证了方法的有效性。

3.2 广域图像拼接

Sentinel-1A一次观测获取的图像在配准时已经完成了图像Burst和子带之间的图像拼接，但是对于超广域(例如省域，全国的图像)结果，需要对多景数据进行拼接，然而多景图像的拼接不是按

图4 基于数据预处理的广域形变探测处理流程

照地理坐标进行简单排列，需要考虑多景图像之间由于轨道误差、滤波误差、解缠绕误差等引入的形变估计结果不一致现象。因此本文提出一种利用重叠区域像素估计结果进行联合加权处理的方法，并且利用“区域增长法”的思想，从平均相干系数较高的区域向平均相干系数较低的区域进行“增长”，以减小拼接误差的传递。广域观测图像多景数据联合拼接的代价函数构造如下：

(10)

式中，,,,分别为第个子图的相干系数均值和形变速度值，,,,分别为第个子图的相干系数均值和形变速度值，以此类推，完成超广域观测场景全域数据拼接。

4 基于Sentinel-1A数据的广域形变探测实验

本节利用Sentinel-1A哨兵IW模式SAR影像，采样率为13.9 m(方位向)×2.3 m(距离向)。选取的地表数字高程(DEM)数据为SRTM数据，网格间距为30 m×30 m。

4.1 云南昭通地区的形变探测实验

如图5(a)所示，昭通地势西南高、东北低，属典型的山地构造地形，山高谷深。汛期输电线周边极易发生地质灾害。将SAR强度图进行地理编码加载到Google Earth上的结果如图5(b)所示。昭通地区卫星数据监测时间为2018年3月—2019年5月内的33幅时序SAR图像。本文选取的时序影像数量已经满足PS-InSAR的数据处理量要求，根据时序图像的时空基线分布情况如图5(c)所示，选取2018年12月7日的SAR影像为主图像，将剩余32幅SAR辅影像与主影像分别按照本文所提的层级配准方法进行配准，获得配准后的SAR影像时间序列。估计得到的形变分布图如图5(d)所示。图中所述形变指的是雷达视线向方向(Line of Sight，LOS)的形变，正值代表朝向雷达运动，负值为背离雷达运动。为便于表达分别表述为抬升、沉降，绿色表示稳定，红色(沉降)和蓝色(抬升)的颜色越深，表示形变越大。图5(d)表明昭通地区的SAR影像(但是需要3个子带拼接)对应的形变处理结果在拼接后没有发生马赛克现象，表明文中所提方法的有效性。

图5 昭通地区时序数据时空基线分布及沉降探测结果

4.2 云南全境超广域形变探测实验

云南省全境界于北纬21°8′～29°15′，东经97°31′～106°11′之间，总面积39.4万平方千米，一次完全覆盖需要51景图像。本文实验数据为2017年5月至2019年5月的Sentinel-1A卫星升降轨获取的2 000余景。由于不能在一景影像内将全省地区覆盖，且每景影像的形变探测起始点不同，因而全省PS-InSAR形变探测结果还需要对多景数据进行拼接才能得到广域形变监测，图6(a)给出了云南全景所需的哨兵数据覆盖图简要示意，而图6(b)给出了云南全境的形变探测结果，形变结果呈缓慢变化趋势，跟图6(a)图像覆盖区域边界对比的实验结果表明，基于“区域增长”的加权拼接方式没有出现由于单景图像分别处理拼接后边界处所存在的“马赛克”现象。

(a) 图像覆盖示意图

4.3 外场试验比对

基于本文完成的云南省全境形变数据，并叠加输电杆塔GIS信息后，筛选出可能受地质风险隐患影响的输电杆塔，抽取其中部分输电杆塔进行了现场勘察。我们选取了某线路#84号杆塔进行现场核查，#84号杆塔及其周边区域分析如图7所示。

(a) #84号杆塔周边区域形变图

图7(a)～(b)为某杆塔及其周边区域的形变图和累积形变量，图7(c)～(d)为杆塔现场核查照片，经过现场对该杆塔4个塔基立柱顶端测量，发现杆塔基础有不同程度沉降。初步分析，造成#84号杆塔沉降的原因为：塔基位置土质为砂质土，孔隙比大，土质相对较软，强度较低，地基均匀性稍差，另外该塔基础为大板基础，埋深为3.4 m，较之掏挖基础稳定性低，加之基础在外部水分浸透等外部因素干扰，出现了不均匀沉降，引起塔材变形，塔身扭曲。现场考察测得的形变量约为2 cm，而处理得到的形变量约为1.5 cm，误差产生的原因主要是PS-InSAR长时间观测时间去相干导致相位精度下降导致的。

5 结束语

针对传统小区域PS-InSAR形变探测方法存在的问题，提出了一种广域形变探测时的数据预处理方法。首先提出了一种基于辅助DEM数据的实时轨道误差校正方法，能有效消除由于轨道误差和定时误差对差分形变相位的影响；并针对广域图像配准偏移量存在空变性的特点，提出了一种层级图像配准方法，解决了广域图像配准偏移量空变问题；并进一步提出了多景影像处理结果加权拼接处理思路，实现了多景图像超广域形变探测。利用仿真实验和云南全境Sentinel-1A覆盖数据的处理结果表明：基于粗DEM辅助的PS-InSAR预处理方法能够实现广域形变探测，对某输电线路杆塔所在位置形变结果的实地考证，验证了本文方法的有效性。