大地测量反演是大地测量研究的主要方式之一，即利用大地测量技术(如GPS观测、水准测量、重力测量等)测得的地表形变来解释地壳深部断层滑动的动力过程。利用多种数据联合反演能更客观、准确地反映地壳深部断层运动特征，进而对其真实运动规律进行研判和解释，国内外学者在此方面已有一定的研究成果。Lisowski等^[1]利用多种大地测量数据对1989年加州洛马普列塔地震同震滑动进行联合反演; 李爽等^[2]研究多种数据联合反演模式及算法; 刘琦等^[3]利用GPS、水准和强震动观测资料联合反演2013年芦山7.0级地震同震滑动分布; Zhao^[4]利用重力和GPS数据对红河断裂带活动进行反演计算。本文将一种新的群智能算法——微粒群优化算法(particle swarm optimization，PSO)引入到大地测量联合反演非线性算法的研究中^[5]，并基于GPS和水准数据对祁连山北缘断裂三维滑动速率进行反演。

1 带惯性权重微粒群算法的原理及方法

微粒群优化算法基本原理来源于人工生命，假设有n+1个微粒，微粒群以特定速率在D维空间中运动，运动速率根据个体微粒运动速率与微粒群运动速率经验进行实时动态调整。

在微粒群进化方程中引入惯性权重因子w以优化微粒群算法，此时微粒i在D维中的位置x与速度v的更新公式为^[6]:

$ \begin{array}{l} {v_{ij}}(k + 1) = w{v_{ij}}(k) + {c_1}{r_{1j}}\left[ {{p_{ij}}(k) - {x_{ij}}(k)} \right] + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{c_2}{r_{2j}}(k)\left[ {{g_j}(k) - {x_{ij}}(k)} \right] \end{array} $

(1)

$ {\mathit{x}_{\mathit{ij}}}{\rm{(}}\mathit{k}{\rm{ + 1) = }}{\mathit{x}_{\mathit{ij}}}{\rm{(}}\mathit{k}{\rm{) + }}{\mathit{v}_{\mathit{ij}}}{\rm{(}}\mathit{k}{\rm{ + 1)}} $

(2)

式中，w为惯性权重，使微粒保持运动惯性，在算法开始时，赋予w一较大正值，用以保证微粒在算法初始阶段能以较大的速度步长在整个空间内搜索优秀种子，w在搜索过程中逐渐线性减小，保证微粒在后期搜索中更容易得到全局最优解^[7]; c₁、c₂为微粒运动的加速度，用以调整全局最优微粒和个体最优微粒的最大运动步长^[8]; r_1j、r_2j的区间范围为[0, 1];p_ij=(p_i1, p_i2, …p_iD)为个体微粒i的最优位置; g_j=(g₁, g₂, …g_D)为微粒群的全局最优位置。V_i=(v_i1, v_i2, …v_iD)为微粒i当前飞行的速度; X_i=(x_i1, x_i2, …x_iD)为微粒i当前的位置。

2 矩形位错理论模型

断层参数(包括断层倾角φ、断层长度L、断层宽度W、断层深度d，矩形断层上盘相对于下盘分别在走向、倾向和法向的滑动量U_i(i=1, 2, 3)^[9])与地面上GPS观测的三维位移速率场之间的关系可通过如下矩形位错理论模型来表示^[10-11]。图 1中，X轴指向断层走向，Z轴指向断层的垂线方向，Y轴构成右手坐标系。

GPS观测得到的地面三维位移场速率u_i(i= x, y, z)可以写为:

$ {u_i} = f\left( {\xi , \eta , \zeta , {U_j}, L, W, d, \varphi } \right) $

(3)

3 算例及结果对比分析 3.1 GPS、水准数据联合反演

为了验证微粒群算法在联合反演中的有效性和稳定性，深入理解实际三维形变场与断层深部活动的关系，利用位错模型引起的地表水平位移和高程变化的解析式，结合祁连山北缘断裂的地质背景及青藏高原东北缘2001~2004年GPS数据和2000~2006年水准数据，对断裂的三维滑动速率进行联合反演。建立2种数据联合反演模型的关键步骤之一就是确定2种数据在反演中所占的权重^[12]。本文采用2种定权方案:方案1为常规定权法，即根据粗略的先验信息，较主观地取λ∈(0, 1)内的某一数值; 方案2为把λ作为未知数，与断层参数一同反演，确定权重。

联合反演模型(目标函数)为:

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;{\lambda _{{\rm{GPS}}}}\sum\limits_{{\rm{GPS}}} {{{\left( {\frac{{{u^{{\rm{obs}}}} - u(m)}}{{{u^{{\rm{obs}}}}}}} \right)}^2}} + \\ \left( {1 - {\lambda _{{\rm{GPS}}}}} \right)\sum\limits_{{\rm{leveling }}} {{{\left( {\frac{{{h^{{\rm{obs}}}} - h(m)}}{{{h^{{\rm{obs}}}}}}} \right)}^2}} = \min \end{array} $

(4)

式中，λ_GPS为GPS观测得到的水平位移在目标函数中所占权重; u^obs、h^obs分别为GPS和水准测量的实测值; u(m)、h(m)分别为利用位错理论模型计算出的地表水平位移和高程变化值。

为使模型更接近真实断裂带，采用微分思想，将复杂断裂看成是多个小的矩形断层组合^[12]，根据断裂带的分布排列等几何特征进行断裂形态的几何分段。

方案1中，给定λ_GPS的值为0.5，利用微粒群算法进行联合反演，反演结果见表 1。

方案2中，将权λ_GPS限定在[0, 1]区间内作为未知参数与断层三维滑动速率一同反演，即要求微粒群在四维空间(U₁, U₂, U₃, λ_GPS)中搜索，最终输出反演结果见表 2。

3.2 其他算法反演

表 3为蚁群算法联合反演断层滑动速率的反演结果。从表 1~3不难发现，祁连山北缘断裂主要滑动特征为右旋走滑和逆冲，且空间分布不均匀。断裂右旋走滑和逆冲2个分量明显大于挤压分量，与前人结论^[13-16]基本一致。比较以上反演结果, 并列于表 4。

由表 4看出，无论是数值还是数据的稳定性和收敛程度，微粒群算法反演结果都优于蚁群算法。而方案2通过灵活调整观测值所占权重，更合理地反映了不同观测值在不同区域的响应大小，其结果优于方案1。

4 结语

1) 微粒群算法收敛速度快、稳定性高，结合经典的矩形位错理论模型，是一种可以有效求解断层三维滑动速率反演问题的优化算法。

2) 祁连山北缘断裂滑动空间分布不均匀，其主要运动特征是以右旋走滑和逆冲为主，同时兼有相互挤压。

3) 在多种数据联合反演计算中，不同类别的测量数据所占权重不同。以水准数据和GPS数据为例，若二者分布大致均衡，在联合反演中，GPS数据所占的权重远大于水准数据，也就意味着以走滑和倾滑为主的断层深部运动带来的地面响应以水平运动为主。同时，顾及权比的方法比常规定权的方法得到的结果与地质资料更吻合。