基金项目:国家自然科学基金(40334040和40474049)和华北水利水电学院高层次人才科研启动项目联合资助.
(1.华北水利水电学院,郑州 450011; 2.中国地震局地球物理勘探中心,郑州 450002; 3. 河南有色金属地质矿产局第三地质大队, 河南南阳 474250)
(1.North China Institate of Water Conservation and Hydroelectric Power,Zhengzhou 450011, Henan, China; 2. Geophysical Exploration Center, CEA, Zhengzhou 450002, Henan, China)(3. The 3rd Team of Henan Nonferrous Metal Geologic Minerals Bureau, Nanyang 474250, Henan, China)
multiscale successive approximation; tomography of seismic traveltime; SAGA; A'nyemaqen suture zone
备注
基金项目:国家自然科学基金(40334040和40474049)和华北水利水电学院高层次人才科研启动项目联合资助.
根据多尺度逐次逼近思想,建立了多尺度逐次逼近退火遗传算法。该方法能有效地解决遗传算法中存在的收敛“早熟”问题。采用该方法对一个高速异常体进行了数值模拟试验,并对青藏高原东北缘阿尼玛卿缝合带东段上部地壳速度结构实际资料进行了处理,结果表明,多尺度逐次逼近退火遗传算法能够较好地应用于地震走时层析成像研究。
SAGA is a combined algorithm of Simulated Annealing(SA)and Genetic Algorithm(GA). The method has the feature of higher searching efficiency than SA or GA alone, but can not radically solve the problem of convergence prematurity of GA. The multiscale SAGA algorithm, which is based on multiscale successive approximation, can effectively solve the problem. We used the method to reconstruct an anomalous high-velocity body in a numerical test model, and to study the upper crustal structure characteristics of the eastern segment of the A'nyemaqen suture zone. Both test and practical application prove that the multiscale SAGA is a good tool for seismic tomography.
引言
地震走时层析成像是一种获取地壳速度结构的重要手段。地震走时和地壳速度之间存在着强非线性关系,对这类问题的早期解决方法主要是局部的线性及线性化反演方法(Bishop等,1985; Scales,1987; White,1989; Zhang,1989; Kosloff等,1996)。这些方法一般需要计算目标函数的一阶或二阶导数,常用的有广义逆方法、奇异值分解法、最速下降法、拟牛顿法、马奎特法、共轭梯度法和广义最小二乘法等。例如,Lutter等(1990)将共轭梯度法应用到PASSCAL Ouachita 试验的折射和反射走时数据上; Zelt(1992)利用阻尼最小二乘法得到一种同时获得速度和界面的走时反演方法; Hole(1992)利用反投影法发展一种层析成像方法。然而,这些局部反演方法不能够完全适应地震走时层析成像的需要,因为对于需要做线性化处理的方法来讲,其结果是不精确的,对于需要解大型方程组的方法,其解的可靠性有待讨论,对于需要求导的方法,若导数不存在则无法应用。尽管局部反演方法计算速度快,可以反演模型空间几百个以上的参数,但这些反演方法只利用局部有限的信息来改进初始模型,因此对初始模型有很大的依赖性,容易陷入解空间的局部极值区而得不到全局极值。
非线性全局优化反演技术在地球物理领域的应用始于20世纪80年代中期,全局方法一般分为枚举算法和随机算法。随机算法的典型代表是蒙特卡洛(Monte Carlo)方法(Jin和Madariaga,1994)、模拟退火算法(SA)(Kirkpatrick,1983; Chen,2005)和遗传算法(GA)(Jin和Madariaga,1993)。Rothman(1986)首先提出利用模拟退火方法来解决反演难度较大的剩余静校正问题,Jin 等(1994)利用两步蒙特卡洛方法反演地壳速度,Stoffa(1991)最先研究利用遗传算法进行地震波形的多参数非线性反演问题。总的说来,全局反演方法简单、有效,具有普遍适应性,不存在局部反演方法所具有的诸多问题。但是全局随机方法计算十分耗时,目前仅仅局限于反演简单的速度模型。相比较而言,模拟退火算法和遗传算法的计算效率都比蒙特卡洛方法高,并且计算与初始值的选择关系不大,不需对问题作线性化处理。两者比较,遗传算法的收敛速度更快。Sen等人(1997)将遗传算法和模拟退火方法相结合得到称作退火遗传算法的方法,崔建文等(2004)也有类似的做法。但遗传算法的“早熟”问题依旧没有得到解决。本文利用退火遗传算法,采用一种多尺度逐次逼近的方法(Bunks等,1995; Jin等,2000; 师学明等,2000)来反演地壳速度结构,一定程度上抑制了遗传算法中“早熟”问题的发生。
1 基本原理
遗传算法的基础是自然选择和自然遗传机制,其基本思想是基于模仿生物界的遗传过程,把问题的参数用基因代表,把问题的解用染色体代表(在计算机中为字符串),从而得到一个由具有不同染色体的个体组成的群体,这个群体在问题特定的环境里生存和竞争,适者有最好的机会生存和产生后代。后代随机地继承了父代的最好特征,并在生存环境的控制和支配下继续这一过程。这样群体的染色体都将逐渐适应环境,不断进化,最后收敛到一族最适应环境的类似个体,即得到问题最优的解。
遗传算法的一个主要缺点是收敛“早熟”问题,这实际上是问题的某种局部优化解。要继续进化,就要在群体中引入新个体,在一定程度上恢复基因的多样性。依据这一想法,根据一定的“早熟”判断标准,如果出现“早熟”现象,则进入模拟退火操作,生成新的个体,替代群体中相对较差的个体,从而在一定程度上解决“早熟”问题。
1.1 遗传算法操作设群体由M个模型Xj(j=1,2,…,M)构成,其初始值由下式随机生成
(Xj)i=(Xmin)i+ξ(Xmax-Xmin)i.(1)
式中,Xmax和Xmin分别为模型参数的最大值和最小值,ξ是在[0,1]区间上均匀分布的随机变量,i表示模型中的第i个参数。从各个目标函数计算其生存系数,确定生存系数最大值和平均值,计算生存概率:
fj=1/(sigma)exp(-Φ(Xj)/T),(2)
sigma=∑Mi=1exp(-Φ(Xj)/T).(3)
式中,fj是个体的生存系数,T是引入模拟退火算法中的温度,它是随计算进程逐渐减小的量,设降温次数为k,温度Tk=0.99kT0,其中T0为初始温度。
平均生存系数由下式定义:
f^-=∑Mj=1fj/M.(4)
收敛判据由下式给出:
min{Φ(Xj)}≤ε1.(5)
式中,Φ(Xj)为第j个模型的目标函数,ε1是给定的某一较小的常数。若无模型满足判据,则用M个模型建造新的子代群体,经历遗传算法的三种核心运算:
(1)“繁殖”
用轮赌法进行选择操作,生成父代种群,令
Ai=∑ij=1fj.(6)
当公式(1)中的ξ满足Ai-1<ξ≤Ai时,将第i个个体选入种群,作M次选择,可产生规模为M的父代种群。种群中的个体可以重复,生存系数较大的个体有可能被多次选中,在种群中形成优势群体; 相反,生存系数较小的个体将以较大的机率被淘汰,这体现了生物界优胜劣汰的规律。为了不丧失已获得的进化成果,采用杰出个体保护策略,用群体中的最优个体替换种群中最差的个体。
(2)“交配”
首先采用等长度二进制编码,将模型中的参数用二进制表示,并将它们依序相连形成字串,其长度用L表示。对群体中的个体随机配对,在一对个体中,随机确定两个交叉位置L1、L2,且L1=ξ1·L小于L2=ξ2·L,然后交换两个个体位于L1~L2之间的代码。
(3)“变异”
对“交配”产生的新个体,随机确定变异位置L3=ξ3·L,按变异概率Pm>ξ4进行变异操作,将L3位置上的码变异的值,由0变为1或由1变为0。
“早熟”判别依据下式:
fmax-fmin≤η.(7)
式中,fmax和fmin分别是群体中所有个体的生存系数的最大、最小值,η为给定常数。如果判据成立,则重新生成新的个体,计算各自的生存系数,并与旧的个体相比较,如果新生成的个体的生存系数较大,那么就依照一定的更新概率来决定是否替代旧的个体。
1.2 模拟退火算法操作采用Metropolise算法,对于适合度满足fj<f^-的模型Xj,随机生成一新模型
Xjg=Xj+ΔXj.(8)
式中ΔXj为模型参数的随机小扰动向量。
计算目标函数Φig=Φ(Xjg),令ΔΦ=Φig-Φj,如果ΔΦ<0,则用Xjg替换Xj; 反之,则计算接受概率p=exp[-ΔΦ/Tk]。当p≥ξ时,用Xjg替换Xj,否则重新生成新模型。然后计算新个体的适合度fj及平均适合度f^-,令k=k+1,最终温度Tk=0.99kT0,其中T0为初始温度,k为降温次数。
逐次逼近的思想是为解决某一个问题,而从一个与该问题的实质内容有着本质联系的较大范围开始进行解决,逐步缩小范围,逐步逼近,最后达到所要求的解。根据多尺度逐次逼近的思想,首先把速度场划分为不同的空间尺度,定义网格节点上的速度作为待反演参数,采用双三次样条函数模型参数化,正问题采用一种快速、准确的有限差分走时计算方法(Zhao,1996),反问题采用多尺度逐次逼近退火遗传算法。首先在较大尺度空间范围内进行全局寻优,将得到的结果作为下一级较小空间尺度全局寻优的初始值,直至满足收敛判别标准。这是因为反问题的目标函数极值点的分布随空间尺度的变化而发生变化,在大尺度上较少而且相互分得很开,因而在大尺度上较易于求解。同时,由于待反演的模型个数减少,“染色体”的长度变短,“基因交换”、“变异”等“遗传”操作更为彻底,有效地克服了基因丢失引起的“早熟”收敛问题。将大尺度的解当作次一级尺度反问题的初始模型集,所增加节点上的速度值经双三次样条插值得到,然后再进行退火遗传反演,如此类推逐次逼近全局最优解,这便是多尺度逐次逼近退火遗传算法的核心。据此,可建立地震走时层析成像流程(图1)。
2 数值模拟试验
图2是一个2维速度模型,沿x方向长30 km,沿z方向深8 km, 地表速度v0=2 km/s,沿z方向速度呈梯度变化,dv/dz=0.5(1/s),在x方向12~18 km和z方向1~7 km的区域内制作一个速度异常区,其速度值是原来数值增加1 km/s。将激发点和接收点均放置在地表。炮点分别位于0、5、10、15 、20、25、30 km处,共7炮。接收点分别位于0、1、2、3、…、30 km处,共31道,对于每炮激发,31道都全部参与接收。正演走时计算采用规则网格1 km×1 km。图2b为其射线分布图。在经历了9、15、25、55、99个节点共5次反演后,对应网格分别为15 km×4 km、7.5 km×4 km、7.5 km×2 km、3 km×2 km和3 km×1 km。对于初始模型的9个节点,其水平节点位置为0、
15、30 km,其垂向节点位于0、4、8 km深度上。给定9个节点处的速度搜索范围均为1.0~7.0 km/s。遗传算法的参数采用交换概率Pc=0.90,变异概率Pm=0.05,模型群体大小Q=16,最大迭代次数为500, 随机种子数为23。图3为不同节点数目的反演结果,随着节点数的增加,高速异常区域越来越得到良好的反映,当节点数大于55以后,反演结果趋于相同。 图4a表明反演速度节点数为55时虚线以上部分的结果是可信的。图4b、c为检测板试验图,随着深度的增加,分辨逐渐变差,速度模型的两端分辨差,中间分辨较好,这是与射线分布有关的。图4d为走时残差随反演节点数的变化,节点9的初始残差是0.748 s,节点55的走时残差减小为0.022 s。之后,尽管速度节点数目增加,走时残差却不再减小。图5是采用多尺度逐次逼近退火遗传算法(a)和直接采用退火遗传算法对279个格点速度反演(b)的结果对比,前者能够反映高速异常区域,而后者经历500代反演,走时残差从初始的0.576 s下降到0.131 s,之后不再减小,此时不能够正确反映高速异常区。分析其原因,后者可能是由于收敛早熟问题而陷入局部极值区所造成的。
3 实际应用
为研究阿尼玛卿缝合带东段及其两侧邻近区域的上部地壳结构特征,2004年5~6月,中国地震局地球物理勘探中心在实施马尔康—古浪约600 km长的宽角折射、反射人工地震测深剖面的同时,在川北甘南交界地区完成了一条长为210 km,走向近南北、与宽角折射、反射人工地震测深剖面完全重合的高分辨折射地震探测剖面(图6),在沿剖面共进行了10次爆破,获得了877个地震初至走时数据,利用前文提到的多尺度逐次逼近退火遗传算法,进行节点为9、27、85、153、297、594、1 166和2 226等一系列的反演,对于初始模型的9个节点,其水平节点位置为190、295、410 km,其垂向节点位于0、10、20 km深度上,0 km深度上节点的速度搜索范围为1.0~8.0,10 km深度上节点的速度搜索范围为4.0~8.0 km/s, 20 km深度上节点的速度搜索范围为6.0~8.0 km/s。遗传算法的参数采用交换概率Pc=0.90,变异概率Pm=0.05,模型群体大小Q=32,随机种子数为23,最大迭代次数为500。
图7为不同节点数目的反演结果,随着节点数的增加,反演出来的速度模型的空间尺度越来越小,直到不变化。
(资料源自马杏垣(1989))F1:库赛湖—玛沁断裂; F2:武都—迭部断裂;
F3:舟曲—两当断裂; F4:阿坝弧形断裂
图8a表示反演速度节点数297时虚线以上部分的结果是可信的,图8b、c为检测板试验图,显示扰动恢复的情况,虚线以上部分得到较好恢复,表明该部分反演结果是可信的。图8e表示走时残差随反演节点数的变化,初始残差是0.867 s,节点增大到297以后,走时残差减小为0.063 s,之后,尽管反演的速度节点增加,但是走时残差却不再减小。
分析认为结晶基底的速度大约是6.0 km/s,在测线南端(190 km桩号),基底深度约3.5 km; 在215~240 km桩号之间,结晶基底变深,最深达6.0 km; 在240~260 km桩号之间,结晶基底深度逐渐变浅至3.0 km; 在阿尼玛卿缝合带下,速度横向变化剧烈,260~280 km桩号之间呈现一个低速带,最大速度值仅为5.4~5.6 km/s,这可能由于强烈的挤压,此处基底已经破碎; 280~290 km桩号之间,速度值突然抬升,说明此处基底完整性较好; 埋深大约4.5 km; 290~300 km桩号之间又出现一个低速条带,但是其规模没有前者大; 300~310 km桩号之间速度急剧升高,基底深度上抬至3 km; 在缝合带内,基底深度整体呈南深北浅的趋势,横向速度和基底深度均发生剧烈变化,缝合带内的基底已经不再完整; 310 km桩号以北的西秦岭褶皱带内,依速度特征大体分为两段,340 km桩号为两者分界线,桩位以南为高速段,以北为低速段; 基底深度在350 km桩号处最浅,约2.0 km,之后自南向北,基底埋深逐渐加深。在研究区有数条断裂通过剖面,在283 km桩号附近,速度从低速剧变为高速,基底深度由深突然变浅,这里正好是库赛湖—玛沁断裂穿过剖面的位置; 从图8a中所示的速度低、高相间的状况来看,该断裂可能是由多条断层组成的断裂带; 在320~330 km桩号之间,速度横向变化亦较大,基底深度从2.2 km突然加深至4.5 km,此处是武都—迭部断裂(F2)的体现; 340 km桩号是低速与高速的分界线,是舟曲—两当断裂(F3)的反映。可见,在基底上速度结构的横向剧烈变化和基底深度的突变都是由于断层的存在引起的,因此我们可以把速度横向剧烈变化和基底深度突变看作是断层存在的标志。
4 结论和讨论
数值模拟试验和实际资料的处理结果表明,多尺度逐次逼近退火遗传算法是反演地壳速度结构的一种有效的非线性反演方法。用它可以较快地逼近全局最优解,避免选取初始模型、计算Jacobi偏导数矩阵和陷入局部极值等缺点。该方法克服了蒙特卡洛方法和模拟退火反演方法中的模型搜索时间较长的弱点,是解决退火遗传算法“早熟”问题的一条有效途径。而且,反演方法本身对观测系统的布设没有要求,可以适用于任意复杂的观测系统。
多尺度逐次逼近退火遗传算法不能给出分辨矩阵,但是,对反演结果的分辨可依靠检测板、误差曲线、射线分布等手段来估计。同线性或线性化反演方法相比较,全局反演方法的收敛速度仍显缓慢,若能将多尺度逐次逼近退火遗传算法同局部反演方法(如单纯形)相结合,则既可保证全局性又具有一定的收敛速度。
王夫运博士、段永红博士、徐朝繁博士、刘志高级工程师和杨卓欣高级工程师等提供了宝贵的意见,表示感谢。
- 崔建文.2004.一种改进的全局优化算法及其在面波频散曲线反演中的应用[J].地球物理学报,47(3):521-527.
- 师学明,王家映,张胜业,等.2000.多尺度逐次逼近遗传算法反演大地电磁资料[J].地球物理学报,43(2):122-130.
- Bishop T, Bube K, Cutler R, et al.1985.Tomographic determination of velocity and depth in laterally varying media[J]. Geophysics, 50, 903-923.
- Bunks C, Salek F E M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion[J]. Geophysics. 60(5):1457-1473.
- Hole J A, Clowes R M, Ellis R M. 1992. Interface inversion using broadside seismic refraction data and three-dimensional travel time calculations[J]. JGR, 97, 3417-3429.
- Improta L, Zollo A, Herrero A, et al. 2002. Seismic imaging of complex structures by non-linear traveltime inversion of dense wide-angle data: application to a thrust belt[J]. Geophys.J.Int. 151, 264-278.
- Jin S, Beydoun W. 2000. 2D multiscale non-linear velocity inversion[J]. Geophysical Prospecting.48, 163-180.
- Jin S, Madariaga R. 1994. Non-linear velocity inversion by a two step Monte-Carlo method[J]. Geophysics 59, 577-590.
- Jin S, Madariaga R. 1993. Background velocity inversion with a genetic algorithm[J]. Geophysical Research Letters, 20, 93-96.
- Kirkpatrick S, Gelatt C D Jr, Vecchi M P. 1983. Optimization by simulated annealing[J]. Science, 220, 671-680.
- Kosloff D, Sherwood J, Koren Z, et al, 1996. Velocity and interface depth determination by tomography of depth migrated gathers[J]. Geophysics, 61, 1511-1523.
- Lutter W J, Nowack R L. 1990. Inversion for crustal structure using reflections from the PASSCAL Ouachita Experiment[J]. JGR, 95, 4633-4646.
- Rothman D H.1986.Automatic estimation of large residual static correction,Geophysics,51(2):322-346.
- Scales J A. 1987. Topographic inversion via the conjugate gradient method[J]. Geophysics, 52,179-185.
- Sen M K, Stoffa P L.1997. Global optimization methods in geophysical inversion[J]. ELSEVIER Science, New York: 140-144.
- Stoffa P L,Sen M K.1991.Nonlinear multiparameter optimization using genetic algorithms:inversion of plane-wave seismogram.Geophysics,56(11):4794-1810.
- White D J. 1989. Two-dimensional seismic refraction tomography[J]. Geophys. J. Int, 97, 223-245.
- Zelt C A, Smith R B. 1992. Seismic traveltime inversion for 2-D crustal velocity structure[J]. Geophys.J.Int.108,16-34.
- Zhang X K.1989 Generalized inversion and its application in inverse scattering problems[J]. Journal of Computational Mathematics, 7(4):374-382.
- Zhao P. 1996. An efficient computer program for wavefront calculation by the finite difference method[J]. Computers & Geosciences, 22(3):239-251.