升降轨PSInSAR观测反演沉降与水平向位移试验

王艳， 张玲， 葛大庆， 张学东， 李曼

(1.中国国土资源航空物探遥感中心,北京 100083; 2.北京建筑大学测绘与城市空间信息学院，北京 100044)

0 引言

雷达干涉测量(InSAR)所获取的形变监测结果为视线向形变量，可依据形变方向与视线向的投影关系，进行不同方向形变量的求解。上下、东西和南北3个方向上形变量求解的前提是具备多维(3个方向)观测值，但实际上，受制于雷达的成像方式，难以同时满足3个方向的独立观测。因而，直接利用InSAR监测结果反演三维形变量存在一定的难度。实际应用中，如地震等大尺度形变场，一般存在不同方向上的形变特征，可以利用升降轨、左右视观测[1-3]，以及图像偏移量跟踪方法等实现不同方向形变量的求解。

依据地表形变三维矢量的构成模型，升降轨InSAR在二维观测量的条件下难以直接反演三维变化，其根本原因在于InSAR升降轨观测量个数不够。同时，受制于InSAR观测对于形变的敏感程度和误差影响，南北方向求解的稳定性最差，其次为东西方向，垂直方向最为准确[4-6]。在3个不同方向的InSAR观测模式下具备了3组观测数据，但解析模型仍然不存在多余观测作为约束条件。尽管满足解方程的条件，其结果的可靠性仍需要其余约束。因而，当前对于三维形变量估算仍以观测量的构成模型，结合地表形变活动模型进行反演。根据雷达视线向实际观测值对各个方向的敏感程度，在PSInSAR观测条件下，对三维模型进行简化可求解水平和垂直方向上的移动量。

本研究中以升降轨观测下的形变量模型构成和相位观测在各个方向上对形变量的敏感程度，将其简化处理以求解上下和东西2个方向的移动量。试验中利用ENVISAT卫星获取的升降轨数据，对PSInSAR处理得到的形变速率进行沉降和水平向移动量反演试验，分别求解了水平移动和垂直变形速率，结果表明试验区内具有垂直移动占主导，水平移动极其微小的特点。

1 升降轨反演地表形变量的简化模型

融合升降轨InSAR观测的形变量反演方法可以区分垂直方向和水平方向(东西方向)形变移动量。考虑到升降轨2组观测是从2个几乎相反的方向获得数据(图1)，对于地表同一目标点而言，东西方向的水平移动在2组观测数据中的表现是基本相反的形变，而垂直方向上则表现为基本相同的形变。视线向形变量(dLOS)的构成[4]为

dLOS=(UNsinφ-UEcosφ)sinθ+UVcosθ。

(1)

图1 升降轨观测模式下的垂向与东西向形变分解示意图

式中:UN，UE和UV分别表示北南、东西和垂直方向的地表移动分量；φ为方位角(正北与传感器飞行方向夹角，自北顺时针为正)，按照习惯降轨方向设为正，升轨方向为负；θ为入射角。对于ENVISAT雷达卫星扫描模式下的IS2数据而言，升降轨条件下其方位角分别为-11.9°和191.9°，由此可近似得到sinφ≈0(为-0.19)，cosφ≈±1(升轨时为0.98，降轨时为-0.98)； 代入式(1)中，视线向形变量表达式可进一步简化为

dLOS≈UEsinθ+UVcosθ

，

(2)

式中右式第1项在降轨时取“+”号，升轨时取“-”号。由于入射角在感兴趣区的变化极其微小，因此可用其均值来估计地表形变的水平和垂直分量。对应于ENVISAT IS2数据，整景覆盖的入射角变化低于6°，一般的监测范围内均小于该值。由此，假定视线向产生1 cm的形变量，而由采用入射角均值所引起的最大误差低于1 mm。

显然，根据式(2)不能区分地表形变的南北方向的水平移动。实际上，南北方向的水平移动对垂向和东西向分量估计的影响非常微小，垂向分量估计几乎不受南北向移动的影响。假设目标在南北方向产生1 cm的水平移动，将其错误估计为垂向移动的误差为0.8 mm，而在东西向的分量错误估计为2 mm。

因而，由式(2)可推出升降轨观测模式下东西和垂直方向的形变量求解公式为

,

(3)

。

(4)

式中:dLOS(desc)表示降轨数据的视线向形变量；dLOS(asc)表示升轨数据的视线向形变量。

利用式(3)(4)分别解算东西和垂直方向的地表移动分量的前提是要求升降轨模式下同一观测区域内具有相同的观测对象，即分别在升轨和降轨解算数据中具有同名像元。这在相干性良好的差分干涉图解算地震形变场移动分量中较为常见。而对于非线性形变场中相干目标的解算，该方法的缺点是，在地形起伏较大的地区，相对于单个数据集可用的像元数目，用来估计东西向和垂向的地表移动分量的像元数目较少。这是由每个像元在升轨和降轨数据集中固有的成像方式及其对相干性的影响所引起的。由于单个目标的解算依赖于对同名相干像元的识别，而实际上升降轨下完全同名点的数目有限，对于估计整体形变场不利。为此,针对升降轨PSInSAR观测值，可将升降轨中的一景作为主影像，另一景为辅影像，以此来插值提取相同位置的变形量，从整体上反演地表形变的水平和垂直位移。

2 PSInSAR数据处理关键步骤

永久散射体干涉测量(PSInSAR)技术的核心是对相干像元，即永久散射体(permanent scatter)的差分干涉相位序列进行分析，根据差分干涉相位各分量的时空特征，估算大气波动影响、DEM误差以及噪声等，将其从差分干涉相位中逐个分离，最终获取每个PS的线性和非线性形变速率以及DEM误差等参数。

干涉像对组合是PSInSAR处理的重要过程，其目的是将干涉像对进行优化组合。由于地形相位误差和相干性均受空间基线的影响，而短基线条件有利于相干性的保持，因而选择空间基线小于给定阈值，如ENVISAT可选300 m的像对生成干涉图。同时，顾及到地表形变速率的大小，对于一般地面沉降监测而言，选择时间间隔小于2 a，大于2～3个卫星重复周期的像对生成干涉纹图，在多数情况下是适用的。

对识别出的相干目标候选点的差分相位序列的互差进行处理是PSInSAR处理的关键所在，其过程是以Delaunay三角网构建不规则格网连接将其进行连接，对邻近点pm和pn的干涉相位时间序列进行二次求差，应用二维周期图估计点间形变速率和高程误差。式(5)为利用点间相位差估计形变速率和高程误差的函数模型，在完成频率估计后，利用最小二乘法实现整个三角网的积分，生成形变速率图，以反映形变场分布特征。

(5)

其中

(6)

式中：γmodel为模型相关系数；pm和pn分别为2景影像同名点对应的第m和第n个点；N为干涉像对个数；φdiff表示差分相位；Ti为第i景干涉纹图的时间基线；φmodel表示模型相位；λ为雷达波波长；γi为第i景干涉纹图的相关系数；bi为第i景干涉纹图的垂直基线；θ为雷达入射角；υ(pm)和υ(pn)分别表示第m点和第n点的形变速度；ε(pm)和ε(pn)分别为第m点和第n点的DEM误差。

3 升降轨PSInSAR形变量反演与分析

3.1 升降轨下PSInSAR观测量

研究需选择同时获取了升轨和降轨模式下ENVISAT-ASAR数据的地区(如图2中红色框和粉色框的重叠区)。本研究选定的图像(图2中的蓝色框区域)覆盖范围为30 km×30 km，位于苏州市区西北周边，分别获取了2006—2010年间升轨(Track-39)和降轨(Track-275)雷达数据24景和27景。

图2 研究区及升降轨数据分布

PSInSAR处理得到的形变速率图为离散点模式，而难以直接提取同名点，因而需对2幅图进行插值处理，生成连续分布速率图。图3所示为升降轨模式下Track-39和Track-275的PSInSAR处理结果。在插值处理之前，先进行坐标系的统一，以消除

图3 升(左)降(右)轨模式下PSInSAR视线向形变速率(空间插值后)

由于解算过程中参考位置不同而产生的基准偏差[7-8]。坐标系的统一采用地形校正方法，然后在相同的地面坐标系下同时完成2幅SAR影像的地理编码。地理编码的误差均小于0.15个像元，2图像之间不再进行二次匹配。参考基准偏差的补偿利用同名点统计整体偏差值，进而补偿形变场整体差异。最终得到的2景图像中各形变中心分布位置相同，监测结果量值相近(图3)。

3.2 地表变形量求解

根据式(2)(3)对入射角进行归一化(以正弦为例)处理，按照同一模式下的入射角进行分析。图4和图5分别为研究区入射角归一化结果和方位角归一化结果。

图4 升(左)降(右)轨模式下归一化入射角

图5 升(左)降(右)轨模式下归一化方位角

在升降轨解算时，对于入射角的考虑为2组数据集里每个像元的入射角，而实际上这种情况下难以进行参数估计。因此，结合本研究区范围较小的实际，对每个像元均按照统一入射角来计算。图4中入射角的最大差别为0.05 rad，图5中方位角的最大差别为0.002 rad，表明二者的变化极其微小，对入射角变化的简化和对方位角变化的忽略处理是正确的。

利用升降轨视线向观测量解算垂向与东西向形变速率(图6)。解算结果表明，研究区垂向变形占主导，即以地面沉降为主，而水平方向变化极其微小。主要沉降变形区(图中A区)水平变形较为显著，这主要是在地势平坦地区由于不均匀沉降引发地表曲率变化，进而产生向沉降中心方向的轻微滑动。对照图6中左图所示，水平变形相对显著的地区与沉降中心位置相同。这种形变类型在矿山沉陷

等大变形中表现较为普遍; 在地形起伏较显著的地区(图6B区)，其水平变形量相对较明显，这主要与地形有关，地形影响或主导了该地区发生形变的活动方向，这类变形在滑坡等单一方向为主导的移动中较为普遍，而在地势平坦地区不显著。

4 结论

试验研究表明，利用简化的二维模型可实现对升降轨模式下InSAR观测量的分解，从而可利用升降轨数据求解地表垂直方向和东西方向的变形。由于模型过于简化，适用性会受到很大的局限。对于文中研究区地势起伏变化小的平原地区，升降轨模式对于沉降监测影响均不明显，其监测结果都能较好地反映地表形变的实际情况，受雷达观测方向影响不明显。这为大范围地面沉降监测提供了依据，即利用统一模式观测可满足沉降监测需要。

实际上，升降轨观测仅提供了二维观测量，仍然不能满足多参数反演的需要。受制于地形和地表形变类型的影响，结合高分辨率雷达观测数据，针对形变模型的参数反演将更为有效。这也是今后需要进一步研究的方向。

参考文献(References)：

[1] Fialko Y,Sandwell D,Simons M,et al.Three-dimensional deformation caused by the Bam,Iran,earthquake and the origin of shallow slip deficit[J].Nature,2005,435(7040):295-299.

[2] Funning G J,Parsons B,Wright T J,et al.Surface displacements and source parameters of the 2003 Bam(Tran)earthquake from Envisat advanced synthetic aperture radar imagery[J].Journal of Geophysical Research,2005,110(B9):B09406.

[3] Bonforte A F,Guglielmino M,Coltelli A,et al.Structural assessment of mount etna volcano from permanent scatterers analysis[J].Geochem Geophys Geosyst,12(2),doi:10.1029/2010GC003213.

[4] Hanssen R F.Radar Interferometry-data Interpretation and Error Analysis[M].New York:Kluwer Academic Publishers,2002:12-65.

[5] 夏耶.巴姆地震地表形变的差分雷达干涉测量[J].地震学报,2005,27(4):423-430.

Xia Y.Bam earthquake:Surface deformation measurement using Radar interferometry[J].Acta Seismologica Sinica,2005,27(4):423-430.

[6] 孙建宝,梁芳,徐锡伟,等.升降轨道ASAR雷达干涉揭示的巴姆地震(Mw6.5)3D同震形变场[J].遥感学报,2006,10(4):489-496.

Sun J B,Liang F,Xu X W,et al.3D co-seismic deformation field of the Bam earthquake(Mw6.5) from ascending and descending pass ASAR Radar interferometry[J].Journal of Remote Sensing,2006,10(4):489-496.

[7] 葛大庆,王艳,张玲,等.德州地面沉降-回弹及地下水位波动的InSAR长时序监测——以德州市为例[J].国土资源遥感,2014,26(1):103-109.

Ge D Q,Wang Y,Zhang L,et al.Seasonal subsidence-rebound and ground water level changes monitoring by using coherent target InSAR technique: A case study in Dezhou,Shandong[J].Remote Sensing for Land and Resources.2014,26(1):103-109.

[8] 王艳,葛大庆,张玲,等.升降轨PSInSAR地面沉降监测互检验与时序融合[J].国土资源遥感,2014,26(4):125-130.

Wang Y,Ge D Q,Zhang L,et al.Inter-comparison and time series fusion of ascending and descending PSInSAR data for land subsidence monitoring[J].Remote Sensing for Land and Resources,2014,26(4):125-130.