[摘 ?要]在勘察內(nèi)蒙阿拉善某高放廢物地質(zhì)處置預(yù)選場(chǎng)址目標(biāo)花崗巖體完整性與穩(wěn)定性的地震勘探過(guò)程中，利用基于MSFM（Multistencils Fast Marching Methods）射線追蹤算法的地震初至波走時(shí)反演方法進(jìn)行了近地表高精度速度建模研究。獲得了準(zhǔn)確、精細(xì)的近地表速度模型成果。該成果不但可被用于解決對(duì)于陸上復(fù)雜地表地震勘探來(lái)說(shuō)至關(guān)重要的層析靜校正問(wèn)題，而且可作為地震解釋工作中重要的參考資料。

[關(guān)鍵詞]高放廢物地質(zhì)處置; 復(fù)雜地表; 近地表速度建模; 射線追蹤; 速度反演

中圖分類號(hào)：P5 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1009-914X（2016）23-0013-02

根據(jù)以往的探測(cè)經(jīng)驗(yàn)可知，近地表速度不準(zhǔn)確會(huì)嚴(yán)重影響地下地質(zhì)結(jié)構(gòu)的地震成像效果，對(duì)于地表?xiàng)l件復(fù)雜的工區(qū)更是如此。在勘察內(nèi)蒙阿拉善某高放廢物地質(zhì)處置預(yù)選場(chǎng)址目標(biāo)花崗巖體的過(guò)程中，為了解決對(duì)于陸上工區(qū)地震勘探至關(guān)重要的層析靜校正問(wèn)題，利用基于MSFM（Multistencils Fast Marching Methods， 多模板快速步進(jìn)）射線追蹤算法的地震初至波走時(shí)反演技術(shù)，建立了沿地震測(cè)線準(zhǔn)確、精細(xì)的近地表速度模型。該成果不但被用于有效解決了后續(xù)的層析靜校正問(wèn)題，而且是最終地震地質(zhì)解釋工作中精確、可靠的參考資料。

1 方法原理

復(fù)雜地表問(wèn)題通常指包含劇烈起伏的崎嶇地表和強(qiáng)速度變化的近地表速度結(jié)構(gòu)。在選擇適用于復(fù)雜地表初至波走時(shí)反演的射線追蹤算法時(shí)，通常需要考慮計(jì)算精度、計(jì)算效率、以及算法對(duì)于模型復(fù)雜度的適應(yīng)能力三個(gè)方面[1]。

1.1 MSFM走時(shí)計(jì)算公式

MSFM（Multistencils Fast Marching Methods）多模板快速步進(jìn)算法是由Hassouna等提出的對(duì)于經(jīng)典FMM方法[2]的改進(jìn)算法[3]。其原理是在標(biāo)準(zhǔn)FMM算法的基礎(chǔ)上通過(guò)坐標(biāo)旋轉(zhuǎn)產(chǎn)生多個(gè)FMM計(jì)算模板，使對(duì)角方向鄰點(diǎn)也能夠參與計(jì)算。該算法相對(duì)標(biāo)準(zhǔn)FMM方法在計(jì)算精度方面有顯著提高，并且更容易推廣到處理高維高階問(wèn)題。

其計(jì)算公式如下（為了簡(jiǎn)化討論，僅列舉橫向和縱向空間步長(zhǎng)相等時(shí)的情況）：

① 相鄰模板（S1）：

對(duì)于方向?qū)?shù)的一階近似，可通過(guò)下式求解：

（1）

對(duì)于方向?qū)?shù)的二階近似，可通過(guò)下式求解：

（2）

② 對(duì)角模板（S2）：

對(duì)于方向?qū)?shù)的一階近似，可通過(guò)下式求解：

（3）

對(duì)于方向?qū)?shù)的二階近似，可通過(guò)下式求解：

（4）

在上述公式（1）～公式（4）中，為待求結(jié)果（也就是理論計(jì)算走時(shí)），為地下介質(zhì)空間慢度的函數(shù)，為空間步長(zhǎng)，為鄰域范圍內(nèi)的最小到達(dá)時(shí)。

1.2 算法實(shí)現(xiàn)過(guò)程

算法的實(shí)現(xiàn)需要利用Dijkstra算法思想[4]，即首先將計(jì)算網(wǎng)格點(diǎn)劃分成三個(gè)子集，分別為到達(dá)時(shí)已確定的{已知點(diǎn)子集（known）}，已經(jīng)被計(jì)算但最終到達(dá)時(shí)待定的{窄帶子集（narrow_band）}和尚未被計(jì)算的{遠(yuǎn)點(diǎn)子集（far）}：

（1） 首先將炮點(diǎn)初始化為{已知點(diǎn)子集}，到達(dá)時(shí)為零;并根據(jù)MSFM模板計(jì)算炮點(diǎn)所有在{遠(yuǎn)點(diǎn)子集}中的鄰點(diǎn)的到達(dá)時(shí)，并將這些點(diǎn)移入{窄帶子集}.

（2） 選取{窄帶子集}中走時(shí)最小的網(wǎng)格點(diǎn)移入{已知點(diǎn)子集}，并利用模板計(jì)算該點(diǎn)所有在{窄帶子集}和{遠(yuǎn)點(diǎn)子集}內(nèi)的鄰點(diǎn)的到達(dá)時(shí);若鄰點(diǎn)在{遠(yuǎn)點(diǎn)子集}中，則將該點(diǎn)移入{窄帶子集};若在{窄帶子集}中，則根據(jù)到達(dá)時(shí)大小選擇更新或者保持原到達(dá)時(shí)（取較小者）。

（3） 重復(fù)實(shí)施步驟（2）直至所有網(wǎng)格點(diǎn)移入{已知點(diǎn)子集}。

1.3 Model94模型測(cè)試

為了檢驗(yàn)MSFM射線追蹤算法的計(jì)算精度及其對(duì)于復(fù)雜模型的適應(yīng)能力，本次研究首先設(shè)計(jì)了針對(duì)Model94復(fù)雜地表模型的測(cè)試試驗(yàn)。具體試驗(yàn)流程為：① 給定觀測(cè)系統(tǒng)，利用時(shí)間2階、空間20階交錯(cuò)網(wǎng)格有限差分法[5]求解聲波方程正演出一套地震記錄。② 利用2階MSFM算法計(jì)算在相同觀測(cè)系統(tǒng)下各檢波點(diǎn)位置的波前面到達(dá)時(shí)（初至波到時(shí)）。③ 將MSFM計(jì)算結(jié)果投影到有限差分法地震模擬記錄上，對(duì)比通過(guò)MSFM方法計(jì)算得到的初至波到達(dá)時(shí)（走時(shí)）與模擬地震錄中初至波的吻合程度。

按照上述流程共進(jìn)行了等間隔分布的277炮數(shù)據(jù)的數(shù)值測(cè)試試驗(yàn)。其中位于模型正中的第139號(hào)單炮數(shù)據(jù)數(shù)值測(cè)試結(jié)果如圖1所示。

可以看到，就Model94模型復(fù)雜度的地質(zhì)問(wèn)題來(lái)說(shuō)MSFM算法的計(jì)算精度是可以滿足要求的。圖中顯示的波前面到時(shí)計(jì)算結(jié)果與通過(guò)高階交錯(cuò)網(wǎng)格有限差分法得到的模擬地震記錄中的初至波基本吻合。

1.4 初至波走時(shí)反演速度建模實(shí)驗(yàn)

相對(duì)于全波形反演（FWI）來(lái)說(shuō)，射線類走時(shí)反演算法具有計(jì)算效率高、在應(yīng)用中穩(wěn)健可靠、人工交互技術(shù)成熟等多方面的優(yōu)勢(shì)。上文已經(jīng)證實(shí)了即使是對(duì)于復(fù)雜地表、復(fù)雜地下結(jié)構(gòu)的模型來(lái)說(shuō)，MSFM射線追蹤算法在計(jì)算精度方面也是能夠滿足要求的。因此接下來(lái)，本研究進(jìn)行了后續(xù)的地震初至波走時(shí)反演速度建模實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果如圖2所示。

可以看到，利用基于MSFM射線追蹤算法的初至波走時(shí)反演方法進(jìn)行復(fù)雜地表模型速度反演，即使是從一個(gè)與真實(shí)模型差別很大的初始模型出發(fā)（如圖2a所示），經(jīng)過(guò)多次迭代（圖2b）最終得到的反演結(jié)果依然能夠保證近地表速度分布的精確重建（圖2c）。

2 復(fù)雜地表大型花崗巖體地震勘探中的應(yīng)用實(shí)效

將前述基于MSFM射線追蹤算法的地震初至波走時(shí)反演方法應(yīng)用于內(nèi)蒙阿拉善崎嶇地表工區(qū)高放廢物地質(zhì)處置預(yù)選場(chǎng)址目標(biāo)花崗巖體地震勘探工作中，人工拾取地震初至波（包括直達(dá)波、回轉(zhuǎn)波、折射波或多個(gè)層折射波的組合）走時(shí)，進(jìn)行速度反演建模，所得的建模結(jié)果如圖3和圖4所示。

圖3和圖4分別為內(nèi)蒙阿拉善1線（16.856 km）和2線（10.196 km）近地表速度建模結(jié)果，其中黑色實(shí)線對(duì)應(yīng)射線密度的底界（可以認(rèn)為底界以上部分的速度反演結(jié)果是非常可靠的）。

首先可以看到，本次地震勘探，有效地震射線的分布是非常密集的，因此所得到的兩條測(cè)線近地表速度建模結(jié)果均應(yīng)是可靠的。速度建模成果首先被應(yīng)用于解決了對(duì)于陸上復(fù)雜地表工區(qū)至關(guān)重要的層析靜校正問(wèn)題，進(jìn)而又被作為地震地質(zhì)解釋工作中重要的參考資料。從圖3與圖4中均可觀察到地震波速度在橫向上存在局部異常。在阿拉善1線（圖3）3500 m、4200 m樁號(hào)的下方以及16000 m樁號(hào)附近的兩側(cè)，阿拉善2線（圖4）6000 m、7000 m樁

號(hào)的下方，均顯示出了特征明顯的速度梯度突變。推測(cè)這些異常、突變特征與沿測(cè)線小型斷層、凹陷構(gòu)造的分布有關(guān)，但限于篇幅擬另文予以詳細(xì)闡述。

3 結(jié)論與討論

經(jīng)過(guò)理論分析、模型試算、以及在內(nèi)蒙阿拉善工區(qū)針對(duì)大型目標(biāo)花崗巖體的地震勘探工作實(shí)踐，可以認(rèn)為采用基于MSFM射線追蹤算法的地震初至波走時(shí)反演方法解決復(fù)雜地表?xiàng)l件下近地表高精度速度建模問(wèn)題是可靠而且有效的。

基于MSFM射線追蹤算法的初至波走時(shí)反演可以作為在我國(guó)西部、南部進(jìn)行地震勘探時(shí)，當(dāng)遇到復(fù)雜地表近地表速度建模問(wèn)題時(shí)的針對(duì)性處理技術(shù)和解決措施。

參考文獻(xiàn)

[1] 張釙，劉洪，李幼銘，2000，射線追蹤方法的發(fā)展現(xiàn)狀[J].地球物理學(xué)進(jìn)展，15（1）：36-45.

[2] Sethian J A， Popovici A M， 1999， 3-D travel time computation using the fast marching method[J]. Geophysics， 64（2）： 516-523.

[3] Hassouna M S， Farag A A， 2007， Multi-stencils fast marching methods： a highly accurate solution to the eikonal equation on Cartesian domains[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 29（9）： 1563-1574.

[4] 盧回憶，劉伊克，常旭. 20013， 基于MSFM的復(fù)雜近地表模型走時(shí)計(jì)算[J]. 地球物理學(xué)報(bào)， 56（9）： 3100-3108.

[5] 王秀明， 張海瀾， 2004， 用于具有不規(guī)則起伏自由表面的介質(zhì)中彈性波模擬的有限差分算法. 中國(guó)科學(xué)（ G 輯）， 34（ 5） ： 481-493.

作者簡(jiǎn)介：李穩(wěn)，男，1983年生，工程師，博士，主要從事地震勘探數(shù)據(jù)處理方法研究和地震資料解釋工作.

資金資助：核設(shè)施退役及放射性廢物治理項(xiàng)目（科工二司[2013]727號(hào)）、中國(guó)地震局物探中心青優(yōu)專項(xiàng)（SFGEC2014006）