三維復(fù)雜黃土塬區(qū)變加密不等距網(wǎng)格初至走時計(jì)算與特征分析

韓令賀， 孫章慶， 胡自多， 劉威， 韓復(fù)興， 田彥燦

1 中國石油勘探開發(fā)研究院西北分院， 蘭州 730020 2 中國石油大學(xué)(北京)地球物理學(xué)院， 北京 102249 3 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026

0 引言

復(fù)雜黃土塬區(qū)在我國有較大的分布范圍，尤其是在我國西北部的鄂爾多斯盆地和塔里木盆地西南部.這兩個區(qū)域均蘊(yùn)藏著豐富的油氣資源.但是，由于經(jīng)歷了長期的風(fēng)化、侵蝕、沖刷、切割，黃土塬區(qū)形成了塬、梁、坡、溝等復(fù)雜地貌，地表高程變化劇烈，黃土層厚度劇烈變化且黃土疏松、彈性差、速度低，表層結(jié)構(gòu)復(fù)雜多變，因此黃土塬區(qū)一直被視為地震勘探的“禁區(qū)”(閻世信和呂其鵬，2002).這些復(fù)雜的地震地質(zhì)條件使得復(fù)雜黃土塬區(qū)也面臨著地震激發(fā)接收困難、地震資料信噪比低、噪聲異常發(fā)育、相干和次生干擾嚴(yán)重、近地表建模困難、靜校正問題嚴(yán)重、中深層成像困難或基本不成像等諸多挑戰(zhàn)，因此對復(fù)雜黃土塬區(qū)地震勘探的相關(guān)問題進(jìn)行研究具有重要的理論和實(shí)際意義(劉寶國等，2017).

黃土塬區(qū)的復(fù)雜地震地質(zhì)條件通常會引起地震波場結(jié)構(gòu)異常復(fù)雜，如：初至扭曲嚴(yán)重、面波散射發(fā)育、體波地形散射發(fā)育、黃土層內(nèi)吸收衰減嚴(yán)重、黃土層內(nèi)多次波發(fā)育、中深層反射能量很弱且非常不連續(xù)等(閻世信和呂其鵬，2002；劉寶國等，2017).對這些復(fù)雜波場結(jié)構(gòu)進(jìn)行深入解析是解決上述復(fù)雜黃土塬區(qū)諸多挑戰(zhàn)的基礎(chǔ).然而，目前專門針對復(fù)雜黃土塬區(qū)地震波傳播規(guī)律與低信噪比數(shù)據(jù)的形成機(jī)制的相關(guān)研究并不多見，這幾乎也是一直以來制約復(fù)雜黃土塬區(qū)地震數(shù)據(jù)處理效果欠佳的一個瓶頸(閻世信和呂其鵬，2002；劉寶國等，2017).面對該難題，很難一蹴而就的全面解析該區(qū)域的復(fù)雜地震波場結(jié)構(gòu).不過，在它們中初至波是相對最重要、最易、最應(yīng)先解析的地震波場成分，因?yàn)樗鼘乇硭俣冉：挽o校正，這兩項(xiàng)位于相對處理前端又極其重要的地震數(shù)據(jù)處理環(huán)節(jié)具有重要的意義.鑒于此，本文將緊密圍繞“三維復(fù)雜黃土塬分區(qū)初至走時計(jì)算與特征分析”這一主題展開，意在從地震波場運(yùn)動學(xué)角度，分析復(fù)雜黃土塬區(qū)初至的基本特征，并為后續(xù)該區(qū)近地表建模和靜校正等建立穩(wěn)定、靈活且兼顧精度和效率的初至走時算法.

綜上所述，為了精細(xì)刻畫三維復(fù)雜黃土塬的復(fù)雜近地表介質(zhì)特征，建立穩(wěn)定、高效、高精度、靈活的初至波走時算法，并深入分析該介質(zhì)條件下初至波的各方面特征，本文提出一種三維復(fù)雜黃土塬分區(qū)變加密不等距網(wǎng)格初至波走時計(jì)算方法，該方法采用三維分區(qū)局部變加密不等距網(wǎng)格精細(xì)剖分復(fù)雜黃土塬模型，采用一種三維迎風(fēng)變網(wǎng)格雙線性插值法精細(xì)計(jì)算初至波走時，最后基于此得出三維復(fù)雜黃土塬區(qū)初至波的一些運(yùn)動學(xué)特征.下面將緊密圍繞這幾個方面的內(nèi)容展開.

1 三維復(fù)雜黃土塬模型的網(wǎng)格剖分

三維黃土塬模型通常非常復(fù)雜，采用怎樣的網(wǎng)格剖分才能精細(xì)地刻畫該模型，并便于后續(xù)初至波走時計(jì)算非常重要.為了同時兼顧計(jì)算精度和效率，本文提出一種三維分區(qū)局部變加密不等距網(wǎng)格剖分三維復(fù)雜黃土塬區(qū)模型，具體步驟如下：①讀入三維地表高程數(shù)據(jù)、計(jì)算空間范圍、計(jì)算網(wǎng)格間距等信息；②根據(jù)步驟①讀入的信息，在計(jì)算范圍內(nèi)，采用網(wǎng)格間距較大的均勻的正方體粗網(wǎng)格剖分整個計(jì)算空間；③對地表不做任何近似，讓它根據(jù)高程直接落在步驟②的均勻的正方體網(wǎng)格線上，并去掉地表以上非必要網(wǎng)格，形成地表附近的不等距網(wǎng)格；④從地表處所在的網(wǎng)格單元開始向地下由淺至深，按照(2n,2n-1,2n-2,…,2)指數(shù)衰減的加密程度對步驟②的單個正方體粗網(wǎng)格，進(jìn)行網(wǎng)格加密,其中n為加密因子，它用于控制網(wǎng)格的加密程度；⑤以初至走時計(jì)算時的震源點(diǎn)為中心點(diǎn)由內(nèi)向外，按照(2n,2n-1,2n-2,…,2)指數(shù)衰減的加密程度，對步驟④的網(wǎng)格，實(shí)施震源附近的網(wǎng)格加密.

在上述復(fù)雜網(wǎng)格生成的過程中，步驟④和⑤最為關(guān)鍵，它們均涉及加密因子n的選取和加密范圍的確定兩個核心問題.其中，在實(shí)施步驟④的地表附近的網(wǎng)格加密時，加密因子n的選取是由最粗的網(wǎng)格間距的大小、地表處需要加密的程度和計(jì)算精度共同決定的，n越大地表附近加密程度越大，計(jì)算精度越高；關(guān)于加密范圍，由于采用的是變加密策略，網(wǎng)格的加密程度是漸變的，數(shù)值實(shí)驗(yàn)表明：由地表向下，上一級加密程度的網(wǎng)格的加密范圍不超過兩個下一級網(wǎng)格的網(wǎng)格間距，即可達(dá)到加密效果和數(shù)值計(jì)算的穩(wěn)定，這也能最大程度的保證計(jì)算效率.此外，在實(shí)施步驟⑤的震源附近的網(wǎng)格加密時，加密因子n的選取是由最粗的網(wǎng)格間距的大小、震源附近需要加密的程度和計(jì)算精度共同決定的；而加密范圍則是根據(jù)震源附近需要控制的計(jì)算誤差來確定的，較低誤差要求的范圍越大則加密的范圍越大，達(dá)到上一級加密程度的誤差控制范圍后，再過渡到下一級加密網(wǎng)格，直到過渡到常規(guī)均勻粗網(wǎng)格.

以上實(shí)現(xiàn)步驟最終獲得了如圖1所示的精細(xì)剖分三維復(fù)雜黃土塬區(qū)的分區(qū)局部變加密不等距網(wǎng)格.與其他網(wǎng)格相比，此處采用的網(wǎng)格有如下優(yōu)缺點(diǎn)：①地表附近的不等距網(wǎng)格，不對地表高程做任何近似，完全準(zhǔn)確的保留地表高程點(diǎn)的實(shí)際位置；②從地表開始向地下由淺入深的變加密網(wǎng)格，能更加細(xì)致的描述復(fù)雜黃土塬的地表形態(tài)和近地表復(fù)雜介質(zhì)分布(尤其是薄風(fēng)化層)，同時還有利于保證復(fù)雜地表處初至波走時計(jì)算的穩(wěn)定性和精度；當(dāng)然，該加密網(wǎng)格一定程度上會增加計(jì)算時的存儲空間和計(jì)算量，但是由于實(shí)施的僅是地表附近的變加密，故僅會增加非常有限的存儲空間和計(jì)算量；③由震源處出發(fā)由內(nèi)而外的變加密網(wǎng)格，能夠在不過多增加計(jì)算量的前提下，很好的解決常規(guī)初至波走時算法在震源附近誤差較大的問題；當(dāng)然，這種誤差的控制是通過大幅度減少網(wǎng)格間距來實(shí)現(xiàn)震源附近截?cái)嗾`差降低的，雖然它不能完完全全從本質(zhì)上解決震源奇異性的問題，但這種變加密策略是一種僅花費(fèi)很少計(jì)算量的增加來大幅度提高計(jì)算精度的高性價比算法；④與曲線網(wǎng)格相比，上述網(wǎng)格剖分方法無需單獨(dú)花費(fèi)網(wǎng)格生成的計(jì)算量而實(shí)現(xiàn)對不規(guī)則邊界不做任何近似處理；⑤在編程實(shí)現(xiàn)時，只要合理的處理好各個加密級別之間的調(diào)用和重新初始化，僅需在地表附近單獨(dú)開辟用于存儲加密節(jié)點(diǎn)的數(shù)組空間，而在震源加密時也僅需在一個數(shù)組中來回采樣和反復(fù)初始化賦值多次，即可完成不同加密程度的變加密運(yùn)算；⑥因?yàn)椴捎镁植考用?，所以算法還能保證后續(xù)初至波走時計(jì)算的絕大部分計(jì)算工作量在規(guī)則均勻的正方體網(wǎng)格中進(jìn)行.

圖1 三維復(fù)雜黃土塬模型的網(wǎng)格剖分Fig.1 Mesh generation of 3D complex loess plateau model

2 面向三維復(fù)雜黃土塬區(qū)的初至走時計(jì)算

上述采用了三維分區(qū)局部變加密不等距網(wǎng)格剖分三維復(fù)雜黃土塬模型，在該復(fù)雜網(wǎng)格中實(shí)現(xiàn)面向三維復(fù)雜黃土塬區(qū)的初至波走時計(jì)算，核心需要解決復(fù)雜網(wǎng)格中的局部走時算法和復(fù)雜網(wǎng)格中的整體波前擴(kuò)展方式這兩個問題，下面將詳細(xì)闡述這兩方面問題.

2.1 復(fù)雜網(wǎng)格中的局部走時算法

如圖2a所示，圖1中的三維分區(qū)變加密不等距網(wǎng)格，實(shí)際上包括了如圖2b—e所示的四種由簡單到復(fù)雜的情況：①網(wǎng)格間距沒有變化的均勻立方體網(wǎng)格(圖2b)；②網(wǎng)格間距有變化的不均勻過渡區(qū)立方體網(wǎng)格(圖2c)；③鄰近地表附近的復(fù)雜網(wǎng)格(圖2d)；④地表處的復(fù)雜網(wǎng)格(圖2e).這4種不同的網(wǎng)格情況，需要采用不同的局部走時計(jì)算公式.目前，關(guān)于這方面的三維走時計(jì)算問題，有學(xué)者曾提出過一種迎風(fēng)混合法(孫章慶等，2017)，該方法采用不等距差分公式來簡潔的處理地表附近的走時計(jì)算問題，但是相比于統(tǒng)一的迎風(fēng)雙線性插值法(孫章慶等，2015)，該方法在適應(yīng)上述變加密不等距網(wǎng)格時相對困難，且地表附近的有限差分法精度相對有限.針對該問題，在此提出適應(yīng)上述復(fù)雜網(wǎng)格的迎風(fēng)變網(wǎng)格雙線性插值法.該方法的核心思想是：首先基于Fermat原理，通過采用無條件穩(wěn)定的迎風(fēng)思想，在封閉被算點(diǎn)的鄰近所有網(wǎng)格單元中，挑選走時分布最小的網(wǎng)格單元，然后再在該網(wǎng)格單元上采用雙線性插值法進(jìn)行局部走時計(jì)算公式的建立.在具體局部走時計(jì)算公式建立時，迎風(fēng)變網(wǎng)格雙線性插值法可以將上述提及的4種情況簡化為如圖2f—h所示的3種情況：

(1)非地表附近的規(guī)則立方體網(wǎng)格(圖2f)中的局部走時計(jì)算公式

如圖2所示，圖2b、c對應(yīng)的局部網(wǎng)格形態(tài)最終都有可能通過基于Fermat原理的迎風(fēng)思想的挑選，簡化為如圖2f所示的規(guī)則立方體網(wǎng)格.網(wǎng)格中G點(diǎn)或N點(diǎn)的走時TG or N為需要計(jì)算的走時值.陰影部分網(wǎng)格單元MACB為基于Fermat原理的迎風(fēng)思想，從圖2b、c中封閉G點(diǎn)或N點(diǎn)的24個正方形網(wǎng)格單元中挑選出的，走時分布為最小的網(wǎng)格單元，其端點(diǎn)M點(diǎn)、A點(diǎn)、C點(diǎn)、B點(diǎn)的走時值TM、TA、TC、TB為已知.在這種情況下，可以直接基于雙線性插值法(孫章慶等，2015)建立局部走時計(jì)算公式：

(1)

其中，s為當(dāng)前計(jì)算網(wǎng)格單元的平均慢度(速度的倒數(shù))，h為正方體網(wǎng)格的網(wǎng)格間距.公式(1)給出了在一個網(wǎng)格單元，其網(wǎng)格端點(diǎn)走時值均為已知且網(wǎng)格單元內(nèi)走時為線性分布假設(shè)的情況下，計(jì)算其鄰近點(diǎn)走時值的局部計(jì)算公式.但是，與常規(guī)的雙線性插值法(孫章慶等，2015)不同，為了無條件穩(wěn)定的適應(yīng)如圖1所示的復(fù)雜變化的網(wǎng)格，在此提出相應(yīng)的迎風(fēng)變網(wǎng)格雙線性插值法，其核心是在復(fù)雜變化的網(wǎng)格中如何確定公式(1)建立時的已知條件，即如何確定如圖2f所示的網(wǎng)格單元MACB及其中的端點(diǎn)M點(diǎn)、A點(diǎn)、C點(diǎn)、B點(diǎn)的走時值TM、TA、TC、TB.確定圖2f所示的網(wǎng)格單元實(shí)際上涉及兩種情況：①如圖2b所示，當(dāng)G點(diǎn)位于均勻立方體網(wǎng)格中時，根據(jù)Fermat原理M點(diǎn)應(yīng)為與G點(diǎn)直接相鄰的6個網(wǎng)格節(jié)點(diǎn)中走時值最小的網(wǎng)格節(jié)點(diǎn)，A點(diǎn)和B點(diǎn)則分別為與M點(diǎn)直接相鄰的另外兩個坐標(biāo)軸方向上走時值較小者，即TM、TA、TB需滿足：

圖2 迎風(fēng)變網(wǎng)格雙線性插值法(a) 三維復(fù)雜地表附近的局部網(wǎng)格； (b—e) 不同位置的網(wǎng)格形態(tài)；(f—h)簡化后的計(jì)算網(wǎng)格.Fig.2 The method of upwind variable grid bilinear interpolation(a) Local grid nearby the 3D complex earth′s surface; (b—e) The shapes of grid at different locations; (f—h) Simplified computational grid.

TM=min(TMi)i=1,2,…,6,

(2a)

TA=min(TMix+,TMix-),TB=min(TMiy+,TMiy-),i=1或6,

(2b)

TA=min(TMiy+,TMiy-),TB=min(TMiz+,TMiz-),i=2或4,

(2c)

TA=min(TMix+,TMix-),TB=min(TMiz+,TMiz-),i=3或5.

(2d)

公式(2a)表述了M點(diǎn)為與G點(diǎn)直接相鄰的6個網(wǎng)格節(jié)點(diǎn)中走時值最小的網(wǎng)格節(jié)點(diǎn)的含義，而公式(2b)—(2d)則表述了A點(diǎn)和B點(diǎn)分別為與M點(diǎn)直接相鄰的另外兩個坐標(biāo)軸方向上走時值較小者(例如：如圖2b所示，若M=M6時，TA和TB則分別可以表示為：TA=min(TM6x+,TM6x-),TB=min(TM6y+,TM6y-)，其中：TM6x+,TM6x-分別為x軸方向上M6鄰近的點(diǎn)M6x+和M6x-的走時值，TM6y+,TM6y-則為y軸方向上M6鄰近的點(diǎn)M6y+和M6y-的走時值)；②如圖2c所示，當(dāng)N點(diǎn)位于網(wǎng)格間距有變化的不均勻立方體網(wǎng)格時，同樣根據(jù)Fermat原理來確定網(wǎng)格單元MACB及其中的端點(diǎn)M點(diǎn)、A點(diǎn)、C點(diǎn)、B點(diǎn)的走時值TM、TA、TC、TB，但是此時情況變得復(fù)雜了很多，其中M點(diǎn)應(yīng)為與N點(diǎn)直接相鄰的10個網(wǎng)格節(jié)點(diǎn)中走時值最小的網(wǎng)格節(jié)點(diǎn)，而A點(diǎn)和B點(diǎn)的確定方法已有所變化，TM、TA、TB需滿足：

TM=min(TMi)i=1,2,…,10,

(3a)

TA=min(TMix+,TMix-),TB=min(TMiy+,TMiy-),i=1或10,

(3b)

TA=min(TMiy+,TMiy-),TB=TM1x+orTM1x-,i=2或4,

(3c)

TA=min(TMix+,TMix-),TB=TM1y+orTM1y-,i=3或5,

(3d)

TA=min(TMiy+,TMiy-),TB=TM10x+orTM10x-,i=6或8,

(3e)

TA=min(TMix+,TMix-),TB=TM10y+orTM10y-,i=7或9.

(3f)

公式(3a)表述了M點(diǎn)為與G點(diǎn)直接相鄰的10個網(wǎng)格節(jié)點(diǎn)中走時值最小的網(wǎng)格節(jié)點(diǎn)的含義；而公式(3b)與公式(2b)類似，它表述了當(dāng)M=M1or M10時，A點(diǎn)和B點(diǎn)分別為與M點(diǎn)直接相鄰的另外兩個坐標(biāo)軸方向上走時值較小者；但是，公式(3c)—(3f)則不同，此時A點(diǎn)和B點(diǎn)分別為與M點(diǎn)直接相鄰的另外兩個坐標(biāo)軸方向上的鄰近點(diǎn)，但是其中一個方向?yàn)閮蓚€鄰近點(diǎn)的走時值較小者，而另外一個方向上則僅有一個選擇(例如：如圖2c所示，若M=M6時，TA和TB則分別可以表示為：TA=min(TM6y+,TM6y-),TB=TM10x+，其中：TM6y+,TM6y-分別為y軸方向上M6鄰近的點(diǎn)M6y+和M6y-的走時值，而在y軸方向上TB只能取TM10x+).

至此就建立了規(guī)則立方體網(wǎng)格(圖2f)中的局部走時計(jì)算公式.該公式看似和常規(guī)雙線性插值走時計(jì)算的公式區(qū)別不大，但是該公式在其已知條件的確定時，則隱含著基于Fermat原理的迎風(fēng)思想；同時，該思想還能靈活的將網(wǎng)格間距變化的復(fù)雜網(wǎng)格處的局部走時計(jì)算策略簡化到規(guī)則立方體網(wǎng)格中，并保證局部算法的無條件穩(wěn)定性.

(2)鄰近地表的不規(guī)則網(wǎng)格(圖2g)中的局部走時計(jì)算公式

如圖2d所示，當(dāng)計(jì)算地表鄰近點(diǎn)H點(diǎn)的走時值時，其也是處于一種復(fù)雜的不規(guī)則網(wǎng)格中.此時同樣首先采用基于Fermat原理的迎風(fēng)思想，挑選封閉H點(diǎn)的所有網(wǎng)格單元中走時分布最小的網(wǎng)格單元.但是，這種情況比圖2b、c兩種情況要復(fù)雜一些，因?yàn)檫@里除了考慮與H點(diǎn)直接相鄰的網(wǎng)格點(diǎn)外，還需要考慮鄰近地表點(diǎn)的影響，所以此處先給出確定M點(diǎn)走時值TM的表達(dá)式，再根據(jù)TM的取值不同，給出該復(fù)雜局部網(wǎng)格中的走時計(jì)算公式.在該復(fù)雜局部網(wǎng)格中，M點(diǎn)的走時值TM應(yīng)該滿足：

TM=min[(TMi)i=1,2,…,6, (TSi)i=1,2,…,8]

(4)

公式(4)實(shí)際上表述了M點(diǎn)應(yīng)為所有與其相鄰的網(wǎng)格點(diǎn)或地表點(diǎn)中走時值最小的節(jié)點(diǎn).如圖2d所示，當(dāng)M點(diǎn)為不同情況時，TH的計(jì)算公式將不一樣：①當(dāng)M=M1時或M=(Mi)i=2,3,4,5且A=M1x+,M1y-,M1x-orM1y+，則圖2d的復(fù)雜網(wǎng)格可以簡化成圖2f的規(guī)則立方體網(wǎng)格中的局部走時計(jì)算公式，此時根據(jù)雙線性插值方法，直接采用公式(1)即可；②當(dāng)M=(Mi)i=2,3,4 or 5且A=S1,S3,S5orS7時，此時圖2d的復(fù)雜網(wǎng)格可以簡化成圖2g所示的不規(guī)則網(wǎng)格，同樣采用雙線性插值法，可得TH的計(jì)算公式：

(5)

其中：s為當(dāng)前計(jì)算網(wǎng)格單元的平均慢度(速度的倒數(shù))，h為正方體網(wǎng)格的網(wǎng)格間距，d為M點(diǎn)到A點(diǎn)的距離.③當(dāng)M=M6or (Si)i=1,2,…,8時，即M點(diǎn)為地表上的點(diǎn)時，則直接根據(jù)Huygens原理，將M點(diǎn)視為子震源點(diǎn)，則有：

TH=TM+sd，

(6)

其中：d為M點(diǎn)到H點(diǎn)的距離.

至此就建立了鄰近地表的不規(guī)則網(wǎng)格(圖2d)處的局部走時計(jì)算公式，建立時首先需要基于Fermat原理的迎風(fēng)思想給出確定M點(diǎn)的公式(4)，然后分別根據(jù)公式(4)的挑選結(jié)果，將此時的復(fù)雜網(wǎng)格簡化為圖2f的規(guī)則立方體網(wǎng)格、圖2g的不規(guī)則網(wǎng)格或地表上點(diǎn)的直達(dá)波這3種情況中的一種，進(jìn)而對應(yīng)的采用公式(1)、(5)或(6)作為計(jì)算地表附近鄰近H點(diǎn)的局部走時計(jì)算公式，該公式的建立也是綜合利用迎風(fēng)思想、Fermat原理和Huygens原理的結(jié)果，其能靈活地適應(yīng)地表附近的復(fù)雜網(wǎng)格并能保證局部算法的無條件穩(wěn)定性.

(3)地表處的不規(guī)則網(wǎng)格(圖2h)中的局部走時計(jì)算公式

如圖2e所示，當(dāng)計(jì)算地表上的S點(diǎn)的走時值時，其也是處于一種復(fù)雜的不規(guī)則網(wǎng)格中.此時同樣首先采用基于Fermat原理的迎風(fēng)思想，挑選封閉S點(diǎn)的所有網(wǎng)格單元中走時分布最小的網(wǎng)格單元.不過，此時情況變得略微簡單點(diǎn)，因?yàn)槿鐖D2e所示，與S點(diǎn)直接相鄰的網(wǎng)格點(diǎn)僅有M1和(Si)i=1,2,…,8點(diǎn).因此首先根據(jù)基于Fermat原理的迎風(fēng)思想，M點(diǎn)的走時值TM應(yīng)該滿足：

TM=min[TM1,(TSi)i=1,2,…,8].

(7)

同樣當(dāng)M點(diǎn)為不同情況時，TS的計(jì)算公式將不一樣：①當(dāng)M=(Si)i=1,2,…,8時，則可以直接根據(jù)Huygens原理將M點(diǎn)視為子震源點(diǎn)，采用公式(6)的形式計(jì)算TS，不同之處僅為此時的d為M點(diǎn)到S點(diǎn)的距離；②當(dāng)M=M1時，則圖2e的復(fù)雜網(wǎng)格可以簡化成圖2h的不規(guī)則網(wǎng)格，同樣采用雙線性插值法，可得TS的計(jì)算公式：

(8)

其中：s為當(dāng)前計(jì)算網(wǎng)格單元的平均慢度(速度的倒數(shù))，h為正方體網(wǎng)格的網(wǎng)格間距，d為M點(diǎn)到S點(diǎn)的距離.

至此就建立了地表處不規(guī)則網(wǎng)格(圖2e)中的局部走時計(jì)算公式，建立時首先需要基于Fermat原理的迎風(fēng)思想給出確定M點(diǎn)的公式(7)，然后分別根據(jù)公式(7)的挑選結(jié)果，將此時的復(fù)雜網(wǎng)格簡化為圖2h的不規(guī)則網(wǎng)格或地表上點(diǎn)的直達(dá)波這2種情況中的一種，進(jìn)而對應(yīng)的采用公式(6)或(8)作為計(jì)算地表附近鄰近S點(diǎn)的局部走時計(jì)算公式，該公式的建立也是綜合利用迎風(fēng)思想、Fermat原理和Huygens原理的結(jié)果，其能靈活地適應(yīng)地表處的復(fù)雜網(wǎng)格并能保證局部算法的無條件穩(wěn)定性.

綜上所述，為了計(jì)算復(fù)雜黃土塬區(qū)的復(fù)雜網(wǎng)格中的初至波走時，提出了一種迎風(fēng)變網(wǎng)格雙線性插值法，該方法首先通過基于Fermat原理的迎風(fēng)思想，將復(fù)雜網(wǎng)格簡化為圖2f—h所示的3種網(wǎng)格，并綜合利用雙線性插值法和Huygens原理建立了這些網(wǎng)格中的局部走時計(jì)算公式.該方法具有如下優(yōu)點(diǎn)：①采用復(fù)雜網(wǎng)格剖分復(fù)雜黃土塬模型，對地表高程沒做任何近似處理，并在地表附近采用加密的網(wǎng)格用于精細(xì)刻畫地表形態(tài)，進(jìn)而提高對不規(guī)則計(jì)算邊界刻畫的精細(xì)程度；②上述局部走時計(jì)算公式的建立方法，能夠靈活的將復(fù)雜網(wǎng)格中的走時計(jì)算問題，在迎風(fēng)思想、Fermat原理和Huygens原理的綜合應(yīng)用下簡化為3種網(wǎng)格中的雙線性插值問題；③算法能保證在復(fù)雜網(wǎng)格中的任意網(wǎng)格位置均采用計(jì)算精度相對較高的雙線性插值法，進(jìn)而保證整個計(jì)算的精度；④迎風(fēng)思想、Fermat原理和Huygens原理的綜合應(yīng)用能夠保證算法的靈活性和無條件穩(wěn)定性；⑤上述變加密網(wǎng)格僅在地表附近和震源附近進(jìn)行網(wǎng)格加密，因此能夠在不增加過多計(jì)算量的前提下很好地解決震源附近和不規(guī)則計(jì)算邊界處誤差較大的問題.

2.2 復(fù)雜網(wǎng)格中的整體波前擴(kuò)展方式

以上給出了局部走時算法，它可以解決在復(fù)雜網(wǎng)格的不同情況下的局部走時計(jì)算問題.但是，要想完成整個計(jì)算區(qū)域的初至走時計(jì)算，還需要一種整體的波前擴(kuò)展方式.為了很好的解決該問題，在此引入群推進(jìn)法(Kim and Folie, 2000).該方法在三維情況下有很高的計(jì)算效率，同時只需稍加改進(jìn)就能靈活的適應(yīng)復(fù)雜地表的問題.如圖3所示，常規(guī)群推進(jìn)法將計(jì)算區(qū)域內(nèi)的網(wǎng)格節(jié)點(diǎn)的屬性(用Gid的數(shù)值表示)分為三類，即Gid=0、1、2三種(圖3中分別用白色、灰色、黑色充填的原點(diǎn)表示).其中，Gid=0表示未開始走時計(jì)算的網(wǎng)格節(jié)點(diǎn)，Gid=1表示當(dāng)前波前上的網(wǎng)格節(jié)點(diǎn)，Gid=2表示已經(jīng)完成走時計(jì)算的網(wǎng)格節(jié)點(diǎn).

圖3 三維復(fù)雜黃土塬條件下修正后的群推進(jìn)法Fig.3 Modified group marching method in 3D complex loess plateau condition

以上群推進(jìn)法的實(shí)現(xiàn)過程，實(shí)際上是一個通過不斷改變網(wǎng)格節(jié)點(diǎn)屬性，來近似模擬波前傳播的過程，不斷納入波前點(diǎn)的運(yùn)算表征著波前不斷向前傳播到達(dá)新的網(wǎng)格節(jié)點(diǎn)的過程，不斷移除波前點(diǎn)的過程表征著波前到達(dá)過的區(qū)域?yàn)樽邥r計(jì)算完成的區(qū)域.上述波前擴(kuò)展過程是同時選取波前上的很多網(wǎng)格節(jié)點(diǎn)作為子震源點(diǎn)，這與Huygens原理是完全吻合的.此外，為了適應(yīng)復(fù)雜地表問題，只需永久性的將地表以上的網(wǎng)格節(jié)點(diǎn)的屬性設(shè)為Gid=2，走時值設(shè)為無窮大，即可實(shí)現(xiàn)地震波無法穿越地空界面進(jìn)入空氣中物理事實(shí).

綜上所述，綜合采用上述復(fù)雜網(wǎng)格中的局部走時算法和整體波前擴(kuò)展方式，即可無條件穩(wěn)定、在不過多增加計(jì)算量的前提下高精度的計(jì)算復(fù)雜黃土塬區(qū)的初至波走時.當(dāng)然，該算法能否被用于分析復(fù)雜黃土塬區(qū)初至波特征，還需要對算法進(jìn)行精度和效率評估后才能決定，因此接下來首先對算法的精度與效率進(jìn)行分析.

3 計(jì)算精度與效率分析

為了對比分析本文提出的新算法的計(jì)算精度與效率，采用如圖4a所示的黃土塬模型.模型的大小為1.6 km×0.8 km×1.2 km，模型的地表部分由兩個半徑為0.4 km的半球形黃土塬(便于獲得精確的地表上任意一點(diǎn)的解析走時值)構(gòu)成，地下為均勻介質(zhì)，地震波在介質(zhì)中的傳播速度為1.0 km·s-1.計(jì)算時，震源點(diǎn)置于(0.8 km，0.4 km，0.4 km)處.圖4b給出了網(wǎng)格加密策略：在離震源最近的半徑為80.0 m的球體內(nèi)網(wǎng)格間距為2.5 m，依次向外進(jìn)行球形擴(kuò)展，在徑向距離分別為80.0～160.0 m和160.0～320.0 m的球環(huán)內(nèi)，計(jì)算的網(wǎng)格間距分別為5.0 m和10.0 m，然后在320.0～1600.0 m的非加密區(qū)域，網(wǎng)格間距過渡為均勻的20.0 m.上述震源附近的加密策略，相當(dāng)于1.2節(jié)中闡述的(2n,2n-1,2n-2,…,2)指數(shù)衰減的加密策略取加密因子n=3的情況.此外，地表附近則直接采用1.2節(jié)闡述的加密策略.

圖4 計(jì)算精度與效率(a) 模型； (b) 網(wǎng)格加密策略.Fig.4 The computational accuracy and efficiency(a) Model; (b) The strategy for increasing the density of grid.

圖5a—c和圖5d分別給出了當(dāng)單獨(dú)采用網(wǎng)格間距20.0 m、10.0 m、5.0 m和采用上述變加密網(wǎng)格剖分整個模型時，地表上所有網(wǎng)格節(jié)點(diǎn)走時計(jì)算相對誤差的分布，對比分析可知：隨著網(wǎng)格間距的減小，走時計(jì)算的相對誤差逐漸減小(圖5a—c)，但是即使當(dāng)網(wǎng)格間距減少到5.0 m時，整個計(jì)算區(qū)域計(jì)算誤差均存在震源點(diǎn)附近誤差較大的問題；不過，通過采用本文的網(wǎng)格加密策略，卻能很好的解決震源附近誤差較大的問題(圖5d).

圖5 地表上的計(jì)算精度(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h為變加密.Fig.5 The computational accuracy of the earth′s surface(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h is variable density.

圖6和圖7分別給出了Y=0.4 km和X=0.8 km剖面上，各個網(wǎng)格點(diǎn)的走時計(jì)算誤差.與地表上分析的結(jié)果類似，在地表以下的剖面上，同樣存在：隨著網(wǎng)格間距的減小，走時計(jì)算的相對誤差逐漸減小(圖6a—c和圖7a—c)；而本文采用的變加密網(wǎng)格方法(圖6d和圖7d)，可以很好的解決震源附近誤差較大的問題.此外，同時分析圖6和圖7的計(jì)算誤差分布，可以發(fā)現(xiàn)一個有趣的現(xiàn)象：走時計(jì)算誤差的分布具有明顯的方向性.其中，在震源向外輻射的0°、45°、90°、135°、180°、225°、270°和315°的八個方向及附近區(qū)域誤差相對較小，而在這八個方向中的相鄰的兩個方向的平分線方向及附近區(qū)域誤差相對較大.產(chǎn)生這種現(xiàn)象的原因，初步分析為計(jì)算局部走時的計(jì)算公式在均勻介質(zhì)特定方向?yàn)榻馕鼋鉄o誤差.

圖6 Y=0.4 km剖面上的計(jì)算精度(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h為變加密.Fig.6 The computational accuracy of Y=0.4 km section(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h is variable density.

圖7 X=0.8 km剖面上的計(jì)算精度(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h為變加密.Fig.7 The computational accuracy of X=0.8 km section(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h is variable density.

表1給出了計(jì)算精度和計(jì)算效率對比的定量化統(tǒng)計(jì)數(shù)據(jù)，對比分析可以發(fā)現(xiàn)：①無論在哪個統(tǒng)計(jì)區(qū)域(地表上、Y=0.4 km剖面或X=0.8 km剖面上)，當(dāng)網(wǎng)格間距減小時平均相對計(jì)算誤差均減??；②采用上述變加密網(wǎng)格策略獲得的計(jì)算精度(平均相對誤差)，介于網(wǎng)格間距單獨(dú)取10.0 m和5.0 m之間；③變網(wǎng)格策略的計(jì)算耗時(計(jì)算的硬件設(shè)備參數(shù)為：Intel(R) Xeon(R) CPU E5-2620 0 @2.00 GHz，32 G運(yùn)行內(nèi)存)僅為h=10.0 m和h=5.0 m時的19.22%和1.53%.綜上所述，本文的變加密網(wǎng)格策略能在僅增加少量計(jì)算量的前提下，大幅提高計(jì)算精度，并能很好的解決震源附近誤差較大的問題.

4 三維復(fù)雜黃土塬區(qū)初至特征分析

黃土塬區(qū)地震地質(zhì)條件非常復(fù)雜，地表高程劇烈變化(溝、坡、梁、塬等復(fù)雜地貌并存)、表層黃土疏松且風(fēng)化嚴(yán)重、黃土層厚度較大且變化劇烈、表層結(jié)構(gòu)復(fù)雜多變、巨厚黃土層內(nèi)存在低速層、降速層和高速層、黃土層內(nèi)隨機(jī)非均質(zhì)性嚴(yán)重、部分區(qū)域存在明顯的礫石層散射.這些復(fù)雜地震地質(zhì)條件給黃土塬區(qū)地震勘探帶來了巨大的挑戰(zhàn).為獲得三維復(fù)雜黃土塬區(qū)地震初至波走時特征，利用本文所提出的方法計(jì)算了三維水平地表層狀模型和三維復(fù)雜黃土塬模型的地震波走時場.

圖8a給出了三維水平地表層狀介質(zhì)模型的形態(tài)和震源點(diǎn)的位置，其中：模型的大小為1.2 km×1.2 km×2.55 km，由淺及深三層介質(zhì)中地表、界面1和界面2的深度分別為0.15 km、0.37 km和1.4 km；三層介質(zhì)中地震波的傳播速度依次為0.8 km·s-1、2.0 km·s-1和3.0 km·s-1，計(jì)算時震源位于(0.0 km，0.0 km，0.15 km)處，震源和地表附近的網(wǎng)格加密策略采用第3節(jié)“計(jì)算精度與效率分析”中的網(wǎng)格加密方式.圖8b—e分別給出了地表上、Y=0.0 km剖面上、X=0.0 km剖面上和X-Y=0.0 km剖面上的走時分布.

分析圖8b—e可以得出在三維層狀介質(zhì)模型中初至波的一些特征：①在地表上接收的初至波可以分為直達(dá)波和折射波兩個區(qū)域(圖8b中半徑約為0.6 km的1/4紅色虛線圓弧分開的兩個區(qū)域)，其中直達(dá)波等時線間隔距離遠(yuǎn)小于折射波等時線間隔距離(圖8b)；②由于模型為層狀介質(zhì)，其在橫向上是均勻分布的，因此初至波等時線在地表上是以同心圓形式存在的(圖8b)；③在地下介質(zhì)中，由于由淺至深地震波速度是逐層增加的，因此在地表以下兩個界面處都會產(chǎn)生折射波，不過由于各層間的速度比不一樣，其產(chǎn)生折射波的臨界角也有所差別，第1個界面處臨界角較小，其先產(chǎn)生折射波.

圖9a給出了三維復(fù)雜黃土塬模型的形態(tài)和震源點(diǎn)的位置，其中：模型的大小為1.2 km×1.2 km×2.55 km，模型的地表起伏形態(tài)是大致根據(jù)一個實(shí)際黃土塬工區(qū)近地表調(diào)查的黃土層厚度生成的(圖9b)；近地表速度結(jié)構(gòu)復(fù)雜，由淺及深，分別由風(fēng)化層、礫石散射層、低速層、降速層和高速層構(gòu)成，其中：風(fēng)化層、低速層和降速層中地震波的傳播速度分別為：0.5 km·s-1、0.9 km·s-1和1.5 km·s-1，礫石散射層由一個隨機(jī)介質(zhì)構(gòu)成，高速層參考了一個塔西南黃土塬區(qū)深層的復(fù)雜構(gòu)造模型.整體上分析，該三維復(fù)雜黃土塬模型在大尺度上考慮了復(fù)雜構(gòu)造的非均勻性，在幾何形態(tài)上考慮了各種復(fù)雜地形，在小尺度上考慮了礫石層的隨機(jī)擾動，因此其全面細(xì)致的描述了三維復(fù)雜黃土塬區(qū)的復(fù)雜地震地質(zhì)條件.計(jì)算時，震源位于(0.0 km，0.0 km，0.35 km)處.圖9c—f分別給出了地表上、Y=0.0 km剖面上、X=0.0 km剖面上和X-Y=0.0 km剖面上的走時分布.

圖9 三維復(fù)雜黃土塬模型中的初至走時分布(a) 模型； (b) 地表高程的等高線； (c) 地表上走時的等值線在平面上投影； (d) Y=0.0 km剖面上的走時分布； (e) X=0.0 km剖面上的走時分布； (f) X-Y=0.0 km剖面上的走時分布.Fig.9 Travel-time distribution of first arrival in 3D complex loess plateau model(a) Model； (b) The contour of the earth surface′s elevation； (c—f) Travel-time distribution on the earth′s surface，on the section of Y=0.0 km，on the section of X=0.0 km，on the section of X-Y=0.0 km.

分析圖9c—f可以得出復(fù)雜黃土塬模型中初至波的一些特征：①同樣，在地表上接收的初至波可以分為直達(dá)波和折射波兩個區(qū)域(圖9c中半徑約為0.06 km的1/4紅色虛線圓弧分開的兩個區(qū)域)，其中直達(dá)波等時線間隔距離遠(yuǎn)小于折射波等時線間隔距離(圖9c)；但是，由于極淺風(fēng)化低速層的存在，地表上直達(dá)波覆蓋的區(qū)域很小，很快地表上就接收到了折射波；②由于復(fù)雜黃土塬區(qū)地形的劇烈起伏，地表高程變化引起了較大的初至走時時差，因此地表上的初至走時等時線出現(xiàn)了很多極值點(diǎn)(圖9c中的一系列近似同心圓的圓心)；③由于近地表極淺層低速風(fēng)化層的存在，該薄層中折射波發(fā)育，等時線的間距很小(圖9d—f地表附近的等時線)；④礫石層中介質(zhì)速度的隨機(jī)擾動，對初至波的等時線有明顯的隨機(jī)改造，等時線有明顯的局部隨機(jī)擾動(圖9d—f中隨機(jī)擾動層的等時線)；⑤同樣，在地下介質(zhì)中，由于至淺而深地震波速度是逐層增加的，因此在地表以下速度劇變的界面處都會產(chǎn)生折射波(圖9d—f)；⑥由于大體上地下介質(zhì)速度是隨深度的增加而增加的，因此初至波等時線的間隔，由地表向下是逐漸增加的(圖9d—f).

對比分析圖8和圖9可以發(fā)現(xiàn)：①與三維水平地表層狀介質(zhì)模型相比，三維復(fù)雜黃土塬模型地表接收到的初至走時的直達(dá)波等時線形態(tài)在近偏移距同樣為同心圓，但是隨著偏移距增加其形態(tài)受近地表黃土層厚度橫向變化快和地下速度結(jié)構(gòu)復(fù)雜的影響較大，其等時線分布非常復(fù)雜，到達(dá)中遠(yuǎn)偏移距折射波的等時線形態(tài)與地形的等高線有一個很好的對應(yīng)關(guān)系：地形高程的高值閉合區(qū)(圖9b)對應(yīng)著等時線的高值閉合區(qū)(圖9c)，這與地震波傳到黃土塬的頂部需要花更多走時的物理事實(shí)是相吻合的；②與三維水平地表層狀介質(zhì)模型相比，由于近地表薄低速風(fēng)化層的存在，復(fù)雜黃土塬模型中地表初至成分里直達(dá)波覆蓋的區(qū)域比水平層狀模型的要小很多，前者的覆蓋范圍僅為半徑約為0.06 km的1/4紅色虛線圓弧圈定的扇形區(qū)(圖9c)，而后者的覆蓋半徑約為0.6 km(圖8b)；③即使是在如圖9a所示的三維復(fù)雜黃土塬模型中，本文算法也能保證很好的穩(wěn)定性，能很好的適應(yīng)地形劇烈起伏、風(fēng)化層、礫石散射層、低速層、降速層和高速層等近地表及地下復(fù)雜介質(zhì)(圖9c—f).

如圖10所示，為了進(jìn)一步對比分析三維復(fù)雜黃土塬模型中的初至特征，分別提取了三維水平地表層狀介質(zhì)模型各個剖面地表上的時距曲線，以及三維復(fù)雜黃土塬模型各個剖面地表上的時距曲線，對比二者可以發(fā)現(xiàn)：①三維水平地表模型的時距曲線是典型的折線形態(tài)，每段折線的斜率與各層介質(zhì)的速度相關(guān)，其中直達(dá)波和折射波明顯分開，但是三維復(fù)雜黃土塬模型的時距曲線則是不規(guī)則的曲線，其中也能明顯看出直達(dá)波的直線段，但是其形態(tài)也一定程度受到了地表形態(tài)和近地表復(fù)雜介質(zhì)的影響；②三維復(fù)雜黃土塬模型地表上的時距曲線形態(tài)雖然為不規(guī)則的曲線，但是它們的形態(tài)和地表高程線的形態(tài)之間有一定的對應(yīng)關(guān)系，而三維水平地表層狀介質(zhì)模型則沒有；③同樣，圖10中紅色的虛線非常明顯的劃分開了直達(dá)波和折射波的時距曲線的區(qū)間，三維復(fù)雜黃土塬模型的直達(dá)波覆蓋的范圍遠(yuǎn)遠(yuǎn)小于三維水平地表層狀介質(zhì)模型，這主要是三維復(fù)雜黃土塬模型近地表薄風(fēng)化層帶來的影響；④二者形態(tài)的對比可以預(yù)料到，與常規(guī)水平地表地下分層很強(qiáng)的區(qū)域相比，三維復(fù)雜黃土塬區(qū)的初至拾取會因?yàn)闀r距曲線的嚴(yán)重不規(guī)則而面臨著很大的挑戰(zhàn)，同時利用該初至?xí)r距曲線進(jìn)行初至層析獲取近地表速度模型時也將面臨著很大的挑戰(zhàn).

圖10 地表上時距曲線的對比(a)、(c)、(e)三維水平地表層狀介質(zhì)模型的Y=0.0 km、X=0.0 km、X-Y=0.0 km剖面的時距曲線；(b)、(d)、(f)三維復(fù)雜黃土塬模型的Y=0.0 km、X=0.0 km、X-Y=0.0 km剖面的時距曲線.Fig.10 The comparison of time-distance curves on the earth′s surface(a),(c),(e) The time-distance curves of the 3D model with horizontal earth′s surface and layered medium on the section of Y=0.0 km，X=0.0 km and X-Y=0.0 km；(b),(d),(f) The time-distance curves of the 3D complex loess plateau model on the section of Y=0.0 km，X=0.0 km and X-Y=0.0 km.

5 結(jié)論與討論

為了快速穩(wěn)定的計(jì)算三維復(fù)雜黃土塬區(qū)初至波走時并分析其特征，提出一種群推進(jìn)迎風(fēng)變網(wǎng)格雙線性插值法；通過對該方法的理論分析、計(jì)算精度和效率以及計(jì)算實(shí)例的分析，可以得出如下結(jié)論：

(1)采用三維分區(qū)局部變加密不等距網(wǎng)格剖分三維復(fù)雜黃土塬模型，具有精細(xì)刻畫地表形態(tài)、提高地表及震源附近計(jì)算精度、在增加很少計(jì)算量的前提下大幅提高整個模型區(qū)域的走時計(jì)算精度、保證整個計(jì)算區(qū)域穩(wěn)定性等優(yōu)勢；

(2)計(jì)算精度和計(jì)算效率對比分析表明：為了計(jì)算三維復(fù)雜黃土塬區(qū)復(fù)雜網(wǎng)格中的初至波走時而提出的一種迎風(fēng)變網(wǎng)格雙線性插值法，具有能無條件穩(wěn)定適應(yīng)復(fù)雜地形和近地表及地下復(fù)雜介質(zhì)、靈活穩(wěn)定適應(yīng)復(fù)雜網(wǎng)格、保證計(jì)算精度、在不過多增加計(jì)算量的前提下解決震源附近誤差較大的問題等優(yōu)勢；

(3)計(jì)算實(shí)例分析表明：本文算法能穩(wěn)定適應(yīng)三維復(fù)雜黃土塬地震地質(zhì)條件，同時還得出了三維復(fù)雜黃土塬區(qū)初至成分組成及分布特征，對于充分拾取和利用復(fù)雜黃土塬區(qū)的初至走時信息具有一定的理論意義和實(shí)際價值.