人类对地球内部物理性质和矿产资源分布的了解,大多来自地表地质和地球物理、地球化学资料的反演和解释(王家映,2000)。在科学研究和工程技术中,在由观测数据反推模型或体系时,都无法回避反演问题。绝大多数地球物理反演问题都属于非线性问题。一部分非线性反演问题可以通过线性化方法迭代求解,如Newton法、拟Newton法及其变形等; 另一类方法则是直接搜索的方法,这类方法又分为非启发式和启发式方法。非启发式的反演方法,如蒙特卡罗法(Monte Carlo method),计算量很大,效率较低。近年来,人们不断发展启发式的优化方法,并用于地球物理反演问题中,如遗传算法(万永革等,1997)、模拟退火算法(曹小林等,2000; 于鹏等,2007)、量子退火算法(魏超等,2006)、差分演化算法(李志伟等,2006)等。

粒子群算法(Particle Swarm Optimization)(Kennedy等,1995; 谢晓锋等,2003)是近年来发展起来的一种新的、有效的现代启发式算法,具有计算简便、收敛速度快、收敛准确等特点,在最优化问题中有一定的应用,但尚未见到用于地球物理反演问题的文献报道。本文中笔者首先介绍了粒子群算法的基本原理和步骤,然后建立利用余震数据拟合主震断层面参数的范数准则下的数学模型,最后根据2003年7月21日云南大姚6.2级地震及其余震的双差定位法定位结果(华卫等,2006),使用粒子群算法求解,通过研究震源点在断层面附近的分布规律,说明模型和算法的有效性。

粒子群优化算法是一种基于群智能的演化计算技术,它是Kennedy等(1995)受人工生命研究结果的启发,于1995年首先提出来的。该算法与遗传算法类似,它也是一种基于群体的优化工具,但与遗传算法相比,该算法具有收敛速度快,容易实现,而且又具有深刻智能背景的优点。

粒子群算法是根据对环境的适应度将群体中的个体移动到好的区域; 将每个个体看作D 维搜索空间中的一个没有体积的粒子(点),所有的粒子都有一个由被优化的函数决定的适应值,主要特点为:① 每一粒子都被赋予了初始随机速度并在解空间中流动; ② 个体具有记忆功能; ③ 个体的进化主要依靠它本身的“飞行经验”以及同伴的“飞行经验”进行动态调整,最终飞向最优位置,即通过迭代找到最优解(Parsopolos等,2002)。

设在D维空间中有m个粒子,第i个粒子表示为x_i=(x_i1,x_i2,…,x_iD)∈R^D,它所经历的最好位置(对应最优的适应值)为p_i=(p_i1,p_i2,…,p_iD),在群体中所有粒子所经过的最好位置的索引号用符号g表示,即p_g,粒子的速度v_i=(v_i1,v_i2,…,v_iD)。对于每一代粒子,其第d(1≤d≤D)维第i个粒子的速度和位置根据如下方程进行变化(Shi等,1998):

v_id=ω·v_id+c₁·rand₁()·(p_id-x_id)

+c₂·rand₂()·(p_gd-x_id); (1)

x_id=x_id+v_id.(2)

其中,rand₁()、rand₂()是[0,1]内的均匀分布的相互独立随机数; c₁和c₂为取值于[0,2]内的加速常数(Suganthn,1999),通常解释为“自身学习率”和“社会学习率”; ω为惯性权重因子,当权值较大时,粒子对未探测空间搜索能力强,局部细调能力弱,反之,权值较小时,算法的局部搜索能力增强,而搜索新空间的能力减弱,利用这样的特性,在计算的初期选取较大的权值以加大搜索区域,而在计算的后期,为了更好地进行局部搜索,选取较小的权值。Shi(2000)提出了一个模糊系统,Chatterjee(2006)提出一个规则来动态调节ω,但方法都过于复杂。一般认为,ω可以在[w_min,w_max]内线性递减取值,即

ω=((Iter_max-Iter)/(Iter_max))(w_max-w_min)+w_min.(3)

其中Iter_max表示最大迭代次数,由所研究的问题确定,Iter表示当前迭代次数。

此外,为了防止搜索发散(爆炸),当前粒子的加速导致它在某一维的速度超过了最大速率,则粒子速度v_id会被一个最大速率v_max所控制,即

v_id={v_max, v_id>v_max;

-v_max, v_id<-v_max.

粒子群算法的基本步骤如下:

(a)初始化设置:种群规模为N,随机初始化它们的位置x_i和速度v_i,最大迭代次数为Iter_max;

(b)设置每个粒子目前最优位置为p_i=x_i,并记p_g=max_ip_i为全局最优位置;

(c)计算每个粒子的适应度值,即对每个粒子计算其目标函数值;

(d)对每个粒子,将其适应度值与其历史最优位置进行比较,更新p_i;

(e)将每个粒子的当前最优位置p_i与群体历史最优位置p_g进行比较,更新p_g,并更新p_g的索引号;

(f)根据粒子群速度和位置更新方程(1)和(2)来调整微粒的速度和位置;

(g)检查终止条件通常为达到最大迭代次数或者足够好的适应值或者最优解停滞不再变化。若上述条件得以满足,则终止迭代,否则返回第2步。

地震震源机制解是描述地震的重要参数,在地震学研究中占有相当重要的位置,该参数又几乎是确定地壳深部应力场的唯一途径。目前地震的震源机制解往往通过P波初动(Reasenberg等,1985)、体波波形(姚振兴等,1985)、面波波形(马淑田等,1999)和大地测量数据(陈运泰等,1979)来测定。

大地震发生后,大量余震在发震断层面及其附近发生。因此余震震源位置的空间分布可以较为精确地勾画出发震断层面的形状和位置。假定地震发震断层可以用一个平面来模拟,且大多数余震发生在这个断层面的附近,则可以通过余震震源位置参数来求解发震断层的走向、倾角及位置。早在1992年,王鸣等(1992)就采用了这种原则,可是他们采用的Gauss-Newton法对初始解的选取比较敏感。万永革等(2008)利用唐山地震序列余震资料双差定位法精确定位的结果,采用模拟退火算法对唐山地震序列破裂面的几何参数进行求解,虽然他们的算法不依赖于初始解,但仍然采用l₂-范数得到评价解的总残差。由于l₂-范数准则在可用数据较少且分布较为散乱时,受野值点的影响较大,本文讨论分别采用l₁-范数和l₂-范数来评价解的总残差的结果。

本文给出根据余震位置确定主震断层面走向、倾角的粒子群算法。如果在地震之后能够快速确定余震的震源位置,就可以给出地震断层面的一种约束。

下面给出根据余震震源位置参数来求解发震断层参数的非线性模型。

设在地理坐标系中,(x_i,y_i,z_i)表示第i个余震的震源位置,地震断层面的走向为φ,倾角为δ,到坐标原点距离为ρ,则断层面法向量为(sinφsinδ,-cosφsinδ,cosδ)(万永革等,2000)。于是断层面在地理坐标系中的方程(王福昌等,2007)为

xsinφsinδ+y(-cosφ)sinδ+zcosδ-ρ=0.(4)

所以震源点到平面的垂直残差为

d_i=x_isinφsinδ+y_i(-cosφ)sinδ+z_icosδ-ρ.(5)

假定共有n个可以利用的余震数据样本,这可以建立以所有余震到断层面垂直距离残差向量d=(d₁,d₂,…,d_n)∈Rⁿ的p范数为目标函数的优化问题,即

E(ρ,φ,δ)=‖d‖_p=[∑ⁿ_i=1(d_i)^p]¹/p,p>0.(6)

我们的目标是找到一组符合地震学解释要求的值(ρ^*,φ^*,δ^*),使得E(ρ,φ,δ)在l_p-范数准则下在此处取最小值。在式(6)中,p=2时,该式即为王鸣等(1992)在研究中使用的最小二乘准则模型; p=1时,可以削弱较大残差的影响(Prugger等,1988; Guitton等,2003),从而在可用数据较少且数据分布较为散乱时,给出一个比最小二乘准则更稳健的参数估计值。

为了验证算法和准则的有效性,我们利用2003年7月21日云南大姚6.2级地震余震双差定

位数据来确定主震断层面参数。使用双差定位数据,在式(6)中,求出一组值(ρ^*,φ^*,δ^*),使得E(ρ,φ,δ)在l_p-范数准则下在此处取到最小值,再根据地震学解释,确定断层面的走向φ和倾角δ,最后确定断层的区域。

由于粒子群算法不要求计算目标函数的导数,因此无论p=1或p=2时,粒子群算法都同样有效。在粒子群算法中,种群大小设为50个,粒子的速率上界不超过20,最大迭代次数200,自身学习率为0.3,社会学习率为0.8,最大惯性权重因子为0.8,最小惯性权重因子为0。表1给出了p=1或p=2时,不同样本情形下的反演结果。

表1 不同样本和准则下的断层面参数估计

根据哈佛大学给出的该地震的震源机制解结果,断层面走向为291°,倾角为84°。我们取424个样本值,在范数1的情况下得到的断层面走向为106.4°,倾角为88.9°,两者的走向相差4.6°,倾角相差6.1°(根据Aki等(1980)关于走向和倾角的定义),结果非常接近,说明我们得到的结果是可靠的。

根据粒子群算法求出的参数,可以确定断层面的几何形态。从图1到图4分别给出了4种情形下余震分别在水平面(图1a,2a,3a,4a)、断层面(图1b,2b,3b,4b)和垂直于断层面的横断面(图1c,2c,3c,4c)上的投影,及余震距断层面距离的分布(图1d,2d,3d,4d)。其中,空

图1 取424个样本时l1-范数最小准则结果(a)余震在水平面上的投影;(b)余震在断层面上的投影; (c)余震在垂直于断层面的横断面上的投影;(d)余震距断层面距离的分布

图2 取424个样本时l2-范数最小准则结果(a)余震在水平面上的投影;(b)余震在断层面上的投影; (c)余震在垂直于断层面的横断面上的投影;(d)余震距断层面距离的分布

图3 取82个样本时l1-范数最小准则结果(a)余震在水平面上的投影;(b)余震在断层面上的投影; (c)余震在垂直于断层面的横断面上的投影;(d)余震距断层面距离的分布

图4 取82个样本时l2-范数最小准则结果(a)余震在水平面上的投影;(b)余震在断层面上的投影; (c)余震在垂直于断层面的横断面上的投影;(d)余震距断层面距离的分布

心圆表示余震震源点的位置,黑色线框所包围的区域表示由算法所确定的断层面区域,直线段表示垂直于横断面的断层投影。

将图1和图2、图3、图4进行对比,可以看出,l₁-范数最小准则下求得的断层面受野值点数据的影响较小(图1c,3c),余震距断层面距离的分布呈现宽尾巴的对阵分布(图1d)或者偏态分布(图3d),而l₂-范数最小准则下求得的断层面受野值点数据的影响较大(图2c,4c),余震距断层面距离的分布呈现近似正态分布(图2d,4d)。由于在取l₂-范数时,野值点对应的残差求平方之后“夸大”了它对断层面的作用,因此当数据样本中野值点较多时,使用l₁-范数更为合理。l₁-范数最小准则是一种比较稳健的方法,以前由于它对应的目标函数不可导,在数学上处理起来不方便,而且由于当时计算机没有普及,导致了它很少被使用。现在随着优化算法的发展和计算能力的提高,l₁-范数最小准则的应用比较容易了,而且人们还提出了其他的稳健估计准则,如Guitton等(2003)使用的可以看成l₁和l₂-范数混合的Huber范数最小准则等。

笔者的计算机配置为Pentium IV3.0G 的芯片,512MB内存,计算软件为MATLAB7.0环境,文中所有结果均在此环境下计算得出。

粒子群算法的优势在于算法的简洁性,它易于实现,参数调整少,无需目标函数的导数信息,是一种有效的非线性全局优化算法,具有很强的搜索能力和很高的准确性与稳定性。由于很多地球物理问题都是非线性的,因此粒子群算法具有一定的适用性。尽管粒子群算法仍然是一种随机搜索算法,无法估计参数的精度,但是这种直接搜索的算法也不存在如Gauss-Newton迭代法过程中的初始点选择、病态矩阵处理和误差积累等问题。因此对于许多非线性、多极值的非线性反演问题,采用这种简洁而又易于实现的启发式随机优化算法求解是一种不错的选择。