劉四新，朱怡諾，王旭東，宋二喬，賀文博

吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

地震波在向介質(zhì)中傳播的過(guò)程中，當(dāng)遇到地層界面下伏地層速度高于上覆地層速度，且入射角大于或等于臨界角時(shí)會(huì)發(fā)生地震波的折射，地震波沿折射界面滑行，并在上覆介質(zhì)中產(chǎn)生折射波。區(qū)別于研究地殼結(jié)構(gòu)的寬角折射勘探(勘探深度一般為幾千米至幾十千米不等，縱波主頻一般為幾赫茲至幾十赫茲不等)，工程地震折射波勘探的深度一般在1 km以?xún)?nèi)，縱波主頻可達(dá)150 Hz以上。工程地震折射波勘探記錄在地表激發(fā)并在近地表地質(zhì)結(jié)構(gòu)中傳播后返回地面的人工地震波，從記錄中獲得波到達(dá)地面的時(shí)間，即波的旅行時(shí)，通過(guò)選取不同的折射波解釋方法來(lái)推算地下不同巖層分界面的埋藏深度、傾角、地層速度等要素，從而了解近地表地層的物理參數(shù)和地質(zhì)形態(tài)。

近年來(lái)，針對(duì)近地表地質(zhì)勘查的工程地球物理方法受到越來(lái)越廣泛的關(guān)注，包括地震勘探、電磁勘探、電法勘探、磁法勘探以及地質(zhì)雷達(dá)等。在淺層地震勘探中，常使用的方法有折射波勘探法、反射波勘探法以及面波勘探法等。折射波勘探法是一種使用較久且成熟的方法，它使用大震源，能快速覆蓋長(zhǎng)距離用以探測(cè)近地表目標(biāo)，且能夠刻畫(huà)存在明顯速度對(duì)比的隆起構(gòu)造，具有操作簡(jiǎn)單、工作效率高、初至波易于識(shí)取、解釋方便等優(yōu)點(diǎn)，這些都是反射波法無(wú)法比擬的。

淺層地震折射波法早在20世紀(jì)30年代就已出現(xiàn)并廣泛應(yīng)用于民用工程勘探。折射波法由石油地震勘探引用而來(lái)，因而在最開(kāi)始的勘測(cè)中，震源和檢波器都沿用原有的設(shè)備。到1940年時(shí)，形成了針對(duì)地震折射波的標(biāo)準(zhǔn)縱測(cè)線(xiàn)觀(guān)測(cè)系統(tǒng)。20世紀(jì)50年代，折射波勘探的各種技術(shù)方法和專(zhuān)用儀器取得了突破性的進(jìn)展。20世紀(jì)70—80年代期間，折射波勘探解釋方法由低精度、低效的人工解釋發(fā)展為高精度、高效的計(jì)算機(jī)自動(dòng)化解釋[1-2]。

目前工程折射波法常用來(lái)探測(cè)覆蓋層厚度[3]，研究基巖面起伏[4]，識(shí)別隱伏斷層，測(cè)量潛水面深度，探測(cè)地下空洞和隧道[5]。Rucker[6]指出地震折射波方法可以為工程地質(zhì)應(yīng)用提供有效的近地表地質(zhì)信息；王洪[7]在實(shí)際勘測(cè)中提出并應(yīng)用了一套工程地震折射波勘探方法，在堤壩病害勘查中取得了很好的效果。

傳統(tǒng)的地震折射解釋方法將近地表介質(zhì)解釋為每層都有一個(gè)單獨(dú)速度的簡(jiǎn)單離散層狀結(jié)構(gòu)，基于這個(gè)假設(shè)，通過(guò)測(cè)量地震折射波的走時(shí)數(shù)據(jù)便可以獲得每層的層厚和地層傾角。近年來(lái)，地震折射波處理解釋方法取得了進(jìn)展，產(chǎn)生了一種新的解釋方法，即折射層析成像技術(shù)(SRT)，該方法通過(guò)對(duì)地震剖面進(jìn)行網(wǎng)格剖分實(shí)現(xiàn)了獲得連續(xù)變化的地震波傳播速度。SRT成為對(duì)地震波傳播的物理特性進(jìn)行三維成像的主要方法之一[8]。近年來(lái)的一些研究證明了SRT在近地表速度模型和地下探測(cè)的高適用性，他們指出在識(shí)別速度縱橫向變化梯度方面，SRT比傳統(tǒng)折射解釋方法取得了更好的效果[9]。

本文介紹了常用折射波解釋的方法原理，著重介紹了折射層析成像方法，并闡述了這些方法各自的優(yōu)缺點(diǎn)及適用性，同時(shí)對(duì)國(guó)內(nèi)外工程地震折射波解釋方法的發(fā)展歷程及現(xiàn)狀進(jìn)行闡述。

1 方法原理

工程地震折射波解釋主要用于研究近地表介質(zhì)的速度、層厚及結(jié)構(gòu)，可分為根據(jù)已知的勘測(cè)結(jié)果估計(jì)地下結(jié)構(gòu)的定性解釋?zhuān)耙罁?jù)定性解釋的結(jié)果通過(guò)具體方法進(jìn)行計(jì)算得到準(zhǔn)確地下結(jié)構(gòu)的定量解釋。

折射波解釋方法也可分為解析法[10-18]和數(shù)值法[19]兩大類(lèi)(圖1)。解析法發(fā)展較早，包括波前法[10]、延遲時(shí)間法(又稱(chēng)時(shí)間項(xiàng)法)[11]、截距時(shí)間法[12]、哈萊斯法[13]、共軛點(diǎn)法[14]、t0差數(shù)法(又稱(chēng)加減法、ABC法)[15]、互換法(即表層剝?nèi)シ?、去表層?[16]、廣義互換法[17]、褶積成像法[18]等；數(shù)值法則較晚，其中應(yīng)用最多的是折射層析成像[19]，一般是先給定初始模型，對(duì)初至波進(jìn)行射線(xiàn)追蹤，采用線(xiàn)性非線(xiàn)性算法進(jìn)行反演確定最接近實(shí)際地層介質(zhì)的屬性結(jié)構(gòu)。折射波解釋的數(shù)值方法取決于地層模型、射線(xiàn)追蹤方法和最優(yōu)化算法。該類(lèi)方法精度更高、效果更好、更針對(duì)復(fù)雜結(jié)構(gòu)。在各類(lèi)地震折射波法中，折射波層析成像是目前的研究熱點(diǎn)和未來(lái)趨勢(shì)。

1.1 波前法

波前法又稱(chēng)時(shí)間場(chǎng)法，是Thornburg[10]根據(jù)惠更斯原理提出的一種基于斯奈爾定理的圖解法。Baumgarte[20]對(duì)波前法進(jìn)行過(guò)論述，認(rèn)為這類(lèi)方法是最精確的折射解釋技術(shù)，但是過(guò)于麻煩費(fèi)時(shí)。經(jīng)Schenck[21]改進(jìn)后，畫(huà)波前面的時(shí)間大大減少，更適用于數(shù)字處理技術(shù)，但實(shí)際應(yīng)用仍較少。

波前法利用波場(chǎng)外推公式分別計(jì)算出地下介質(zhì)中各網(wǎng)格點(diǎn)的地震波到時(shí)重建地震波走時(shí)場(chǎng)，并可獲得各個(gè)時(shí)刻的波前位置。設(shè)t1(x,z)和t2(x,z)分別表示正反向走時(shí)場(chǎng)，T表示相遇觀(guān)測(cè)系統(tǒng)的互換時(shí)間，折射界面上任意點(diǎn)R到時(shí)滿(mǎn)足

t1(xR,zR)+t2(xR,zR)=T。

(1)

式中：xR表示點(diǎn)R橫坐標(biāo)；zR表示點(diǎn)R縱坐標(biāo)。則可據(jù)此找出折射界面的位置。折射層的速度可由重構(gòu)的波前傳播時(shí)間估算，即用界面上兩點(diǎn)之間的距離除以?xún)牲c(diǎn)間重構(gòu)的波前時(shí)間差值。

圖1 折射波解釋方法發(fā)展歷程Fig.1 Progress of refraction interpretation method

原則上這個(gè)方法適用于任意多層情況，通過(guò)遞推的方法順序求出各層參數(shù)。但如果目的層以上的地層較為簡(jiǎn)單，那么波前法和截距時(shí)間法、加減法同樣能解決問(wèn)題。當(dāng)存在隱蔽層和速度劇烈變化的地層時(shí)，圓弧形波前不適用上述情況，因此波前法精度低于理論值。波前法是t0差數(shù)法和哈萊斯法的基礎(chǔ)。

1.2 延遲時(shí)間法

延遲時(shí)間法又稱(chēng)時(shí)間項(xiàng)法[11]。延遲時(shí)間是指觀(guān)測(cè)臨界折射時(shí)間(校正到一個(gè)基準(zhǔn)面上的時(shí)間)與假設(shè)該折射面移至地表或基準(zhǔn)面時(shí)所應(yīng)有的折射時(shí)間之差。主要步驟為：1)校正全部波至?xí)r間并給定折射波速度；2)計(jì)算收發(fā)點(diǎn)的延遲時(shí)間并進(jìn)行時(shí)深轉(zhuǎn)換。該方法操作和校正簡(jiǎn)便，但由于求取速度方法始終不佳及對(duì)傾角敏感等原因，近年來(lái)已經(jīng)不常使用。

1.3 截距時(shí)間法

截距時(shí)間法由Ewing等[12]首先提出，后來(lái)Adachi[22]、Mota[23]、Jonhson[24]又用模型，即深度和射線(xiàn)路徑參數(shù)給予了說(shuō)明。

根據(jù)截距時(shí)間法折射波時(shí)距曲線(xiàn)圖(圖2)可知，直達(dá)波時(shí)距曲線(xiàn)斜率的倒數(shù)為覆蓋層速度v1，折射波時(shí)距曲線(xiàn)斜率的倒數(shù)為折射層速度v2，S1、S2分別為炮點(diǎn)O1、O2產(chǎn)生的初至波時(shí)距曲線(xiàn)。將折射波時(shí)距曲線(xiàn)的截距t01與t02定義為截距時(shí)間，由此可算出t01、t02處折射界面的法線(xiàn)深度h1與h2。

(2)

再分別以O(shè)i為圓心，以求出的hi為半徑做弧，則這兩個(gè)弧的公切線(xiàn)即為折射界面。

圖2 截距時(shí)間法原理圖Fig.2 Principle of intercept time method

該方法僅適用于地層界面為水平或傾斜的平界面、地表無(wú)起伏且速度無(wú)橫向變化(即每層速度為定值)的淺層地震解釋。對(duì)于淺層界面，該方法誤差較小，但對(duì)于深層界面或有薄夾層時(shí)，誤差較大。然而實(shí)際地形常有起伏，折射界面深度與炮點(diǎn)處深度不盡相同，截距時(shí)間法近年來(lái)在實(shí)際中應(yīng)用很少；因而發(fā)展出了加減法、互換法等適用于起伏界面的方法。

1.4 哈萊斯法

哈萊斯法通過(guò)作圖計(jì)算臨界折射時(shí)間，再求取折射界面深度進(jìn)行構(gòu)造解釋[13]。

如圖3所示，O1、O2為炮點(diǎn)，M、N為接收點(diǎn)，P為折射界面上任意點(diǎn)，界面傾角為φ，O1APN和O2BPM為相遇折射波射線(xiàn)路徑，過(guò)P點(diǎn)作界面的垂線(xiàn)，交地面于O點(diǎn)，過(guò)M、N分別做PO的垂線(xiàn)交PO于E、F，以MN為底邊，以臨界折射角ic為底角作等腰三角形MNQ，可以證明Q點(diǎn)在PO的延長(zhǎng)線(xiàn)上。以一系列Q點(diǎn)為圓心，相應(yīng)的PQ為半徑作弧，這些弧的包絡(luò)面即為所求界面。

哈萊斯法適合地表平整、折射界面起伏劇烈、下伏地層速度存在差異的情況，王奇[25]將其應(yīng)用于實(shí)測(cè)中取得較好的效果(圖4)，但必須已知上覆地層速度；當(dāng)存在隱蔽層、變速層時(shí)效果不理想，且不易于作圖。

1.5 共軛點(diǎn)法

針對(duì)哈萊斯法的問(wèn)題，發(fā)展出了共軛點(diǎn)法[14]，即用共軛點(diǎn)的方法求取M、N，再與哈萊斯法結(jié)合進(jìn)行計(jì)算。陳滋康[26]用其對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行解釋?zhuān)?shí)現(xiàn)了計(jì)算機(jī)自動(dòng)處理。該方法無(wú)需已知覆蓋層和折射層速度，只根據(jù)一組相遇時(shí)距曲線(xiàn)即可求解，然而要求在剖面上至少有一個(gè)點(diǎn)的界面深度是已知的。

圖3 哈萊斯法原理圖Fig.3 Principle of Hales method

1.6 t0差數(shù)法

t0差數(shù)法在最初稱(chēng)為加減法，又稱(chēng)ABC法，是由Hagedoorn[15]提出的基于幾何地震學(xué)原理來(lái)求取傾角較小均勻?qū)咏缑娴囊环N方法，其實(shí)質(zhì)是在每一個(gè)檢波器處使用截距時(shí)間法。該方法定義了加、減兩個(gè)時(shí)間值，加時(shí)間線(xiàn)描繪了折射層界面形狀，減時(shí)間線(xiàn)可用于求取折射層速度。

如圖5所示，i為入射角，3對(duì)炮檢距射線(xiàn)路徑分別為O1ABS、O2DCS和O1ADO2，定義t+和t-為：

t+=tO1ABS+tO2DCS-tO1ADO2，

(3)

t-=tO1ABS-tO2DCS+tO1ADO2。

(4)

式中，tO1ABS、tO2DCS和tO1ADO2分別是對(duì)O1ABS、O2DCS和O1ADO2射線(xiàn)路徑測(cè)得的地震波初至?xí)r間。則折射層速度v2為

(5)

S點(diǎn)處界面埋深h為

(6)

依次類(lèi)推可獲得測(cè)線(xiàn)上全部點(diǎn)的界面埋深。

該方法在界面的曲率半徑遠(yuǎn)大于其埋深的情況下效果很好，適用于簡(jiǎn)單不規(guī)則起伏界面和界面兩側(cè)地震波速度變化較大的界面，但不適用于傾角大于10°的情況，且橫向速度不能有明顯變化。t0差數(shù)法比波前法速度快、使用成本低，是目前使用最為廣泛的淺層折射波勘探方法。熊章強(qiáng)等[27]將該方法應(yīng)用于隧道勘查，得到了準(zhǔn)確的基巖埋深和各地層速度信息；彭驍?shù)萚28]在候家梁隧道工區(qū)探測(cè)中取得了較好的效果；包勛等[29]成功利用t0差數(shù)法探測(cè)到隱伏斷層的位置和分布情況(圖6)。

據(jù)文獻(xiàn)[25]。圖4 遼寧鞍山境內(nèi)某隧道哈萊斯法解釋圖Fig.4 Interpretation of a tunnel in Anshan, Liaoning Province

圖5 t0差數(shù)法原理圖Fig.5 Principle of t0difference method

1.7 互換法

互換法又稱(chēng)表層剝?nèi)シ?、去表層法，?yīng)用起來(lái)最為簡(jiǎn)潔并能估算折射層結(jié)構(gòu)和速度情況[16]。

該方法引入了延遲時(shí)的概念。如圖7所示，將A點(diǎn)延遲時(shí)定義為烈時(shí)，可得到更準(zhǔn)確的深度；但當(dāng)存在大構(gòu)造時(shí)，解釋界面因虛假速度變得圓滑。該方法與t0差數(shù)法結(jié)合使用，可求取界面深度小于25 m的地層。王奇[25]將其應(yīng)用于隧道探測(cè)中，與t0差數(shù)法進(jìn)行比較取得了相似的效果(圖8、9)。

(7)

則A點(diǎn)下折射界面深度為

(8)

根據(jù)幾何方法，可得D點(diǎn)延遲時(shí)間：

(9)

可見(jiàn)，延遲時(shí)tD等于t0差數(shù)法t0的一半，也是截距時(shí)間的一半。當(dāng)深度無(wú)較大起伏且速度差異劇

據(jù)文獻(xiàn)[29]。圖6 某水庫(kù)折射波綜合時(shí)距曲線(xiàn)及t0差數(shù)法解釋剖面結(jié)果圖Fig.6 A reservoir refraction synthetical traveltime and profile result for t0 difference method

圖7 互換法幾何路徑圖Fig.7 Geometric path diagram of reciprocal method

1.8 廣義互換法

廣義互換法是互換法的拓展，也需要正反兩個(gè)時(shí)距曲線(xiàn)，其計(jì)算速度和時(shí)間深度的公式與傳統(tǒng)互換法相似，但增加了一個(gè)補(bǔ)償項(xiàng)進(jìn)行偏移校正處理[17]。

如圖10所示，X、Y為接收點(diǎn)，如果地表上的D點(diǎn)與接收點(diǎn)不重合，則D點(diǎn)延遲時(shí)為

(10)

(11)

廣義互換法大致可分為3個(gè)過(guò)程：1)求取速度分析函數(shù)；2)求取時(shí)深偏移函數(shù)；3)選取最佳XY值。該方法通過(guò)選取最佳XY值并帶入步驟1)和2)中的函數(shù)可求出折射層速度和深度。當(dāng)XY=0時(shí)，廣義互換法即為通常意義上的互換法，可能會(huì)導(dǎo)致虛假折射速度大量出現(xiàn)，并使得有起伏折射面的細(xì)節(jié)被圓滑掉。該方法處理傾角小于20°的界面效果一般，而應(yīng)用于大傾角、存在橫向速度變化或隱蔽帶時(shí)效果較好[31]。Sj?gren[32]嘗試將廣義互換法用于淺層工程調(diào)查中并與哈萊斯法進(jìn)行比較，驗(yàn)證了其適應(yīng)性，并總結(jié)了勘查過(guò)程中可能會(huì)遇到的問(wèn)題(圖11)。

據(jù)文獻(xiàn)[25]。圖8 遼寧鞍山境內(nèi)某隧道折射波互換法解釋圖Fig.8 Interpretation of reciprocal method of a tunnel in Anshan, Liaoning Province

據(jù)文獻(xiàn)[25]。圖9 遼寧鞍山境內(nèi)某隧道折射波t0差數(shù)法解釋圖Fig.9 Interpretation of t0difference method of a tunnel in Anshan, Liaoning Province

圖10 廣義互換法確定延遲時(shí)原理圖Fig.10 Principle of generalized reciprocal method to calculate delay time

1.9 褶積成像法

Palmer[18]提出了運(yùn)用褶積的淺層地震折射波處理方法，同年考慮了幾何擴(kuò)散的影響，對(duì)折射波褶積成像法做出了改進(jìn)。該方法不需要拾取初至?xí)r間，可以直接將同一點(diǎn)處的正向和反向走時(shí)數(shù)據(jù)進(jìn)行褶積得到走時(shí)數(shù)據(jù)的交叉時(shí)，再利用速度求出該點(diǎn)的垂向深度。Franco[33]提出了利用褶積和互相關(guān)運(yùn)算對(duì)數(shù)據(jù)疊加來(lái)提高信噪比的方法，并對(duì)多個(gè)折射界面成像，是常規(guī)解釋方法的一大突破(圖12)。

1.10 折射波層析成像

層析成像技術(shù)最先用于醫(yī)學(xué)，隨后拓展到地球物理學(xué)界。Aki等[34]最早提出地震走時(shí)層析成像方法，并應(yīng)用于美國(guó)某地區(qū)的地殼結(jié)構(gòu)探測(cè)。Hearn等[19]利用6年的天然地震數(shù)據(jù)首次實(shí)現(xiàn)了折射波層析成像。

在工程地震層析成像方法中，目前應(yīng)用較多的是基于射線(xiàn)理論的折射波走時(shí)層析成像。與傳統(tǒng)折射波法相比，折射波層析成像的優(yōu)勢(shì)在于對(duì)速度縱橫向有較大變化的地層、大傾角地層、隱伏層、界面起伏層等情況反演結(jié)果很好。

地震折射波層析成像計(jì)算流程主要包括初至走時(shí)提取、速度模型建立、正演、迭代反演4個(gè)步驟。

1.10.1 初至走時(shí)提取

地震勘探中把經(jīng)過(guò)地下介質(zhì)首先到達(dá)檢波器的波稱(chēng)為初至波，這些波的旅行時(shí)包含了淺層速度信息，具有易于獲取、能量強(qiáng)、可追蹤性好的優(yōu)點(diǎn)[35]。初至波走時(shí)的準(zhǔn)確拾取構(gòu)成了折射層析成像的基礎(chǔ)和前提。

a. 哈萊斯法；b. 廣義互換法。據(jù)文獻(xiàn)[32]。圖11 速度和深度解釋結(jié)果圖Fig.11 Results of speed and depth interpretation

據(jù)文獻(xiàn)[33]。圖12 意大利北部克盧索內(nèi)盆地實(shí)測(cè)數(shù)據(jù)Fig.12 Application to real data of Clusone basin area in Northern Italy

1.10.2 速度模型建立

走時(shí)層析成像的效果和收斂速度依賴(lài)于初始速度模型，建立近地表初始速度模型的方法一般有3種：1)人工給定常速度初始速度模型；2)根據(jù)工區(qū)已有的速度資料確定初始速度模型；3)利用直達(dá)波和折射波初至?xí)r間自動(dòng)求取初始速度模型(該模型每層速度都是常數(shù)，在此基礎(chǔ)上可以進(jìn)行層間速度的反距離加權(quán)差值計(jì)算，獲得縱向漸變的初始速度模型)[36]。段心標(biāo)等[37]在此基礎(chǔ)上提出一種生成模型網(wǎng)格節(jié)點(diǎn)的初始速度方法。劉玉柱等[38]通過(guò)研究波形反演目標(biāo)函數(shù)性態(tài)分析，驗(yàn)證了初至波走時(shí)層析成像對(duì)初始模型的依賴(lài)性，討論了簡(jiǎn)單、高效的層析初始模型選取問(wèn)題。

1.10.3 正演

地震層析成像對(duì)正演計(jì)算要求很高，其精度和速度對(duì)成像的分辨率和可靠程度具有很大影響。常用兩種數(shù)值模擬：其一是以射線(xiàn)理論為基礎(chǔ)的射線(xiàn)追蹤法；其二是以波動(dòng)理論為基礎(chǔ)的波動(dòng)方程數(shù)值模擬。由于計(jì)算上的優(yōu)勢(shì)，目前折射波層析成像主要為射線(xiàn)理論層析成像方法。

1)射線(xiàn)追蹤法

基于射線(xiàn)理論的地震波走時(shí)層析成像方法是因其直觀(guān)性、高效率、模型強(qiáng)適應(yīng)性而被應(yīng)用的一種波場(chǎng)近似算法，主要有兩種經(jīng)典的算法：試射法(打靶法)[39]和彎曲法[40]。在此基礎(chǔ)上發(fā)展演化出了兩大類(lèi)現(xiàn)代算法：一是基于全局算法的最短路徑法[41]和線(xiàn)性走時(shí)插值算法[42]，在計(jì)算過(guò)程中同時(shí)考慮所有離散點(diǎn)上的走時(shí)和射線(xiàn)路徑；二是基于局部算法，在射線(xiàn)追蹤過(guò)程中只考慮兩點(diǎn)間的走時(shí)和射線(xiàn)路徑。

試射法是根據(jù)斯奈爾定理，通過(guò)不斷調(diào)整震源點(diǎn)處的射線(xiàn)角來(lái)完成對(duì)接收點(diǎn)的追蹤，進(jìn)而確定射線(xiàn)路徑。這種方法能精細(xì)地處理彎曲界面，適用于一般各向異性介質(zhì)的多值走時(shí)計(jì)算，但對(duì)于復(fù)雜的地質(zhì)模型存在陰影問(wèn)題。

彎曲法是根據(jù)費(fèi)馬原理，通過(guò)調(diào)整假定參數(shù)以達(dá)到震源與檢波點(diǎn)間走時(shí)最小的一種方法。這種方法雖然不存在陰影區(qū)，但用于局部復(fù)雜的介質(zhì)模型會(huì)收斂到局部最小值，且計(jì)算效率和精度較低。

最短路徑法又稱(chēng)最小走時(shí)樹(shù)法，由Nakanishi等[41]提出，由Moser[43]、Cheng[44]改進(jìn)。該方法基于惠更斯原理和網(wǎng)格理論，計(jì)算從震源到模型中所有網(wǎng)格點(diǎn)的走時(shí)和射線(xiàn)路徑，給出震源到接收點(diǎn)所有可能的最小走時(shí)節(jié)點(diǎn)和射線(xiàn)路徑，保證了全局收斂性，簡(jiǎn)單且穩(wěn)定，不會(huì)錯(cuò)過(guò)任何一個(gè)接收點(diǎn)。

線(xiàn)性走時(shí)插值法由Asakawa等[42]提出，由Cardarelli等[45]、張東等[46]改進(jìn)，包括兩個(gè)步驟：正向計(jì)算所有網(wǎng)格節(jié)點(diǎn)的走時(shí)；反向追蹤所有收發(fā)排列的射線(xiàn)路徑。這種算法的精度和計(jì)算效率都很高。

基于局部算法的射線(xiàn)追蹤技術(shù)中最具代表性的是基于程函方程有限差分解的射線(xiàn)追蹤方法[47]，其思想是把射線(xiàn)追蹤問(wèn)題分為正向計(jì)算波前走時(shí)和反向追蹤射線(xiàn)路徑兩個(gè)步驟。該方法計(jì)算速度快，不存在傳統(tǒng)方法具有的陰影問(wèn)題；但當(dāng)局部速度變化過(guò)大時(shí)計(jì)算可能出現(xiàn)負(fù)數(shù)開(kāi)方問(wèn)題，導(dǎo)致算法不穩(wěn)定，且其擴(kuò)展方式不符合波前傳播的物理規(guī)律。基于此，一些學(xué)者研究發(fā)展出了快速掃描算法(fast sweeping method)[48]、快速推進(jìn)算法(fast marching method)[49]、高精度快速推進(jìn)算法[50]、基于多模塊快速推進(jìn)算法(multistencils fast marching method)[51]。

2)波動(dòng)方程數(shù)值模擬

射線(xiàn)路徑無(wú)法被準(zhǔn)確地追蹤到，其原因可能是構(gòu)造較復(fù)雜，或存在間斷面；而波動(dòng)方程法正演具有波場(chǎng)齊全、信息豐富的特點(diǎn)，缺點(diǎn)在于費(fèi)時(shí)。地震波場(chǎng)的正演模擬是利用不同算法求解變系數(shù)偏微分方程的過(guò)程。目前求解方法包括有限差分法[52]、有限元法[53]、偽譜法等。在實(shí)際應(yīng)用中通常根據(jù)實(shí)際情況綜合各種方法達(dá)到計(jì)算目的。

1.10.4 迭代反演

迭代反演通過(guò)求解震源和介質(zhì)參數(shù)的方程來(lái)獲得模型參數(shù)，將每次反演出的數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)相比較，不斷迭代直到兩者誤差達(dá)到精度要求。反演方法可以分為兩類(lèi)：第一類(lèi)是基于算子的線(xiàn)性反演方法，包括代數(shù)重建技術(shù)(ART)[54]、聯(lián)合迭代重建法(SIRT)[55]、奇異值分解法(SVD)[56]、最速下降法、共軛梯度法[57]、最小二乘法(如最小二乘正交分解法(LSQR)[58]、極小殘量法(GMRES)[59]和雙穩(wěn)定共軛梯度法(BICGSTAB)[60]等)、阻尼最小二乘法[33]等；第二類(lèi)是基于模型的完全非線(xiàn)性反演方法，包括遺傳算法[61]、模擬退火法[62]和神經(jīng)網(wǎng)絡(luò)法等。非線(xiàn)性反演方法雖然對(duì)于復(fù)雜問(wèn)題效果更為明顯，但是計(jì)算效率低，因而目前在地震體波層析成像中，線(xiàn)性反演方法應(yīng)用更為廣泛。Liu等[63]通過(guò)結(jié)合線(xiàn)性和非線(xiàn)性方法，提出了Two-step MCMC方法，效果很好。

趙昭等[64]同時(shí)使用反射和折射2種走時(shí)資料，在改善層析成像水平和垂向分辨率上取得了較好的效果。崔巖等[65]提出了用于反演的初至波走時(shí)層析成像的Tikhonov正則化模型和一種層析成像的梯度優(yōu)化算法，數(shù)值實(shí)驗(yàn)表明，該方法適用于速度差別的任意模型，且求解穩(wěn)定易于實(shí)現(xiàn)。Pegah[66]將折射波層析成像與面波多通道分析相結(jié)合，進(jìn)行近地表巖土工程研究，取得了較好的效果。

2 發(fā)展前景

目前，工程地震折射波法更多是用來(lái)對(duì)未知地區(qū)做作初步調(diào)查，在極淺層(即勘探深度為幾米至幾十米)和硬巖勘探中有著較大優(yōu)勢(shì)，如要成為被公認(rèn)的勘探技術(shù)，仍需要結(jié)合各個(gè)反演解釋方法的利弊發(fā)展成一套綜合的解釋方法，使之適應(yīng)大多數(shù)實(shí)際情況，而不是分類(lèi)列舉各種巧妙技術(shù)。這些方法的發(fā)展表明折射波法的發(fā)展趨勢(shì)是勘探深度不斷增加、應(yīng)用范圍不斷推廣和復(fù)雜程度不斷提高?？碧教攸c(diǎn)趨向于近地表非均質(zhì)極強(qiáng)、各向異性明顯、地形起伏大、勘探深度小和精度要求高，要求工程地震折射波法有更高的分辨率、信噪比、保真度。

由于淺層橫波勘探分辨率高，對(duì)地層結(jié)構(gòu)成像效果好但不易激發(fā)，而縱波勘探深度大且易于激發(fā)，但其成像效果差。因此，利用二者特點(diǎn)，研究多波聯(lián)合處理解釋對(duì)增加工程地震折射波勘探的精度和分辨率，以及減少反演結(jié)果的不確定性具有意義重大。

以折射層析成像研究為代表的解釋方法研究和觀(guān)測(cè)系統(tǒng)研究，是增加多層折射介質(zhì)成像精度的重點(diǎn)方向。

3 結(jié)論

本文概述了現(xiàn)有的幾種折射波反演解釋方法，當(dāng)探測(cè)深度較淺(5 m左右)且分界面足夠平整、界面傾角較小時(shí)，最方便的解釋方法是截距時(shí)間法；當(dāng)勘探目標(biāo)的深度達(dá)到25 m時(shí)，t0差數(shù)法最為適用；對(duì)于勘探目標(biāo)埋深超過(guò)25 m時(shí)，為了不抹平界面的不規(guī)則性、避免產(chǎn)生虛假的折射速度，應(yīng)當(dāng)使用廣義互換法；波前法、哈萊斯法、延遲時(shí)間法等或者原理與之類(lèi)似，或者效率和方便程度等綜合效果不如以上3種方法。以上的常規(guī)折射波解釋方法都需要求取正向和反向的旅行時(shí)之和，且實(shí)際數(shù)據(jù)存在大量噪音，人工提取走時(shí)存在困難，準(zhǔn)確率不夠高；而折射波褶積成像避免了拾取走時(shí)的過(guò)程，節(jié)約時(shí)間且效果較好。隨著勘探精度要求越來(lái)越高，使用折射波走時(shí)層析成像技術(shù)可以滿(mǎn)足近地表介質(zhì)速度縱橫向變化的地層、大傾角地層、隱伏層、界面起伏層等情況。

[1] 北京鈾礦地質(zhì)研究所淺層地震組. 淺層地震探測(cè)方法與技術(shù)[M]. 北京：原子能出版社，1982:41-84.

Shallow Seismic Group of Beijing Institute of Uranium Geology. Methods and Techniques for Shallow Seismic Detection[M]. Beijing：Atomic Energy Publishing House, 1982: 41-84.

[2] 熊章強(qiáng)，方根顯. 淺層地震勘探[M]. 北京：地震出版社，2002:79-90.

Xiong Zhangqiang, Fang Genxian. Shallow Seismic Exploration[M]. Beijing: Seismological Press, 2002: 79-90.

[3] 葛雙成，李小平，邵長(zhǎng)云，等. 地震折射和電阻率法在水庫(kù)壩址勘察中的應(yīng)用[J]. 地球物理學(xué)進(jìn)展，2008，23(4): 1299-1303.

Ge Shuangcheng, Li Xiaoping, Shao Changyun, et al. Application of Seismic Refraction and Resistivity for Exploration of Reservoir Dam Site[J]. Progress in Geophysics, 2008, 23(4): 1299-1303.

[4] 喻岳鈺，曹代勇，邢春穎，等. 黃土塬區(qū)模型約束折射靜校正技術(shù)應(yīng)用研究[J]. 地球物理學(xué)進(jìn)展，2016，31(3): 1373-1380.

Yu Yueyu, Cao Daiyong, Xing Chunying, et al. Application of the Technology of Model Constraint Refraction Statics on Loess Area[J]. Progress in Geophysics, 2016, 31(3): 1373-1380.

[5] Steven D S, Jeffery J N, Seth W, et al. Using Near-Surface Seismic Refraction Tomography and Multichannel Analysis of Surface Waves to Detect Shallow Tunnels: A Feasibility Study[J]. Journal of Applied Geophysics, 2013, 99: 60-65.

[6] Rucker M L. Integrating Seismic Refraction and Sur-face Wave Data Collection and Interpretation for Geotechnical Site Characterization[C]// Conference on Applied Geophysics. Missouri: Geophysics, 2006.

[7] 王洪. 折射波法在堤壩病害探查中的應(yīng)用研究[D]. 長(zhǎng)沙：中南大學(xué)，2007.

Wang Hong. Study on the Application of Refraction Wave Method in Disease Detection of Embankments and Dams[D]. Changsha: Central South University, 2007.

[8]Thurber C, Ritsema J. Theory and Observations: Seismic Tomography and Inverse Methods[J]. Treatise on Geophysics, 2007, 1: 323-360.

[9] Sheehan, Jacob R, William E, et al. An Evaluation of Methods and Available Software for Seismic Refraction Tomography Analysis[J]. Journal of Environmental and Engineering Geophysics, 2005, 10(1): 21-34.

[10] Thornburgh H R. Wave-Front Diagrams in Seismic Interpretation[J]. American Association of Petroleum Geologists Bulletin, 1930, 14: 185-200.

[11] Gardner L W. AnAreal Plan of Mapping Subsurface Structure by Refraction Shooting[J]. Geophysics, 1939, 4(4): 247.

[12] Ewing M,Woollard, Vine A C, et al. Geophysical Investigations in the Emerged and Submerged Atlantic Coastal Plain: Part IV: Cape May, New Jersey, Section[J]. Geological Society of America Bulletin, 1940, 51: 1821-1840.

[13] Hales F W. An Accurate Graphical Method for Inter-preting Seismic Refraction Lines[J]. Geophysical Prospecting, 1958, 6(3): 285-294.

[14]Гамбурцев Г A. Основы Сейсморазведки Гостоптех- Издаг Москва[M]. Матемсб: Гостоптехиздат, 1959.

[15]Hagedoorn J G. The Plus-Minus Method of Inter-preting Seismic Refraction Sections[J]. Geophysical Prospecting, 1959, 7(2): 158-182.

[16]Hawkins L V. The Reciprocal Method of Routine Shallow Seismic Refraction Investigations[J]. Geophysics, 1961, 26(6): 806.

[17]Palmer D. An Introductionto the Generalized Reci-procal Method of Seismic Refraction Interpretation[J]. Geophysics, 1981, 46(11): 1508.

[18] Palmer D. Imaging Refractors with the Convolution Section[J]. Geophysics, 2001, 66: 1582-1589.

[19] Hearn T M, Clayton R W. Lateral Velocity Vari-ations in Southern California: II: Results for the Lower Crust fromPnWaves[J]. Bull Seism Soc Am, 1986, 76: 511-520.

[20]Baumgarte J. Konstruktive Darstellungvon Seismi-chen Horizon Tenunter Berock Sichtigung Der Strahlen Brechung ImRaum[J]. Geophysical Prospecting, 1955, 3(2): 126-162.

[21]Schenck F L. Refraction Solutions by Wavefront Targeting: Abstract[J]. Aapg Bulletin, 1963, 47(9): 1775-1776.

[22]Adachi R. On aProof of Fundamental Formura Concerning Refraction Method of Geophysical Prospecting and Some Remarks[J]. Kumamoto Journal of Science L Ser A: Mathematics Physics & Chemistry, 1954, 2: 18-23.

[23]Mota L. Determination of Dips and Depths of Geological Layers by the Seismic Refraction Method[J]. Geophysics, 1954, 19(2): 242.

[24] Johnson S H. Interpretation of Split-Spread Refrac-tion Data in Terms of Plane Dipping Layers[J]. Geophysics, 1976, 41(3): 418.

[25] 王奇. 淺層地震氣折射波法在工程應(yīng)用中的研究[D]. 成都：成都理工大學(xué)，2010.

Wang Qi. Study of Shallow Seismic Gas Refraction Wave Method in Engineering Application[D]. Chengdu: Chengdu University of Technology, 2010.

[26] 陳滋康. 用共軛點(diǎn)法解釋折射時(shí)距曲線(xiàn)及其計(jì)算程序[J]. 物探與化探，1986，10(5): 350-355.

Chen Zikang. Explain Refraction Time Distance Curve and Its Calculation Procedure Using Conjugate Point Method[J]. Geophysical and Geochemical Exploration, 1986, 10(5): 350-355.

[27] 熊章強(qiáng)，方根顯. 淺層折射波法在隧道工程地質(zhì)調(diào)查中的應(yīng)用[J]. 水文地質(zhì)工程地質(zhì)，2001，28(1): 63-66.

Xiong Zhangqiang, Fang Genxian. Application of Refractive Method in Tunnel-Engineering Exploration[J]. Hydrogeology and Engineering Geology, 2001, 28(1): 63-66.

[28] 彭驍，王鵬. 淺層地震折射波法在候家梁隧道地質(zhì)勘查中的應(yīng)用[J]. 工程地球物理學(xué)報(bào)，2013，10(5):631-636.

Peng Xiao, Wang Peng. The Application of Shallow Seismic Refraction to Houjialiang Tunnel Exploration[J]. Chinese Journal of Engineering Geophysics, 2013,10(5): 631-636.

[29] 包勛，湯浩，朱照拔. 綜合物探技術(shù)在探測(cè)隱伏斷層中的應(yīng)用[J]. 人民珠江，2016，37(12): 29-32.

Bao Xun, Tang Hao, Zhu Zhaoba. Application of Comprehensive Geophysical Prospecting Method to Buried Fault Detection[J]. Pearl River, 2016, 37(12): 29-32.

[30] 張恒超，李學(xué)聰. 模型法、擴(kuò)展廣義互換法(EGRM)、沙丘曲線(xiàn)法靜校正對(duì)比研究及在ZGE沙漠資料處理中的應(yīng)用[J]. 地球物理學(xué)進(jìn)展，2010，25(6)：2193-2198.

Zhang Hengchao, Li Xuecong. Study of Static Corrections by the Model Method, Extended Generalized Reciprocal Method (EGRM) and Sand Dune Curve Method and Their Application to Seismic Data Processing in the ZGE Area[J]. Progress in Geophysics, 2010, 25(6): 2193-2198.

[31] Tak M L. Controls of Traveltime Data and Problems of the Generalized Reciprocal Method[J]. Geophysics, 2003, 68(5): 1626-1632.

[32]Sj?gren B. A Brief Study of Applications of the Generalized Reciprocal Method and of Some Limitations of the Method[J]. Geophysical Prospecting, 2000, 48(5): 815-834.

[33] Franco R D. Multi-Refractor Imaging with Stacked Refraction Convolution Section[J]. Geophysical Prospecting, 2005, 53(3): 335-348.

[34]Aki K, Lee W H K. Determination of Three‐Dimensional Velocity Anomalies Under a Seismic Array Using First P Arrival Times from Local Earthquakes: 1: A Homogeneous Initial Model[J]. Journal of Geophysical Research Atmospheres, 1976, 81(23): 4381-4399.

[35] 景月紅. 地震初至波走時(shí)層析成像與近地表速度建模[D]. 西安：長(zhǎng)安大學(xué)，2009.

Jing Yuehong. Seismic First Break Travel-Time Tomography and Its Application in Near-surface Velocity Model Building[D]. Xi’an: Chang’an University, 2009.

[36] 陳國(guó)金，高志凌，吳永栓，等. 井間地震層析成像中自動(dòng)生成初始速度模型的方法研究[J]. 石油物探，2005，44(4):339-342.

Chen Guojin, Gao Zhiling, Wu Yongshuan, et al. Automatically Building the Initial Velocity Model in the Seismic Crosshole Tomography[J]. Geophysical Prospecting for Petroleum, 2005, 44(4): 339-342.

[37] 段心標(biāo)，金維浚. 井間地震層析成像初始速度模型[J].地球物理學(xué)進(jìn)展，2007，22(6):1831-1835.

Duan Xinbiao, Jin Weijun. Initial Velocity Model of Crosswell Seismic Tomography[J]. Progress in Geophysics, 2007, 22(6): 1831-1835.

[38] 劉玉柱，丁孔蕓，董良國(guó). 初至波走時(shí)層析成像對(duì)初始模型的依賴(lài)性[J]. 石油地球物理勘探，2010，45(4):502-511.

Liu Yuzhu, Ding Kongyun, Dong Liangguo. The Dependence of First Arrival Travel Time Tomography on Initial Model[J]. Oil Geophysical Prospecting, 2010, 45(4): 502-511.

[39] Engdahl E R. Relocation of Intermediate Depth Earth-quakes in the Central Aleutians by Seismic Ray Tracing[J]. Nature, 1973, 245(141): 23-25.

[40] Julian B R,Gubbins D. Three Dimensional Seismic Ray Tracing[J]. J Geophysics,1977, 43: 95-113.

[41] Nakanishi I, Yamaguchi K. A Numerical Experiment on Nonlinear Image Reconstruction from First-Arrival Times for Two-Dimensional Island Arc Structure[J]. Earth Planets and Space, 1986, 34(2):195-201.

[42] Asakawa E,Kawanaka T. Seismic Ray Tracing Using Linear Traveltime Interpolation[J]. Geophysical Prospecting, 2010, 41(1): 99-111.

[43] Moser T J. Shortest Path Calculation of Seismic Rays[J]. Geophysics, 1991, 56(1): 59-67.

[44]Cheng N. Minimum Traveltime Calculation in 3-D Graph Theory[J]. Geophysics, 1996, 61(6): 1895-1898.

[45]Cardarelli E,Cerreto A. Ray Tracing in Elliptical Anisotropic Media Using the Linear Traveltime Interpolation (LTI) Method Applied to Traveltime Seismic Tomography[J]. Geophysical Prospecting, 2002, 50(1): 55-72.

[46] 張東，謝寶蓮，楊艷，等. 一種改進(jìn)的線(xiàn)性走時(shí)插值射線(xiàn)追蹤算法[J]. 地球物理學(xué)報(bào)，2009，52(1):200-205.

Zhang Dong, Xie Baolian, Yang Yan, et al. A Ray Tracing Method Based on Improved Linear Traveltime Interpolation[J]. Chinese Journal of Geophysics, 2009, 52(1): 200-205.

[47] Vidale J. Finite-Difference Calculation of Travel Time in Three Dimensions[J]. Geophysics, 1990, 55: 521-526.

[48]Tsai Y H R, Cheng L T,Osher S, et al. Fast Sweeping Algorithms for a Class of Hamilton-Jacobi Equations[J]. SIAM Journal on Numerical Analysis, 2004, 41(2): 659-672.

[49] Sethian J A. A Fast Marching Level Set Method for Monotonically Advancing Fronts[J]. Proceedings of the National Academy of Science of the United States of America, 1996, 93(4): 1591-1595.

[50] Sethian J A. Level Set Methods and Fast Marching Method[M]. Cambridge: Cambridge University Press, 1999: 400.

[51] 王飛. 跨孔雷達(dá)走時(shí)層析成像反演方法的研究[D]. 長(zhǎng)春：吉林大學(xué)，2014.

Wang Fei. Research onCrosshole Radar Traveltime Tomography[D]. Changchun: Jilin University, 2014.

[52]Yee K S. Numerical Solution of Initial Boundary Value Problems Involving Maxwell Equations in Isotropic Media[J]. IEEE Trans Antennas Propagat, 1966, 14(3): 302-307.

[53] Yoon Y S, Yeh W G. Parameter Identification in an Inhomogeneous Medium with the Finite-Element Method[J]. Society of Petroleum Engineers Journal, 1976, 16(4): 217-226.

[54] Herman G T, Brouw W N. Book-Review: Image Re-construction from Projections: The Fundamentals of Computerized Tomography[J]. Space Science Reviews, 1980, 32(4): 56-61.

[55]Humphreys E, Clayton R W. Adaptation of Back Projection Tomography to Seismic Travel Time Problems[J]. Journal of Geophysical Research Solid Earth, 1988, 93(B2): 1073-1085.

[56] Golub G H, Reinsch C. Singular Value Decomposi-tion and Least Squares Solutions[J]. Numerische Mathematik, 1970, 14(5): 403-420.

[57] Scales J A. Tomographic Inversion via the Conjugate Gradient Method [J]. Geophysics,1987, 52(2): 179.

[58] Paige C C, Saunders M A. LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares[J]. ACM Transactions on Mathematical Software, 1982, 8(1): 43-71.

[59] Saad Y, Schultz M H. GMRES: A Generalized Mi-nimal Residual Algorithm for Solving Nonsymmetric Linear Systems[J]. SIAM Journal on Scientific and Statistical Computing, 1986, 7(3): 856-869.

[60]Sleijpen G L G, Fokkema D R. BiCGstab(ell) for Linear Equations Involving Unsymmetric Matrices with Complex Spectrum[J]. Electronic Transactions on Numerical Analysis Etna, 1993, 1(11): 1160.

[61] Amoon C J, Vidale J E. Tomography Without Rays[J]. Bulletin of the Seismological Society of America, 1993, 83(2): 509-528.

[62] Kirkpatrck S, Jr G C, Vecchi M P, et al. Optimi-zationby Simulated Annealing[J]. Science, 1983, 220: 671-680.

[63] Liu J, Braile L. Two-Step MCMC Seismic Tomog-raphy[C]// AGU Fall Meeting. San Francisco: AGU, 2004.

[64] 趙昭，趙志新，徐紀(jì)人，等. 高精度人工地震波層析成像技術(shù)勘探地下構(gòu)造[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2011，41(增刊1):362-367.

Zhao Zhao, Zhao Zhixin, Xu Jiren, et al. Surveying Structures with High-Resolution Seismic Tomography[J]. Journal of Jilin University (Earth Science Edition), 2011, 41(Sup. 1): 362-367.

[65] 崔巖，王彥飛. 基于初至波走時(shí)層析成像的Tikhonov正則化與梯度優(yōu)化算法[J]. 地球物理學(xué)報(bào)，2015，58(4):1367-1377.

Cui Yan, Wang Yanfei. Tikhonov Regularization and Gradient Descent Algorithms for Tomography Using First-Arrival Seismic Traveltimes[J]. Chinese Journal of Geophysics, 2015, 58(4): 1367-1377.

[66] Pegah E, Liu H. Application of Near-Surface Seismic Refraction Tomography and Multichannel Analysis of Surface Waves for Geotechnical Site Characterizations: A Case Study[J]. Engineering Geology, 2016, 208: 100-113.