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

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

0 引言

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

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

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

1 三維復雜黃土塬模型的網格剖分

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

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

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

圖1 三維復雜黃土塬模型的網格剖分Fig.1 Mesh generation of 3D complex loess plateau model

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

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

2.1 復雜網格中的局部走時算法

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

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

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

(1)

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

圖2 迎風變網格雙線性插值法(a) 三維復雜地表附近的局部網格； (b—e) 不同位置的網格形態(tài)；(f—h)簡化后的計算網格.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點為與G點直接相鄰的6個網格節(jié)點中走時值最小的網格節(jié)點的含義，而公式(2b)—(2d)則表述了A點和B點分別為與M點直接相鄰的另外兩個坐標軸方向上走時值較小者(例如：如圖2b所示，若M=M6時，TA和TB則分別可以表示為：TA=min(TM6x+,TM6x-),TB=min(TM6y+,TM6y-)，其中：TM6x+,TM6x-分別為x軸方向上M6鄰近的點M6x+和M6x-的走時值，TM6y+,TM6y-則為y軸方向上M6鄰近的點M6y+和M6y-的走時值)；②如圖2c所示，當N點位于網格間距有變化的不均勻立方體網格時，同樣根據Fermat原理來確定網格單元MACB及其中的端點M點、A點、C點、B點的走時值TM、TA、TC、TB，但是此時情況變得復雜了很多，其中M點應為與N點直接相鄰的10個網格節(jié)點中走時值最小的網格節(jié)點，而A點和B點的確定方法已有所變化，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點為與G點直接相鄰的10個網格節(jié)點中走時值最小的網格節(jié)點的含義；而公式(3b)與公式(2b)類似，它表述了當M=M1or M10時，A點和B點分別為與M點直接相鄰的另外兩個坐標軸方向上走時值較小者；但是，公式(3c)—(3f)則不同，此時A點和B點分別為與M點直接相鄰的另外兩個坐標軸方向上的鄰近點，但是其中一個方向為兩個鄰近點的走時值較小者，而另外一個方向上則僅有一個選擇(例如：如圖2c所示，若M=M6時，TA和TB則分別可以表示為：TA=min(TM6y+,TM6y-),TB=TM10x+，其中：TM6y+,TM6y-分別為y軸方向上M6鄰近的點M6y+和M6y-的走時值，而在y軸方向上TB只能取TM10x+).

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

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

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

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

(4)

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

(5)

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

TH=TM+sd，

(6)

其中：d為M點到H點的距離.

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

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

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

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

(7)

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

(8)

其中：s為當前計算網格單元的平均慢度(速度的倒數)，h為正方體網格的網格間距，d為M點到S點的距離.

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

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

2.2 復雜網格中的整體波前擴展方式

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

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

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

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

3 計算精度與效率分析

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

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

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

圖5 地表上的計算精度(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剖面上，各個網格點的走時計算誤差.與地表上分析的結果類似，在地表以下的剖面上，同樣存在：隨著網格間距的減小，走時計算的相對誤差逐漸減小(圖6a—c和圖7a—c)；而本文采用的變加密網格方法(圖6d和圖7d)，可以很好的解決震源附近誤差較大的問題.此外，同時分析圖6和圖7的計算誤差分布，可以發(fā)現一個有趣的現象：走時計算誤差的分布具有明顯的方向性.其中，在震源向外輻射的0°、45°、90°、135°、180°、225°、270°和315°的八個方向及附近區(qū)域誤差相對較小，而在這八個方向中的相鄰的兩個方向的平分線方向及附近區(qū)域誤差相對較大.產生這種現象的原因，初步分析為計算局部走時的計算公式在均勻介質特定方向為解析解無誤差.

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

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

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

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

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

圖9a給出了三維復雜黃土塬模型的形態(tài)和震源點的位置，其中：模型的大小為1.2 km×1.2 km×2.55 km，模型的地表起伏形態(tài)是大致根據一個實際黃土塬工區(qū)近地表調查的黃土層厚度生成的(圖9b)；近地表速度結構復雜，由淺及深，分別由風化層、礫石散射層、低速層、降速層和高速層構成，其中：風化層、低速層和降速層中地震波的傳播速度分別為：0.5 km·s-1、0.9 km·s-1和1.5 km·s-1，礫石散射層由一個隨機介質構成，高速層參考了一個塔西南黃土塬區(qū)深層的復雜構造模型.整體上分析，該三維復雜黃土塬模型在大尺度上考慮了復雜構造的非均勻性，在幾何形態(tài)上考慮了各種復雜地形，在小尺度上考慮了礫石層的隨機擾動，因此其全面細致的描述了三維復雜黃土塬區(qū)的復雜地震地質條件.計算時，震源位于(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 三維復雜黃土塬模型中的初至走時分布(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可以得出復雜黃土塬模型中初至波的一些特征：①同樣，在地表上接收的初至波可以分為直達波和折射波兩個區(qū)域(圖9c中半徑約為0.06 km的1/4紅色虛線圓弧分開的兩個區(qū)域)，其中直達波等時線間隔距離遠小于折射波等時線間隔距離(圖9c)；但是，由于極淺風化低速層的存在，地表上直達波覆蓋的區(qū)域很小，很快地表上就接收到了折射波；②由于復雜黃土塬區(qū)地形的劇烈起伏，地表高程變化引起了較大的初至走時時差，因此地表上的初至走時等時線出現了很多極值點(圖9c中的一系列近似同心圓的圓心)；③由于近地表極淺層低速風化層的存在，該薄層中折射波發(fā)育，等時線的間距很小(圖9d—f地表附近的等時線)；④礫石層中介質速度的隨機擾動，對初至波的等時線有明顯的隨機改造，等時線有明顯的局部隨機擾動(圖9d—f中隨機擾動層的等時線)；⑤同樣，在地下介質中，由于至淺而深地震波速度是逐層增加的，因此在地表以下速度劇變的界面處都會產生折射波(圖9d—f)；⑥由于大體上地下介質速度是隨深度的增加而增加的，因此初至波等時線的間隔，由地表向下是逐漸增加的(圖9d—f).

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

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

圖10 地表上時距曲線的對比(a)、(c)、(e)三維水平地表層狀介質模型的Y=0.0 km、X=0.0 km、X-Y=0.0 km剖面的時距曲線；(b)、(d)、(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 結論與討論

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

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

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

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