工程地震折射波解釋方法研究進(jìn)展

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

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

0 引言

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

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

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

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

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

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

1 方法原理

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

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

1.1 波前法

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

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

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)之間的距離除以兩點(diǎn)間重構(gòu)的波前時(shí)間差值。

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

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

1.2 延遲時(shí)間法

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

1.3 截距時(shí)間法

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

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

(2)

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

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

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

1.4 哈萊斯法

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

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

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

1.5 共軛點(diǎn)法

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

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

1.6 t0差數(shù)法

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

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

t+=tO1ABS+tO2DCS-tO1ADO2，

(3)

t-=tO1ABS-tO2DCS+tO1ADO2。

(4)

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

(5)

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

(6)

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

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

互換法又稱表層剝?nèi)シā⑷ケ韺臃?，?yīng)用起來最為簡潔并能估算折射層結(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)用于隧道探測中，與t0差數(shù)法進(jìn)行比較取得了相似的效果(圖8、9)。

(7)

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

(8)

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

(9)

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

據(jù)文獻(xiàn)[29]。圖6 某水庫折射波綜合時(shí)距曲線及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 廣義互換法

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

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

(10)

(11)

廣義互換法大致可分為3個過程：1)求取速度分析函數(shù)；2)求取時(shí)深偏移函數(shù)；3)選取最佳XY值。該方法通過選取最佳XY值并帶入步驟1)和2)中的函數(shù)可求出折射層速度和深度。當(dāng)XY=0時(shí)，廣義互換法即為通常意義上的互換法，可能會導(dǎo)致虛假折射速度大量出現(xiàn)，并使得有起伏折射面的細(xì)節(jié)被圓滑掉。該方法處理傾角小于20°的界面效果一般，而應(yīng)用于大傾角、存在橫向速度變化或隱蔽帶時(shí)效果較好[31]。Sj?gren[32]嘗試將廣義互換法用于淺層工程調(diào)查中并與哈萊斯法進(jìn)行比較，驗(yàn)證了其適應(yīng)性，并總結(jié)了勘查過程中可能會遇到的問題(圖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ò)散的影響，對折射波褶積成像法做出了改進(jìn)。該方法不需要拾取初至?xí)r間，可以直接將同一點(diǎn)處的正向和反向走時(shí)數(shù)據(jù)進(jìn)行褶積得到走時(shí)數(shù)據(jù)的交叉時(shí)，再利用速度求出該點(diǎn)的垂向深度。Franco[33]提出了利用褶積和互相關(guān)運(yùn)算對數(shù)據(jù)疊加來提高信噪比的方法，并對多個折射界面成像，是常規(guī)解釋方法的一大突破(圖12)。

1.10 折射波層析成像

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

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

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

1.10.1 初至走時(shí)提取

地震勘探中把經(jīng)過地下介質(zhì)首先到達(dá)檢波器的波稱為初至波，這些波的旅行時(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í)測數(shù)據(jù)Fig.12 Application to real data of Clusone basin area in Northern Italy

1.10.2 速度模型建立

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

1.10.3 正演

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

1)射線追蹤法

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

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

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

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

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

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

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

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

1.10.4 迭代反演

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

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

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

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

3 結(jié)論

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

[1] 北京鈾礦地質(zhì)研究所淺層地震組. 淺層地震探測方法與技術(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] 葛雙成，李小平，邵長云，等. 地震折射和電阻率法在水庫壩址勘察中的應(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]. 長沙：中南大學(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í)距曲線及其計(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ù)在探測隱伏斷層中的應(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)、沙丘曲線法靜校正對比研究及在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]. 西安：長安大學(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] 陳國金，高志凌，吳永栓，等. 井間地震層析成像中自動生成初始速度模型的方法研究[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] 劉玉柱，丁孔蕓，董良國. 初至波走時(shí)層析成像對初始模型的依賴性[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)的線性走時(shí)插值射線追蹤算法[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]. 長春：吉林大學(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.