工程地震折射波解釋方法研究進展

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

吉林大學地球探測科學與技術學院，長春 130026

0 引言

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

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

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

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

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

本文介紹了常用折射波解釋的方法原理，著重介紹了折射層析成像方法，并闡述了這些方法各自的優缺點及適用性，同時對國內外工程地震折射波解釋方法的發展歷程及現狀進行闡述。

1 方法原理

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

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

1.1 波前法

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

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

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

(1)

式中：xR表示點R橫坐標；zR表示點R縱坐標。則可據此找出折射界面的位置。折射層的速度可由重構的波前傳播時間估算，即用界面上兩點之間的距離除以兩點間重構的波前時間差值。

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

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

1.2 延遲時間法

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

1.3 截距時間法

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

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

(2)

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

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

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

1.4 哈萊斯法

哈萊斯法通過作圖計算臨界折射時間，再求取折射界面深度進行構造解釋[13]。

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

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

1.5 共軛點法

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

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

1.6 t0差數法

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

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

t+=tO1ABS+tO2DCS-tO1ADO2，

(3)

t-=tO1ABS-tO2DCS+tO1ADO2。

(4)

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

(5)

S點處界面埋深h為

(6)

依次類推可獲得測線上全部點的界面埋深。

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

據文獻[25]。圖4 遼寧鞍山境內某隧道哈萊斯法解釋圖Fig.4 Interpretation of a tunnel in Anshan, Liaoning Province

圖5 t0差數法原理圖Fig.5 Principle of t0difference method

1.7 互換法

互換法又稱表層剝去法、去表層法，應用起來最為簡潔并能估算折射層結構和速度情況[16]。

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

(7)

則A點下折射界面深度為

(8)

根據幾何方法，可得D點延遲時間：

(9)

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

據文獻[29]。圖6 某水庫折射波綜合時距曲線及t0差數法解釋剖面結果圖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 廣義互換法

廣義互換法是互換法的拓展，也需要正反兩個時距曲線，其計算速度和時間深度的公式與傳統互換法相似，但增加了一個補償項進行偏移校正處理[17]。

如圖10所示，X、Y為接收點，如果地表上的D點與接收點不重合，則D點延遲時為

(10)

(11)

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

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

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

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

1.9 褶積成像法

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

1.10 折射波層析成像

層析成像技術最先用于醫學，隨后拓展到地球物理學界。Aki等[34]最早提出地震走時層析成像方法，并應用于美國某地區的地殼結構探測。Hearn等[19]利用6年的天然地震數據首次實現了折射波層析成像。

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

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

1.10.1 初至走時提取

地震勘探中把經過地下介質首先到達檢波器的波稱為初至波，這些波的旅行時包含了淺層速度信息，具有易于獲取、能量強、可追蹤性好的優點[35]。初至波走時的準確拾取構成了折射層析成像的基礎和前提。

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

據文獻[33]。圖12 意大利北部克盧索內盆地實測數據Fig.12 Application to real data of Clusone basin area in Northern Italy

1.10.2 速度模型建立

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

1.10.3 正演

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

1)射線追蹤法

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

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

彎曲法是根據費馬原理，通過調整假定參數以達到震源與檢波點間走時最小的一種方法。這種方法雖然不存在陰影區，但用于局部復雜的介質模型會收斂到局部最小值，且計算效率和精度較低。

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

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

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

2)波動方程數值模擬

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

1.10.4 迭代反演

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

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

2 發展前景

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

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

以折射層析成像研究為代表的解釋方法研究和觀測系統研究，是增加多層折射介質成像精度的重點方向。

3 結論

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

[1] 北京鈾礦地質研究所淺層地震組. 淺層地震探測方法與技術[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] 熊章強，方根顯. 淺層地震勘探[M]. 北京：地震出版社，2002:79-90.

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

[3] 葛雙成，李小平，邵長云，等. 地震折射和電阻率法在水庫壩址勘察中的應用[J]. 地球物理學進展，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] 喻岳鈺，曹代勇，邢春穎，等. 黃土塬區模型約束折射靜校正技術應用研究[J]. 地球物理學進展，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] 王洪. 折射波法在堤壩病害探查中的應用研究[D]. 長沙：中南大學，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] 王奇. 淺層地震氣折射波法在工程應用中的研究[D]. 成都：成都理工大學，2010.

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

[26] 陳滋康. 用共軛點法解釋折射時距曲線及其計算程序[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] 熊章強，方根顯. 淺層折射波法在隧道工程地質調查中的應用[J]. 水文地質工程地質，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] 彭驍，王鵬. 淺層地震折射波法在候家梁隧道地質勘查中的應用[J]. 工程地球物理學報，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] 包勛，湯浩，朱照拔. 綜合物探技術在探測隱伏斷層中的應用[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] 張恒超，李學聰. 模型法、擴展廣義互換法(EGRM)、沙丘曲線法靜校正對比研究及在ZGE沙漠資料處理中的應用[J]. 地球物理學進展，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] 景月紅. 地震初至波走時層析成像與近地表速度建模[D]. 西安：長安大學，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] 段心標，金維浚. 井間地震層析成像初始速度模型[J].地球物理學進展，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] 劉玉柱，丁孔蕓，董良國. 初至波走時層析成像對初始模型的依賴性[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]. 地球物理學報，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]. 長春：吉林大學，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]. 吉林大學學報(地球科學版)，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] 崔巖，王彥飛. 基于初至波走時層析成像的Tikhonov正則化與梯度優化算法[J]. 地球物理學報，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.