InSAR三维同震地表形变监测——窗口优化的SM-VCE算法

刘计洪，胡 俊，李志伟，朱建军

中南大学地球科学与信息物理学院，湖南 长沙 410083

干涉合成孔径雷达(interferometric synthetic aperture radar,InSAR)技术具有高空间分辨率，高精度等优点，已被广泛地应用于地形测量[1-3]和地表形变监测领域[4-9]。然而，传统差分InSAR(differential InSAR,DInSAR)技术仅可以获得真实地表形变沿视线向(line of sight,LOS)的一维形变投影，在地质灾害解译过程中具有极大的局限性[10-11]。为了克服这一缺点，文献[12—14]分别提出了像素偏移量跟踪(pixel offset-tracking,POT)、多孔径雷达(multi-aperture InSAR，MAI)技术和哨兵数据方位向子带重叠区域干涉(burst overlap InSAR，BOI)技术来获取地表真实形变沿雷达卫星方位向(azimuth,AZI)的投影形变。结合升降轨DInSAR和MAI/POT技术，利用加权最小二乘(weighted least square,WLS)即可解算地震、火山、冰川、滑坡等导致的真实三维地表形变场[5,15-24]。此外，全球导航卫星系统(global navigation satellite system，GNSS)观测数据与InSAR数据在空间分辨率和测量精度方面优势互补，二者融合也可实现高精度、高空间分辨率三维地表形变场的获取[25-29]。

然而，上述研究均基于单个点进行三维形变解算，并未考虑临近点三维形变之间的力学关系。文献[30]基于弹性理论提出了一种SISTEM(simultaneous and integrated strain tensor estima-tion from geodetic and satellite deformation measure- ments)方法，首次将应力应变模型(strain model,SM)引入InSAR与GNSS数据相结合的地表三维形变求解当中。但是，现实情况往往很难获取大范围的GNSS观测，文献[31]对SISTEM方法进行了扩展，仅需InSAR观测值即可进行高精度三维形变解算。为了方便起见，本文将DInSAR、MAI、POT和BOI等技术获取的形变观测数据统称为InSAR观测值。然而，在解算真实三维地表形变时，必须要融合不同方向的多源异质观测数据，因此精确确定不同数据之间的权重比例具有重要意义[6,32]。为了实现此目的，科研人员研究利用一个固定的滑动窗口求解窗口内数据的标准差作为先验信息进行定权[17]，这种定权方式是基于数据各态历经性的假设，即假设观测数据在某一数值上下一定范围内是随机分布的，但是对于地震等导致的大梯度形变而言，往往难以满足。也有相关学者在假设远场区域不包含形变信号的前提下，将远场一定区域数据的均方根值作为该类观测值整个区域的先验方差[33-34]。这种方法忽略了InSAR观测噪声在空间上的差异性，因此也会影响三维地表形变的求解精度。文献[25]在求解时序数据的平均形变速率时引入了方差分量估计(variance component estimation,VCE)进行定权，验证了VCE在多源InSAR观测数据定权方面的优越性。但是，对于瞬时形变场的求解，由于缺少多余观测，应用VCE进行定权受到了极大限制。

基于此，文献[35]提出了一种基于地表应力应变模型和方差分量估计的InSAR三维地表形变监测方法(SM-VCE)。该方法基于一定窗口范围内不同像素点三维形变之间的力学关系(SM)建立函数模型，显著增加了冗余观测个数，进而可利用VCE实现不同类InSAR观测值之间的精确定权。研究表明，基于窗口的SM-VCE方法比传统单点解算的加权最小二乘方法得到的三维地表形变精度更高[36]。在原始SM-VCE方法中，求解不同目标点的三维形变时，窗口大小是固定不变的，并且窗口内的InSAR观测值被认为是等精度的。对于同震形变场而言，在地震破裂带附近往往难以获得有效的InSAR观测值，使得固定窗口大小的SM-VCE方法难以获取完整的三维同震地表形变。即使在断层附近有InSAR观测值的情况下，窗口内往往会包含断层两侧的异质InSAR观测值，使得原始SM-VCE方法难以得到可靠的近场三维形变。而近场地表形变结果对于约束断层浅部滑动又是十分重要的[23,37-38]。因此，有必要重建地震破裂带附近三维地表形变场。同时，InSAR观测值极易受到失相干等噪声影响，窗口内同一类观测值中不同像素点的形变测量值精度往往各不相同，因此，忽略窗口内InSAR观测值精度差异的做法将会降低三维地表形变精度。

鉴于此，本文在SM-VCE方法的框架下，提出一种顾及断层破裂带及自适应变化窗口的观测值选取策略，进而基于邻近的InSAR有效观测值即可重建地震破裂带附近的三维地表形变场，同时，引入迭代加权最小二乘(iterative weighted least squares，IWLS)方法[39]实现窗口内同一类InSAR观测值的内部自适应定权。本文首先利用模拟试验验证了本文中窗口优化SM-VCE方法的可靠性和优势，进而基于升降轨哨兵1号卫星数据，获取了2019年7.1级Ridgecrest地震的高精度三维地表形变场。

1 SM-VCE方法

文献[35]较为系统地介绍了SM-VCE方法，本文在此基础上，提出了顾及断裂线和自适应变化窗口的观测值选取策略和基于IWLS方法实现窗口内同一类观测值内部的自适应定权方法，旨在基于SM-VCE方法的框架，获取精度更高、更为完整的三维同震地表形变场。

1.1 基于SM建立观测方程

dk=H·Δk+d0

(1)

(2)

文献[30]将矩阵H表示为一个对称矩阵和一个非对称矩阵之和

H=S+R

(3)

(4)

(5)

式中，ξ和ω代表地表应力应变模型中的应变参数和旋转参数。在原始SM-VCE方法中，H用式(3)表示，进而可将式(1)写成

(6)

(7)

(8)

式中

(9)

综合式(6)、式(8)，可得

(10)

式中

若假设点P0周围有Kj个第j类InSAR观测值，则最终可得

L=B·l

(11)

式中

此时，建立了综合考虑InSAR观测值与三维形变之间的几何关系，以及地表邻近点三维形变之间的力学关系的InSAR三维形变估计函数模型。

1.2 利用VCE实现不同类观测值之间的精确定权

(12)

式中，Wj代表第j类InSAR观测的初始权矩阵，一般取值为单位矩阵。进而根据方差分量估计算法[41]可得各类InSAR观测值的单位权中误差向量为

(13)

根据方差分量估计算法可知，当各类观测值单位权中误差近似相等时，即

(14)

(15)

当上述迭代过程结束时，可根据各个单位权中误差及其相应权阵计算各类观测值对应的方差矩阵为

(16)

同时，通过式(12)即可得到三维地表形变的最优估值。

在观测值方差已知的前提下，根据误差传播定律可以计算得到三维形变的方差。在使用误差传播定律时，应尽可能地保证观测值之间相互独立。因此，在计算三维地表形变的方差时，本文将观测方程式(11)简化为

(17)

即仅考虑了观测值与三维地表形变之间的几何关系。此时，三维地表形变的方差矩阵为

(18)

针对不同像元重复利用本章介绍的算法，直至所有像素点的三维形变解算完成。

1.3 顾及断裂线和自适应变化窗口的观测值选取策略

在利用SM-VCE方法进行三维地表形变求解时，其中的关键一步是选取一定大小窗口范围内的像素点建立观测方程。原始SM-VCE方法选取固定大小窗口(如15×15像素)内的所有点来建立观测方程[35]。这种选取方法在地表形变未发生跳变或者原始形变观测值相干性较高时可以得到较高精度的三维形变[36]。然而，地震导致的地表形变在震中区域往往会发生跳变，且精度较高相位观测值(如DInSAR获取的形变观测值)会发生失相干现象。鉴于此，本文提出顾及断裂线和自适应变化窗口的观测值选取策略，旨在获取可靠的近场三维地表形变结果。

根据光学影像或者POT获取的形变数据人工勾勒出地震发生后导致的主要断裂线，基于该断裂线即可排除窗口内与中心点不在断层同一侧的像素点。在设置初始窗口大小为S×S像素的情况下，当解算震中某像素点的三维形变时，S×S像素大小的窗口内的失相干像素将会被掩膜，以保证解算结果的精度和可靠性。但是，在这种情况下，精度较高的相位观测值的像素点个数可能少于某一阈值Kthr，进而无法得到可靠的三维形变结果。本文提出的窗口大小自适应变化的观测值选取策略将会在Kj

结合已有的地表应力应变模型在火山[35]、地震[31,40]方面的研究发现，参与SM建模的像素点个数大于200之后，即可得到较为稳健的未知参数平差值。基于此，本文设置窗口大小的初值S=15像素，像素点个数阈值Kthr=200。为了确定Sthr的取值大小，本文基于InSAR观测值进行变异函数拟合[40,42]

C(D)=C0·e-D/D0

(19)

式中，C(D)表示距离为D的像素之间的协方差；C0代表区域内的方差；D0即为相关距离。事实上，当两点间距离大于D0时，根据相关距离及像素空间分辨率即可确定不同InSAR观测值对应的阈值Sthr,j。最终的阈值Sthr为

Sthr=min([Sthr,1Sthr,2…Sthr,j…Sthr,J])

(20)

式中，min(■)代表取最小值。本文取一系列Sthr,j的最小值作为最终阈值Sthr的取值，从而保证每一类观测中距离为Sthr的两点对应的观测值都是相关的。在本文中，利用式(19)、式(20)对真实数据进行拟合，计算得到Sthr=63像素。

1.4 基于迭代加权最小二乘的窗口内同一类观测值自适应定权方法

Lk=[∂L/∂xe,∂L/∂xn,∂L/∂xu]·Δk+L0

(21)

为了简单起见，这里省略了下标j。顾及周围的K个点，则可得

L=Bsm,L·lL

(22)

式中

lL=[∂L/∂xe∂L/∂xn∂L/∂xuL0]T

(23)

(24)

(25)

式中，∘代表哈达玛运算；G为对角矩阵，其对角线元素可认为是降权因子。G的第k个对角线元素gk可基于以下权函数确定

(26)

(27)

(28)

2 模拟试验

本文利用2019年Ridgecrest地震的已有的断层几何和滑动分布数据[45]正演了三维地表形变场，整个形变场的大小为890×1120像素，像素分辨率约为100 m×100 m。根据真实数据的卫星成像几何参数(表1)计算得到InSAR视线向与方位向的形变信息。在此基础上，加入不同量级的高斯噪声(升轨DInSAR视线向、升轨POT方位向、升轨POT视线向、降轨DInSAR视线向、降轨POT方位向、降轨POT视线向的高斯噪声均方差分别为5、300、100、5、300和100 mm)，同时在不同的观测值上随机挑选了5%的像素点加入了粗差。需要说明的是，模拟试验中粗差的模拟方法为：将整个影像上随机选取一定比例p=5%的像素作为粗差点，粗差点上的数值为整个形变场中最大形变值的±m倍，其中m的取值与预设的M值有关，不同粗差点处m的取值不一定相同，但所有粗差点处m的取值服从[M-5,M+5]区间内的均匀分布，并且不同粗差点处m取值的正负也是随机的。在此次模拟试验中M=5。另外，为了更加符合真实数据中的InSAR噪声源，本文在升降轨InSAR视线向的模拟数据中加入了模拟的大气噪声，其中大气噪声是利用代尔夫特理工大学开发的Matlab工具包中的分形函数fracsurf进行模拟的，分形维度为2.2，大气噪声的范围为[-2.8 cm,2.8 cm]，在C波段哨兵数据的干涉图上约为2个条纹。通过以上步骤即可得到用于模拟实验的InSAR观测数据(图2)。

表1 真实试验中所用到的哨兵1号卫星SAR数据基本信息

注：红色加号的行列号为(454, 357)，是图4中15×15窗口内中心像素的位置。黑色折线为断裂线。图2 模拟试验中所用的InSAR观测值Fig.2 The InSAR measurements in the simulation

随后，分别利用原始SM-VCE方法和本文的SM-VCE方法计算了模拟数据的三维形变，如图3(d)—(f)和(j)—(l)所示。整体而言，两种方法得到的三维形变与模拟数据图3(a)—(c)较为一致。图3(g)—(i)和图3(m)—(o)分别给出了两种方法得到的三维形变残差图。可以看出，在断层破裂区域，原始SM-VCE方法仅能利用精度较低的POT方法进行三维形变解算，且未区分窗口内断层两侧的异质点，使得该区域的三维形变结果过于平滑、残差值较大。本文的SM-VCE方法则通过窗口大小自适应变化使得高精度的InSAR形变观测值也可以参与解算，同时通过顾及断裂线位置剔除了窗口内的异质点，最终得到的断层破裂区域的三维形变结果残差更小、精度更高。通过两种方法在断层破裂区域的对比，说明了本文提出的顾及断裂线和自适应变化窗口的观测值选取策略的优越性。

注：黑色折线为断裂线。图3 模拟试验中不同方法得到的三维形变场Fig.3 The 3D deformation obtained by different methods in the simulation

表2 模拟试验中不同方法得到的三维地表形变均方根误差对比

此外，鉴于IWLS主要是为了抑制粗差观测值对三维形变结果的影响，本文进一步利用模拟试验研究了不同粗差量级(即不同的M取值)和不同粗差个数(即不同的p取值)情况下，利用本文的IWLS方法获取三维形变结果的精度(图6)。为了避免大气噪声的影响，此处的模拟试验中未添加大气噪声。不难发现，随着观测值中粗差比例的增加，得到的三维形变精度在不断下降。当粗差比例小于等于40%时，不同粗差量级的模拟试验均可得到较为可靠的形变结果。而且，当粗差比例小于等于40%时，粗差量级越大得到的三维形变结果精度越高。这是因为粗差量级越大，IWLS方法识别出粗差的正确率就会越高。但是，当粗差观测值比例过大时(例如大于40%)，本文提及的IWLS方法也无法十分有效地抑制粗差观测值对三维形变结果的影响。此外，南北向和垂直向形变结果均方根误差呈现先减小后增大的趋势(拐点出现在粗差比例约为20%时)。这是因为南北向形变主要是由POT方位向观测值计算而来，并且在计算三维形变时，垂直向形变结果受南北向形变影响较大[46]，而POT方位向观测值本身的方差就比较大，当粗差比例较小时，IWLS方法对粗差的敏感度并不高，随着粗差比例的增大，使得观测值噪声的分布更加非高斯性，进而IWLS可以更为准确地识别出粗差、提高结果精度。而当粗差比例过大时，形变结果的精度会不断降低。为了验证上述观点，本文将模拟试验中DInSAR和POT观测值中的高斯噪声方差均设置为5 mm，得到的南北向和垂直向均方根误差变化曲线和东西向变化曲线较为一致，未发现局部极小值的现象。

需要说明的是，SM-VCE方法在获取三维同震地表形变时的优势主要有两点：①基于地表应力应变模型构建了邻近点三维形变之间的模型关系，有效抑制了空间高频噪声，使得引入力学模型的InSAR三维地表形变估计函数模型更加符合实际情况；②顾及地表应力应变模型建立的InSAR三维地表形变估计函数模型显著增加了多余观测个数，使得利用方差分量估计进行验后迭代定权成为可能，相比于传统的先验定权方法，方差分量估计算法可以显著提高多源观测数据的定权精度，进而可以获得更加精确、可靠的三维地表形变场。相对于传统仅利用观测值与三维形变之间几何关系进行求解的方法[17]而言，SM-VCE方法通过顾及形变在空间上的关联关系以及观测值之间的准确定权，可以得到更为精确的三维地表形变结果。

注：观测值的权重是无量纲。图4 行列号(454, 357)处15×15窗口内(a)—(f)6类观测值数据及(g)—(l)IWLS得到的每一类观测值内部的相对权重Fig.4 (a)—(f) The six kinds of measurements and (g)—(l) the corresponding weights from the IWLS method in the 15×15 window of the pixel (454, 357)

另外，基于地表应力应变模型构建形变在空间上的关联关系在一定程度上可以认为是一个滤波过程，进而会导致SM-VCE方法对一些空间上的局部高频形变信号不敏感。但是，地表应力应变模型作为一个力学模型，用于三维形变求解时，肯定比纯数学的滤波效果更好。此外，传统逐像素求解的WLS方法，由于缺少多余观测数据，难以利用VCE算法进行验后精确定权，对于SM-VCE方法而言，正是因为顾及地表应力应变模型建立的函数模型具有较多的多余观测数据，才能够利用方差分量估计算法进行观测值精确定权。

同时，SM-VCE方法中的窗口大小是一个必不可少的参数。理想情况下，窗口范围内的形变应满足地表应力应变模型，即形变梯度应为一个常数。因此，窗口的范围应远小于感兴趣的形变场尺度。对于高空间分辨率影像，或者大尺度形变场而言，本文中所述的15×15像素，甚至63×63像素均可以得到可靠的三维地表形变。对于本文提出的顾及断裂线和自适应变化窗口的观测值选取策略，当面对浅部发生的中小型地震时，由于断裂带小范围产生的复杂形变特征造成窗口内的观测值呈现异质性，此时应减小窗口阈值Sthr以保证近场三维地表形变结果的可靠性。

注：σj的单位为m。图5 行列号(454, 357)处方差分量估计迭代过程中的单位权中误差收敛Fig.5 Convergence of the variance component estimation at the pixel (454, 357)

注：均方根误差的单位是m。图6 不同粗差量级(即不同的M取值)和不同粗差比例(即不同的p取值)情况下，利用本文SM-VCE方法获取的模拟试验三维形变结果的精度对比Fig.6 The root mean squares errors of the three-dimensional deformations estimated by the enhanced SM-VCE method in this paper with the simulated experiments under different magnitude and different proportion of the gross error

在观测值噪声服从理想的高斯分布时，方差分量估计可以较为准确地得到各类观测值的验后权重。但是，真实的形变观测值往往会包含各种复杂的误差源信号(例如大气噪声、地形残差、轨道误差等)，在个别情况下，多源误差信号整体上表现为严重的非高斯性，进而使得方差分量估计难以收敛，甚至出现负方差，因此，需要针对个别情况进行人为干预的先验定权，或者引入更为稳健的算法来抑制相关误差源的影响。

3 2019年Ridgecrest地震的InSAR三维地表形变场

3.1 研究区域及所用数据

2019年7月6日(UTC),美国加州Ridgecrest地区(震中35.770°N，117.599°W)发生了M7.1级地震，震源深度约为8.0 km。Ridgecrest地震的发震断层位于东加利福尼亚剪切带上(east California shear zone，ECSZ)。该区域是圣安德烈亚斯断层(San Andreas fault，SAF)南段以东的一个地震活跃区域，与莫哈韦沙漠大致重合，具有多个与SAF平行的右旋走滑断层。因此，对加州地震展开地表形变监测可为研究该地震，ECSZ及SAF的构造运动提供基础资料。

针对2019年Ridgecrest地震，本次研究使用了升降轨哨兵1号卫星SAR数据(如表1和图7)。基于瑞士的GAMMA软件，分别利用DInSAR和POT技术进行了数据处理。其中，SRTM(shuttle radar topography mission) 30 m分辨率的DEM数据被用于去除地形相位及地形起伏引起的配准误差。为了抑制失相干噪声、提高计算效率，采用了32×8的多视因子(距离向×方位向)。在DInSAR处理过程中，利用基于最小二乘的滤波方法对多视后的干涉图进行滤波处理[47]，基于最小费用流方法对滤波后的干涉图进行解缠[48]。在POT的数据处理过程中，哨兵数据的匹配窗口为128×128像素，过采样因子为2。

图7 研究区域与数据覆盖情况Fig.7 Shaded relief map of the study area

图8为本文获取的6个InSAR形变观测数据。可以看出，Ridgecrest地震震中的形变可达米级。由于震中地表破坏较为严重，DInSAR获取的形变观测结果在断裂带附近失相干现象较为严重。相反，POT技术可以较好地抵制失相干现象，获得了较为完整的地表形变场。但是，由于POT技术是根据强度匹配获取地表形变，该技术往往比DInSAR技术获取的形变精度低。因此，尽管融合升降轨POT数据即可获得较为完整的三维形变，本文仍有必要通过窗口大小自适应变化的观测值选取策略来使得更高精度的DInSAR形变观测数据参与断裂带附近的三维形变解算。此外，图8(b)、(c)、(e)和(f)显示，基于POT获取的地表形变观测数据噪声较大。图8(b)、(c)和(e)中的形变场左侧还出现了明显的条带信号，通过对比光学影像发现，条带信号与光学影像中的公路较为相似，并且公路走向与卫星飞行方向越相近，POT方位向观测值中的条带信号越强，进而推测此条带信号是因为公路对电磁波反射信号过强，导致影像匹配时公路上的强反射信号沿公路走向较为相似，进而引起了配准错误。本文则基于IWLS方法，通过测量平差数据处理过程中的权重因子来抑制相关噪声(包括图8(b)、(c)和(e)中的条带信号)对SM-VCE方法解算三维形变的影响。

3.2 InSAR三维地表形变场

与模拟试验类似，本文分别利用原始SM-VCE方法和本文SM-VCE方法获取了该地震的三维地表形变场(图9)。可以看出，不同方法得到的三维地表形变具有高度一致性。其中，水平向和垂直向的最大形变分别约为2.1 m和0.38 m。然而，由于原始SM-VCE方法在断层附近无法兼顾高精度的相位观测数据，且未区分窗口内断层两侧的异质点，图9(d)—(f)中，由该方法得到的三维地表形变在断层附近的三维形变场过于平滑，严重偏离真实情况。从图9(g)—(i)，本文SM-VCE方法和原始SM-VCE方法得到的三维形变场差异图中可以看出，两种方法得到的三维形变在大部分区域是基本相同的，在断裂带附近差异最大，说明了本文SM-VCE方法获取可靠断裂带附近三维形变的优越性。同时，在远离断裂带区域(如图9(h)中主断裂带的西南侧)，本文SM-VCE方法可以较好地抑制一些局部误差信号对三维地表形变结果的影响，说明了IWLS在抑制粗差观测值方面的优势。

此外，为了更加突出窗口自适应选点方法的优越性，本文仅基于DInSAR获取的视线向和POT获取的方位向4个形变观测值，利用本文的SM-VCE方法获取了该地震的三维地表形变结果，如图10(a)—(c)所示。与图9(a)—(c)基于6个观测值的结果对比发现，二者的差异几乎可以忽略，如图10(d)—(f)，尤其在断层破裂带区域，在DInSAR完全失相干的情况下，本文的SM-VCE方法也可得到较为可靠的三维地表形变结果。然而，在DInSAR失相干区域，原始SM-VCE方法在15×15像素固定窗口内仅有两个观测值(即升降轨获取的POT方位向的形变观测值)，导致最终的三维地表形变场图10(g)—(i)在DInSAR失相干区域必定会缺失，并且在断裂带附近，由于15×15像素固定窗口内的像素点个数较少，三维形变结果会出现较多的跳变信号。尽管在窗口大小为49×49像素情况下，原始SM-VCE方法可以得到完整的三维形变场，如图10(j)—(l)，但是在断裂带附近，由于窗口内包含了较多的异质点，原始SM-VCE方法得到的三维形变结果会出现较多的错误形变信号，进一步说明了本文方法在获取高精度、完整三维同震形变场的优势。

注：正值代表地面朝向卫星运动，或沿卫星飞行方向运动,紫色曲线代表断层线。图8 不同数据不同方法获取的Ridgecrest地震InSAR形变观测值Fig.8 The InSAR measurements of the Ridgecrest earthquake

图9 在使用升轨DInSAR视线向、升轨POT方位向、升轨POT视线向、降轨DInSAR视线向、降轨POT方位向、降轨POT视线向等6个观测值情况下，不同方法得到的Ridgecrest地震三维地表形变Fig.9 The 3D deformation of the Ridgecrest earthquake obtained by different methods based on the six measurements in Fig.8

图10 在使用升轨DInSAR视线向、升轨POT方位向、降轨DInSAR视线向、降轨POT方位向等4个观测值情况下，不同方法得到的Ridgecrest地震三维地表形变Fig.10 The 3D deformation of the Ridgecrest earthquake obtained by different methods base on the four measurements in Fig.8(a), (b), (d), and (e)

需要说明的是，本文所提出的自适应选点方法主要是针对位于发震断层处的地震破裂带区域。一般情况下，结合历史资料、实地调研及大地测量等数据可以较为准确地获取发震断层断裂线数据，进而利用本文的自适应选点方法可以获取可靠的近场三维同震形变场。实际地震中，部分非发震断层在受到应力挤压、拉张等作用下，也可能表现出形变场不连续的现象。与发震断层处不同的是，InSAR技术在这些区域往往可以获取完整的地表形变场，进而基于断裂线的人工识别即可得到可靠的三维地表形变场。对于Ridgecrest地震而言，真实地表破裂较多[49]，在本文中为了简单起见，仅使用了两个主要的发震断层破裂线计算三维同震形变场来说明本文方法的优越性。

此外，本文利用Nevada Geodetic Laboratory(http:∥geodesy.unr.edu/)免费公开的5个近场GNSS站点数据(站点分布如图7所示)对两种方法得到的三维地表形变进行了精度评估。原始SM-VCE方法在东西向、南北向和垂直向形变精度分别为7.8、6.5和2.6 cm，而本文SM-VCE方法则分别为7.9、5.6和2.6 cm。可以发现，本文方法相比于原始方法仅在南北向分量有一定的改善效果。与模拟试验对比发现，真实观测数据中的粗差量级和个数相对较小，仅在POT方位向形变观测值中具有较为明显的粗差观测值信号，因此，在与GNSS对比中，本文SM-VCE仅在南北向分量改进较为明显。此外，InSAR数据和GNSS数据的时空分辨率均不相同，GNSS数据本身就包含误差，因此，利用GNSS进行InSAR精度评估并不完全可靠[16]。

4 总 结

融合多源InSAR形变观测数据是实现InSAR三维同震地表形变监测的基本途径，其中三维同震地表形变求解过程中建立的函数模型与随机模型对三维形变结果的精度影响至关重要。本文详细介绍了一种基于地表应力应变模型(SM)和方差分量估计(VCE)的InSAR三维地表形变监测方法(SM-VCE)。该方法通过顾及一定窗口范围内不同像素点三维形变之间的力学关系(即SM)建立函数模型，进而利用VCE算法实现多源InSAR形变观测数据随机模型的准确确定。相比于传统的单点求解方法，SM-VCE方法可得到更为准确的三维地表形变结果。

为了使得该方法可以得到更为精确的、完整的三维同震地表形变场，本文进一步提出了顾及断裂线位置和窗口大小自适应变化的观测值选取策略和一种基于迭代加权最小二乘的窗口内同一类观测值相对定权方法。其中，观测值选取策略是为了解决断裂带附近InSAR形变观测数据易缺失且不连续导致的SM-VCE方法精度较低和形变结果缺失问题，而同一类观测值内部相对定权方法是通过迭代加权最小二乘方法抑制窗口内的异常形变观测值，进而提高三维形变的解算精度。

模拟试验和基于升降轨哨兵数据2019年Ridgecrest地震三维形变的研究表明，本文提及的方法可得到更为精确的三维同震形变结果，尤其在断层破裂带附近，窗口优化的SM-VCE方法可在数据缺失的情况下，利用邻近的有效InSAR形变观测值恢复得到可靠的三维地表形变场。同时，本文通过模拟试验验证了观测值中不同的粗差比例对三维形变结果的影响。结果表明，当粗差比例小于等于40%时，窗口优化的SM-VCE方法通过迭代加权最小二乘可显著抑制粗差观测值对结果的影响。

本文介绍的SM-VCE方法也可联合星载、地基InSAR面观测数据和GNSS、水准等传统大地测量点观测数据实现多种地质灾害(例如火山喷发、滑坡等)的三维地表形变准确测量。同时，在独立观测数据较少的情况下，例如只有升降轨的DInSAR观测数据，也可利用SM-VCE方法实现二维地表形变测量，或者，通过引入可靠的先验信息约束(例如坡度约束、地下开采沉陷模型约束)实现三维地表形变测量。此外，本文介绍的窗口大小自适应变化的准则是窗口内的有效观测值个数，只是从数值解算的稳定性方面进行考虑，未来可发展基于数据本身或者局部实际地表形变梯度的窗口大小自适应确定方法，以提高三维地表形变的精度及可靠性。