在同震滑动分布反演中，地表同震形变量与断层滑动量之间为线性关系^[17-18]，可表示为：

式中，$\mathit{d}$表示观测向量（如InSAR视线向位移、GPS位移场等）；$\mathit{G}$表示格林函数矩阵；$\mathit{m}$为滑动参数；$\mathit{\epsilon }$为InSAR数据或GPS数据等观测值的随机观测误差。

对于多类相互独立的地表形变观测数据，则式（1）可以表示为：

在同震滑动分布反演中，为了得到更精细的滑动，通常将断层面划分为n个矩形单元，基于矩形位错模型计算矩形单元的单位走滑量和单位倾滑量引起地表k个观测点的形变值，构成k×2n的格林函数矩阵$\mathit{G}$。断层面离散化处理后，会导致$\mathit{G}$严重秩亏，此时式（1）的法矩阵病态，从而使滑动分布解不稳定。为了避免这种情况，可对断层的位错施加一定的约束，一方面增加虚拟观测数据避免系数矩阵的秩亏性，另一方面抑制相邻断层单元出现大的梯度变化^[17]。

本文采用拉普拉斯约束对断层单元进行平滑约束，使：

式中，$\mathit{H}$为拉普拉斯二阶平滑矩阵。联合式（2）和式（3），得到滑动分布联合反演公式为：

式中，$\mathit{P}$为联合反演的权阵；${\mathit{D}}_i$为各类数据的方差阵；$\mathrm{1}/{\sigma }_{\mathrm{1}}^{\mathrm{2}}$、$\mathrm{1}/{\sigma }_{\mathrm{2}}^{\mathrm{2}}$$\cdots$$\mathrm{1}/{\sigma }_n^{\mathrm{2}}$表示各类观测数据的相对权重；$\mathrm{1}/{\alpha }^{\mathrm{2}}$表示模型平滑度与数据拟合度之间的相对权比。式（5）的最小二乘解为：

采用方差分量估计法对式（5）进行迭代计算，可以求解各类观测值的方差分量，从而确定联合反演的相对权比与正则化参数。然而在滑动分布反演中，利用方差分量估计方法确定相对权比时需要初始值进行迭代，而未加入平滑约束的滑动分布解是不可靠的，故将有偏的滑动分布解作为初始值进行迭代运算。而此时需要考虑其偏差的影响，且已有文献并没考虑滑动分布解是有偏的，因此本文利用偏差改正的方差分量估计进行求解。这里仅列出GPS数据与InSAR数据联合反演的公式，误差方程可表示为：

所求方差为：

其中，

根据方差分量估计公式，可得：

式中，$E$（）表示期望运算；$\mathrm{t}\mathrm{r}$表示矩阵的求迹运算。

同理可得到$E({\mathit{v}}_{\mathrm{2}}^{\mathrm{T}}{\mathit{P}}_{\mathrm{2}}{\mathit{v}}_{\mathrm{2}})$与$E({\mathit{v}}_{\mathrm{3}}^{\mathrm{T}}{\mathit{P}}_{\mathrm{3}}{\mathit{v}}_{\mathrm{3}})$的表达式，将数学期望的符号去掉，则可以求出方差分量的估值，其矩阵形式为：

式中，

其中，$A=n_{\mathrm{1}}+\mathrm{t}\mathrm{r}({\mathit{N}}^{-\mathrm{1}}{\mathit{N}}_{\mathrm{1}}{\mathit{N}}^{-\mathrm{1}}{\mathit{N}}_{\mathrm{1}})-\mathrm{2}\mathrm{t}\mathrm{r}({\mathit{N}}^{-\mathrm{1}}{\mathit{N}}_{\mathrm{1}})$；$B=n_{\mathrm{2}}+\mathrm{t}\mathrm{r}({\mathit{N}}^{-\mathrm{1}}{\mathit{N}}_{\mathrm{2}}{\mathit{N}}^{-\mathrm{1}}{\mathit{N}}_{\mathrm{2}})-\mathrm{2}\mathrm{t}\mathrm{r}({\mathit{N}}^{-\mathrm{1}}{\mathit{N}}_{\mathrm{2}})$；C=$n_{\mathrm{3}}+\mathrm{t}\mathrm{r}({\mathit{N}}^{-\mathrm{1}}{\mathit{N}}_{\mathrm{3}}{\mathit{N}}^{-\mathrm{1}}{\mathit{N}}_{\mathrm{3}})-\mathrm{2}\mathrm{t}\mathrm{r}({\mathit{N}}^{-\mathrm{1}}{\mathit{N}}_{\mathrm{3}})$。

滑动分布解的偏差为：

式中，$\widehat{\mathit{m}}$为无偏估计值，由于其无法获取，故采用正则化解替代，从而得到滑动分布解的偏差计算如下：

由于每次迭代计算时，所代入的滑动分布解都是有偏的，因此需要进行偏差改正，改正后的滑动分布解公式为：

将偏差改正的滑动分布解作为初值进行迭代，每次迭代后都进行偏差改正，从而得到偏差改正的方差分量估计。

文献[11]提出对参数的残差进行偏差改正，其中参数表达式为：

得到参数残差的偏差为：

式中，$\lambda$表示正则化参数；$\mathit{R}={\mathit{H}}^{\mathrm{T}}\mathit{H}$为正则化矩阵。

根据方差分量估计公式，可得：

同理可得$E({\mathit{v}}_{\mathrm{2}}^{\mathrm{T}}{\mathit{P}}_{\mathrm{2}}{\mathit{v}}_{\mathrm{2}})$表达式。将式（16）数学期望的符号去掉，则可以求出方差分量的估值，得到其矩阵形式为：

式中，

其中，$D=n_{\mathrm{1}}-\mathrm{2}\mathrm{t}\mathrm{r}({\mathit{N}}_{\lambda }^{-\mathrm{1}}{\mathit{N}}_{\mathrm{1}})+\mathrm{t}\mathrm{r}({\mathit{N}}_{\lambda }^{-\mathrm{1}}{\mathit{N}}_{\mathrm{1}}{)}^{\mathrm{2}}$；E=$n_{\mathrm{2}}-\mathrm{2}\mathrm{t}\mathrm{r}({\mathit{N}}_{\lambda }^{-\mathrm{1}}{\mathit{N}}_{\mathrm{2}})+\mathrm{t}\mathrm{r}({\mathit{N}}_{\lambda }^{-\mathrm{1}}{\mathit{N}}_{\mathrm{2}}{)}^{\mathrm{2}}$。

考虑到正则化参数会放大方差分量，采用MSE法求解正则化参数，并通过迭代求解得到最后的正则化参数和方差分量，迭代步骤如下：（1）采用MSE法确定滑动分布反演的正则化参数；（2）利用获得的正则化参数求得滑动分布解的有偏估计量；（3）计算出残差的偏差，对残差进行偏差改正，利用偏差改正后方差分量估计求得各类数据的方差分量；（4）根据确定的方差分量得到具有不同权比的滑动分布解，以此作为初值重复步骤（1）~（3），直到均方误差收敛，迭代停止。

为了更好地比较不同方差分量估计的适用性和合理性，采用3种方案进行模拟实验，如表 1所示。

为了验证基于偏差改正方差的分量估计方法的合理性与优势，本文进行了模拟实验，模拟一个双主滑动区的断层进行滑动分布反演。设置断层几何中心位置为X=0 km，Y=0 km，长度为30 km，宽度为30 km，走向角为50°，倾角为45°，滑动角为50°，断层顶深为0 km。首先模拟了断层的真实滑动分布情况，然后利用均匀弹性半空间地表位移与矩形位错公式^[19]计算了每个断层单元的滑动对地表417个站台的地表位移之和，并给每个台站的位移施加观测误差N（0，3² mm²），最终得到模拟形变数据。采用相同的方法获得592个InSAR形变位移，考虑到InSAR数据质量较差，对InSAR形变位移施加误差N（0，1² cm²），得到InSAR模拟形变数据。GPS与InSAR形变数据如图 1所示。在滑动分布反演中，将断层离散成20×20个1.5 km×1.5 km的断层单元，并利用偏差改正方差分量估计进行滑动分布反演。

图  1  模拟地表观测点

Figure 1.  Simulated Surface Observation Points

在模拟实验中，分别采用表 1中的3种方案确定多类数据的相对权比，方案1、2、3确定的GPS、InSAR、正则化参数的相对权比分别为1∶0.210 8∶0.001 5、1∶0.060 9∶0.001、1∶0.209 3∶0.000 9。图2（a）、2（b）、2（c）分别为采用方案1、2、3反演得到的滑动分布结果；图2（d）、2（e）、2（f）分别为方案1、2、3反演得到的滑动分布结果与真实滑移的差值分布。

图  2  模拟实验的滑移分布结果及其残差

Figure 2.  Slip Distribution Results and Residuals in the Simulation Experiment

由图 2可以看出，3种方案都能很好地反演出两个主滑动区，浅部滑动区的滑动能够较好地反演出来，而深部滑动区滑动量存在明显的低估，这与数据难以分辨出地下深层的滑动有关。对于两个主滑动区之间的区域，滑动量均被高估，方案3较方案2更接近模拟滑移，方案1结果最差；对于深部滑动区，方案1难以较好地反演出实际滑移，而方案2较方案1有所改进，方案3相比方案1、2，反演结果更接近于模拟滑动情况。

表 2统计了模拟实验中3种方案确定的相对权比以及反演得到的各参数结果和模拟值。从表 2可以看出，方案1反演的最大滑动量最小，方案2反演的最大滑动量要略高于方案1，方案3结果略高于方案2，与滑动分布图结果一致。对于最大滑动角，3个方案反演结果差别不大，而方案2结果更接近实际滑动角。对比3个方案反演的地震矩，方案1至方案3依次增大，这与最大滑动量情况一致。均方根（root mean square，RMS）误差作为大地测量反演中的一种衡量反演结果指标，3种方案反演结果均为6.3 mm，从RMS角度分析，3个方案的反演精度相当，但方案2和方案3考虑了反演中偏差的影响，理论更为严密。

2016-10-26意大利发生Mw 6.1地震，根据美国地震调查局（United States Geological Survey，USGS）网站发布的震源机制，其震中位置为（13.067°E，42.956°N），震源深度为10 km。127个Visso地震GPS台站同震形变数据来自网站http://ring.gm.ingv.it/，814个InSAR升轨位移数据来自文献[20]，利用获取的GPS和InSAR数据进行联合反演。首先采用多粒子群算法 ^[21]进行非线性反演，得到意大利地震断层的走向角为170.3°，倾角为40°，以及断层中心位置等参数，并根据断层参数来构建断层。将断层沿走向/倾向把断层长度、宽度分别拓展至27 km、15 km，并将断层面延伸至地表。将断层面均匀剖分成1 km × 1 km大小的矩形单元，共得到405个矩形单元。分别采用表 1的3种方案确定GPS和InSAR的相对权比以及正则化参数。

在Visso地震反演实验中，方案1、2、3得到的相对权比分别为1∶2.073 4∶0.002 6、1∶0.000 83∶0.001、1∶2.888 1∶0.001。图 3为利用3种方案反演Visso地震滑动分布结果，其结果与近年来部分学者反演Visso地震滑动分布各参数结果见表 3。其中，文献[22]联合了Sentinel-1和ALOS-2差分干涉数据，及震区GPS三维形变数据，反演得到滑动分布的最大滑动量为0.900 0 m；文献[20]利用Sentinel-1多景影像数据和127个GPS台站形变数据进行联合反演，得到最大滑动量为0.840 0 m。

图  3  3种方案反演的Visso地震滑动分布结果

Figure 3.  Results of the Slip Distribution of the Visso Earthquake for Three Schemes

由图 3可知，采用本文3种方案反演Visso地震时，其反演结果均显示Visso地震同震滑动分布区域主要分布在地下0~6 km，其最大滑动区域为4 km，与文献[20]研究结果较一致，方案1~3反演的最大滑动量逐渐增加，但相差不大，且滑动已延伸至地表，这与现场勘查的地表破裂情况吻合。由表 3可以看出，方案1、2、3确定的最大滑动量分别为0.605 6 m、0.623 9 m和0.692 6 m，均低于文献[20，22]的结果。这可能是由于两者数据源差造成的，也表明了多类数据能够更好地约束断层，提高反演精度。结合文献[20，22]结果发现，方案3较方案1、2能够更好地反演出最大滑动量。方案1、2、3反演的平均滑动角分别为-82.481 5°、-68.742 8°和-85.693 2°，均在各学者的反演结果范围，且方案2反演结果与文献[20]较为接近。方案1、2、3反演矩震级分别为Mw 6.216 2、Mw 6.189 9和Mw 6.242 9，其中方案2结果与文献[20]较为一致，而方案1和3均高于其他文献结果。方案1、2、3反演的RMS分别为7.2 mm、16.2 mm、6.8 mm，方案2结果较差，方案3要略优于方案1。综合以上分析，方案3在方案1的基础上进行偏差改正，两者反演结果较为一致，且方案3反演精度更高。而方案2采用了MSE法求解正则化参数，仅利用方差分量估计求解相对权比，反演结果与方案1和3存在差异，反演参数与文献[20]较为一致。

2016-10-30意大利继Visso地震之后又发生了Mw 6.6地震，根据USGS网站发布的震源机制，其震中位置为（13.096°E，42.862°N），震源深度为8 km。本文采用的Norcia地震114个GPS台站同震形变数据来自网站http://ring.gm.ingv.it/，采用的1 820个InSAR升轨数据来自文献[20]，并利用获得的GPS和InSAR进行同震滑动联合反演。同样采用§3.1中的非线性反演算法获得Norcia地震断层的走向角、倾角及断层中心位置参数。构建好断层之后，将断层沿走向/倾向把断层长度、宽度分别拓展至27 km、15 km，并将破裂面延伸至地表；将断层面均匀剖分成1 km×1 km大小的矩形单元，共得到405个矩形单元。分别采用表 3的方案确定GPS和InSAR的相对权比以及正则化参数。

在Norcia地震反演实验中，方案1、2、3确定的相对权比分别为1∶5.268 4∶0.001 1、1∶0.005 8∶0.01和1∶4.196 4∶0.000 39。图 4为利用3种方案反演Norcia地震滑动分布结果。

图  4  3种方案反演的Norcia地震的滑动分布结果

Figure 4.  Results of the Slip Distribution of the Norcia Earthquake for Three Schemes

由图 4可知，利用3种方案反演Norcia地震时，其反演结果均显示Norcia地震同震滑动分布区域主要分布在地下0~6 km，且地震破裂到地表，与文献[20]结果一致；3种方案反演结果较为一致，方案2反演的最大滑动量最大，其次为方案3和方案1。

本文统计了3种方案反演Norcia地震滑动分布各参数以及近年来部分学者反演Norcia地震滑动分布各参数，结果见表 4。由表 4可以看出，方案1、2、3反演的最大滑动量分别为2.013 6 m、2.484 4 m和2.267 3 m，均小于文献[20]结果，其中方案2与文献[22]结果接近，可能是因为融合了多类数据，能够更好地约束断层。方案1、2、3反演的平均滑动角分别为-87.342 6°、-84.770 7°和-84.894 7°，其中方案2、3反演的结果与USGS较为接近。方案1、2、3反演的矩震级分别为Mw 6.584 4、Mw 6.599 6、Mw 6.589 3，反演的地震矩分别为8.446 2×10¹⁸ N·m、8.900 8×10¹⁸ N·m、8.590 2×10¹⁸ N·m，略低于USGS公布的结果，方案1、3与文献[24]反演结果较为一致，方案2与文献[20，22]较为接近。方案1、2、3反演的RMS分别为13.4 mm、22.3 mm、12.7 mm，从RMS角度分析，方案2结果较差，方案3略优于方案1。

本文采用了3种方案对模拟实验、Visso地震和Norcia地震进行滑动分布反演，并对3种方案进行对比分析，综合实验结果发现，本文3种方案均能很好地反演出合理的相对权比与正则化参数。但无论是对参数进行偏差改正还是对残差进行偏差改正，都能够有效地提高最大滑动量，因为两者均考虑到反演中偏差的影响，但两种方案针对不同的实际情况，改进的程度有所不同。方案1、3反演的滑动分布较为一致，而方案2反演结果在主滑动区较为一致，在周围区域存在差别，这可能与方案2采用MSE法求解正则化参数相关。对于3种方案反演的矩震级，在模拟实验中，3种方案都能较准确地反演出实际的矩震级，而在本文所采用的两个实际地震中，方案2的反演结果更为准确，Visso地震反演实验中，方案1、3的结果要略高于其他文献结果。对于RMS值，方案3的结果明显优于其他两种方案，在两个实际地震中，方案2的RMS值明显高于方案1、3，由此也可以说明方案1、3反演的滑动分布更为合理。

本文分别采用了方差分量估计、对残差进行偏差改正的方差分量估计和对参数进行偏差改正方差分量估计法3种方法确定联合反演的相对权比和正则化参数，验证了3种方法的合理性与可行性，且进行偏差改正后的结果优于未进行偏差改正的方差分量估计法。针对不同偏差改正的方式在不同地震中的应用，其反演结果存在差异。针对残差进行偏差改正的初衷在于参数的偏差可能大于参数本身，但通过本文的实验表明，对参数进行偏差改正时，其反演结果改正较为显著，这也表明了偏差改正方差分量估计法在滑动分布反演中具有一定优势。