多道瞬態面波在復雜地形條件下巖層劃分中的應用研究

席超強， 張平松， 李建寧， 丁美青

(安徽理工大學 地球與環境學院， 淮南 232001)

多道瞬態面波在復雜地形條件下巖層劃分中的應用研究

席超強， 張平松， 李建寧， 丁美青

(安徽理工大學 地球與環境學院， 淮南 232001)

在復雜地表條件下，利用瑞雷面波進行勘探時，需要考慮地形對瑞雷波場傳播地影響。針對傾斜地表兩層地質模型，通過設計不同展布方向的檢波器排列進行瑞雷波信號采集，利用高精度交錯網格有限差分法，討論經過快速傅里葉變換后的地震記錄F-V頻率譜，通過對比信號頻率譜與理論多階頻散曲線，分析在不同地形條件下瑞雷波傳播能量的變化特征，并提取疊加頻散曲線；運用遺傳算法對采集到的地震記錄進行基階與高階瑞雷波聯合反演得到目的層厚與初始模型層厚進行誤差對比；并對理論正演初始模型的頻散曲線與各地質模型利用F-K法提取的頻散曲線進行全局擬合誤差計算。分析表明：傾斜地形進行多道瞬態瑞雷面波勘探，相對于下部激發、上部激發會產生高階頻散趨勢，因此對下部激發的地震記錄進行反演產生的厚度誤差相對較小。

瑞雷波； 傾斜； 地層劃分； 聯合反演

0 引言

近年來，瑞雷波在淺層地質勘查和環境調查方面地應用越來越廣泛，然而在實際勘察工作中往往會遇到起伏的地形條件，在這種復雜的地形條件中采集到的瑞雷波具有不同的傳播特征。為研究瑞雷波在層狀介質中的頻散特征，大量學者進行了深入的正演研究。從Haskell[1]提出計算瑞雷波頻散曲線的Haskell算法開始，許多學者提出多種切實可行的頻散曲線數值運算方法，Abo-Zena[2]通過反對稱矩陣的循環運算得到瑞雷波頻散方程；張碧星[3]采用柱坐標系研究瑞雷波頻散問題，簡化計算；凡友華[4]在柱坐標系下基于快速矢量傳遞矩陣提出快速矢量傳遞算法；潘佳鐵[5]利用快速廣義反射投射系數法計算瑞雷面波的群速度，并用數值實驗的方法合成了兩個典型模型的群速度頻散曲線；袁臘梅[6]通過將廣義反射-投射系數算法無量綱化得到瑞雷波頻散函數的公式體系。

隨著山區、滑坡等地形復雜地區地震勘探的興起，大量專家學者基于彈性交錯網格有限差分法，對起伏地形條件下的瑞雷波地震波場進行正演模擬。Hestholin等[7]對起伏地表下彈性波傳播有限差分數值模擬進行研究；Robertsson[8]給出了含起伏地表彈性和粘彈性介質中地震波傳播的有限差分模擬方法；張大洲[9]應用2×12階高精度交錯網格有限差分法,建立了震源位于自由表面時模擬瑞雷波的邊界條件，對均勻半空間模型模擬；熊章強[10]模擬分析小凸起小凹陷模型條件下的瑞雷波場傳播特征；江凡[11]模擬起伏地表條件下的瑞雷波全波場正演，通過坐標變換的方式實現起伏地表自由邊界條件；琚長輝[12]進行了傾斜和塌陷地形的彈性波場交

錯網格有限差分數值模擬。以上對瑞雷波在數值模擬上的分析主要集中在完善對瑞雷波有限差分模擬地實現與瑞雷波傳播受地形影響方面，并未對地形、對瑞雷波頻散特征影響進行分析。

研究瑞雷波面波在傾斜起伏地形條件下的傳播特征和頻散曲線，具有重要的現實意義，通過交錯網格有限差分法，模擬傾斜地表下瑞雷波傳播，并對地表凸起和凹陷以及平直地形模式進行波場模擬，在此基礎上分析傾斜地形下瑞雷波頻散特征，討論“之”形特征，根據遺傳算法對模擬得到的頻散曲線進行反演，并與模型參數進行對比，獲得了傾斜地形對巖土分層的影響。

1 數值模擬實驗

根據傾斜地表的構造簡化為雙層介質模型，在模型表面進行瑞雷波數值模擬，通過頻散特征對地層進行劃分并與排列中點的層厚進行對比。我們設計四種雙層均勻介質模型如圖1所示，凹陷和凸起地形的寬度也是100 m，運用交錯網格有限差分法使用彈性波動方程，模擬瑞雷波在模型中傳播。模型的尺度為100 m×100 m，網格點數為1 000×1 000，模型介質參數為：上部地層用綠色表示，下部地層用藍色表示。表1為兩層傾斜地質模型的物性參數。

表1 兩層傾斜地質模型參數

數值模擬使用的震源為主頻20 Hz 的Ricker子波。如圖1所示，每個模型從左到右依次有三個炮點，炮點1在模型上部激發，x=20 m處激發，26個檢波器，道間距為2 m，檢波器位置x=22 m～72 m；炮點3在模型下部激發，x=74 m處激發，26個檢波器，道間距為2 m，檢波器位置x=22 m～72 m。采樣間隔為0.2 ms，采樣長度為0.5 s。炮點2在模型中部x=47 m處激發，檢波器沿地層垂直于炮點1、炮點3測線的方向展開，使用26個檢波器，道間距為2 m。

平直、凸起、凹陷三種等厚表層模型的表層豎直深度均為10 m，在傾斜地表上部和下部進行震源激發；由于表層具有傾斜角度和起伏變化，排列中點相對于地層界面的距離不等于豎直方向深度。表2為四類模型三種激發方式排列中，點到沿垂直排列切線方向到地層分界面的距離。

表2 兩層傾斜地質模型參數

2 頻散特征分析

2.1 頻散譜分析

對高信噪比的瑞雷波地震記錄采用有效的信號技術進行處理，實現瑞雷波地震波場分離提高頻散曲線的提取精度和可靠性。目前主要的頻散曲線提取算法有相移法、Prony法和拉冬變換算法等。筆者使用f-k域變換法對瑞雷波地震記錄進行快速傅里葉變換，利用Vr=f/k的關系轉換到f-v域進行頻散曲線提取，可以實現多模式頻散曲線的提取。

2.1.1 單傾直線等厚表層模型

從圖2(a)可見，基階波在低頻占主導地位，在高頻處于次要地位，在10 Hz～30 Hz頻帶范圍內基階波能量為主導，在35 Hz～45 Hz的頻帶范圍內第二高階能量為主導，45 Hz～60 Hz頻帶范圍內第一高階能量為主導。此時地震記錄頻散曲線為基階面波與各高階瑞雷波共同耦合作用的結果。拾取的瑞雷波會出現“之”字形拐點。

圖2(b)頻譜是根據震源為炮點2在檢波點沿模型表面水平方向展開得到的模型模擬地震記錄經過傅里葉變換得到。從圖2(b)可見，在0 Hz～8 Hz的頻帶范圍內，如圖中箭頭所示，相速度值急劇增加，與實際模型并不相符，因此在拾取頻散曲線速度的時候在低頻范圍內需要謹慎，不得根據能量“脊”線盲目拾取，在8 Hz～30 Hz頻帶范圍內，基階能量占主導地位，在30 Hz～60 Hz，能量介于基階能量和第一高階能量之間，此時頻散曲線應該是基階波與第一高階波相互疊加共同耦合的結果，此時疊加后的頻散曲線也會出現“之”字形拐點。通過出現“之”字形拐點，與瑞雷波頻散有關的參數有地層密度、地層厚度、縱橫波速度、密度、泊松比。模型的兩層地層的密度和泊松比均為0.45 g/cm3和1.8 g/cm3，說明利用“之”字形的這一特征可以為地層的速度分層提供參考依據。通過給定的頻率和速度可計算得到確定的深度值，經計算發生拐點的深度為6.2 m，與實際模型的10 m并不相符，說明不能直接根據跳點處的深度，來確定巖性分層的深度，這與文獻[13]關于“之”字形的描述相同。

圖1 兩層傾斜地形模型Fig.1 Two inclined geological model(a)單傾平直等厚表層模型；(b)單傾平直不等厚表層模型;(c)單傾凸起等厚表層模型;(d)單傾凹陷等厚表層模型

圖2 單傾直線等厚表層模型 瑞雷波F-V頻率譜Fig.2 F-V Rayleigh wave frequency spectrum of pour equal thickness single linear surface model(a)炮點1頻譜；(b)炮點2頻譜；(c)炮點3頻譜

從圖2(c)可見，在8 Hz～30 Hz頻帶范圍內，基階能量占主導地位，在30 Hz～60 Hz，能量第一高階能量占主導地位。

通過對比炮點1和炮點3位于模型表層上、下位置的頻譜，上部激發的地震記錄頻散譜第二高階能量分離較清晰；下部激發的地震記錄頻散譜清新分離第一高階能量。炮點2頻譜低頻部分基階能量主導，高頻能量基階能量和第一高階共同耦合。

2.1.2 平直不等厚表層模型

從圖3(a)可見，5 Hz～32 Hz頻譜范圍內能量主要屬于基階能量，30 Hz～60 Hz能量為第一高階能量占主導地位。由圖3(b)可見，排列為水平方向順層激發，不受地形影響，在8 Hz～30 Hz頻帶范圍內，能量主要屬于基階能量，在30 Hz～60 Hz，能量范圍介于基階能量和第一高階能量之間，與前兩組水平激發模型頻譜規律相同。由圖3(c)可見，在5 Hz～35 Hz頻帶范圍內，基階能量占主導地位，在30 Hz～60 Hz，能量第一高階能量占主導地位，與第一組模型下部激發頻散特征相同，45 Hz～50 Hz頻帶范圍內面波能量受第一高階和第二高階能量共同耦合作用影響，有向第二高階分離的趨勢。

2.1.3 單傾凸起等厚表層模型

從圖4(a)可見，8 Hz～26 Hz頻譜范圍內能量主要屬于基階能量，30 Hz～38 Hz能量為第一高階和第二高階共同耦合結果，38 Hz～51 Hz能量為第一高階能量，與平直等厚模型頻散特征相似，疊加頻散曲線呈“基階-第二高階-第一高階”的趨勢。

從圖4(b)可見，排列為水平方向順層激發，不受地形影響，在8 Hz～30 Hz頻帶范圍內，能量主要屬于基階能量，在30 Hz～60 Hz，能量范圍介于基階能量和第一高階能量之間，與前兩組水平激發模型頻譜規律相同。

由圖4(c)可見，在8 Hz～30 Hz頻帶范圍內，基階能量占主導地位，在30 Hz～60 Hz頻帶范圍內，面波能量模式為第一高階，未出現趨向第二高階的特征，與第一組模型下部激發頻散特征相同。

2.1.4 單傾凹陷等厚表層模型

從圖5(a)可見，8 Hz～27 Hz和34 Hz～36 Hz頻譜范圍內能量主要屬于基階能量，28 Hz～32 Hz和35 Hz～47 Hz頻譜范圍能量為第一高階能量，48 Hz～60 Hz能量為第二高階能量，與平直等厚模型和凸起等厚模型頻散特征不同，疊加頻譜曲線呈“基階-第一高階-第二高階”的趨勢。

圖3 單傾平直不等厚表層模型 瑞雷波F-V頻率譜Fig.3 F-V Rayleigh wave frequency spectrum of pour straight single uniform thick surface model(a)炮點1頻譜；(b)炮點2頻譜；(c)炮點3頻譜

圖4 單傾凸起等厚表層模型 瑞雷波F-V頻率譜Fig.4 F-V Rayleigh wave frequency spectrum of single pour projections surface model(a)炮點1頻譜；(b)炮點2頻譜；(c)炮點3頻譜

圖5 單傾凹陷等厚表層模型 瑞雷波F-V頻率譜Fig.5 F-V Rayleigh wave frequency of single pour depression surface model spectrum(a)炮點1頻譜；(b)炮點2頻譜；(c)炮點3頻譜

從圖5(b)可見，排列為水平方向順層激發，不受地形影響。在8 Hz～30 Hz頻帶范圍內，能量主要屬于基階能量，在30 Hz～60 Hz，能量范圍介于基階能量和第一高階能量之間。與前三組水平激發模型頻譜規律相同。基于有限差分的瑞雷波彈性波動方程數值模擬，在相同參數的運算結果較為穩定。

從圖5(c)可見，在8 Hz～30 Hz頻帶范圍內，基階能量占主導地位，30 Hz～35 Hz面波能量呈現向第二高階發育的趨勢，在35 Hz～60 Hz頻率范圍內，第一高階能量占主導地位。

通過對四個模型的頻散曲線對比得到：①上部激發相對于下部激發和水平順層激發方式更易發生高階波；②在40 Hz以上的高頻部分，模型中基階波與高階波的速度差異不大，判斷頻散的階數需謹慎；③傾斜地形的凹陷凸起狀態影響高頻頻帶范圍內高階波的產生，凹陷地形，高頻范圍高階波較發育。

2.2 頻散擬合誤差分析

根據模型的各層介質深度、縱橫波速度、密度等參數，通過頻散曲線正演公式計算可得到理論多模式頻散曲線。由于在頻散譜中部分高階能量為兩種模式波的共同耦合作用，通過手動定義前疊加頻散曲線各部分的階數，將疊加頻散曲線與初始模型的各頻帶范圍的多模正演曲線進行擬合誤差運算。采用計算均方根誤差的方法，來衡量疊加頻散曲線與多模頻譜的偏差。均方根值越小，表明全局相對擬合精度越高。

表3為4種模型各炮點所計算的疊加頻散曲線與多模式理論頻散曲線的均方根誤差。

圖6 理論多模頻散與實測疊加頻散Fig.6 Multi-mode dispersion theory and measured superimposed dispersion(a)10 m厚度表層的理論3階頻散；(b)模態分離的疊加頻散與理論頻散

模型類型炮點1均方根誤差/%炮點2均方根誤差/%炮點3均方根誤差/%平直等厚表層模型7.23.604.10平直不等厚表層模型17.643.078.14凸起等厚表層模型9.053.844.42凹陷等厚表層模型9.693.776.22

通過對比模型各檢波點排列均方根誤差得到結論：①傾斜地形條件下，炮點3(下部激發)的數據記錄提取的頻散曲線與檢波點排列范圍中點的理論頻散線的的擬合誤差最小；②傾斜地形條件下，凸起地形的頻散擬合誤差小于凹陷地形。

3 層厚反演分析

本次建模使用四類簡化的傾斜雙層地質模型，使用沿地形傾向和順層水平展布的三種檢波器排列。多道瞬態瑞雷面波法對采集記錄進行處理得到頻散特征，視為得到所使用排列中心點深度范圍內的頻散特征，并根據頻散曲線進行反演得到與橫波速度有關的相應巖土體地質特征。

對疊加頻散進行反演，建立初始模型，兩層速度模型上層橫波速度為200 m/s，下層橫波速度為400 m/s，以地層深度為反演目標通過遺傳算法運算，遺傳算法反演參數：種群大小64，迭代次數30，交叉概率90%，變異概率2%，種子數80，指數拉伸函數，采用二進制編碼，獲得全局非線性最優解。表4為四種模型的不同激震位置的地震記錄。

將反演深度與目標層深度的兩者之間的差值與模型的界面深度做差值，計算各炮點反演深度值的絕對誤差(表5)。

將反演深度與目標層深度的兩者之間的差值與模型的界面深度做比值，計算各炮點反演深度值的相對誤差(表6)。

圖7為四種模型各采集排列方式發演深度與初始模型深度的對比圖。

通過對比遺傳算法的深度反演結果與模型實際參數的對比，得到以下結論：①傾斜地形條件下，下部激震的記錄提取的頻散曲線反演的深度相對于上部激震更接近初始模型；②等厚模型的反演深度比較，上部和下部激震方式的反演相對誤差：平直地表>凸起地表>凹陷地表；③不等厚表層模型的上、下層激震的反演深度均小于初始模型表層深度。

表4 模型地層反演深度

表5 模型地層反演深度絕對誤差

表6 模型地層反演深度相對誤差

圖7 模型反演層厚對比圖Fig.7 Comparison of layer thickness model inversion(a)單傾平直等厚表層模型;(b)單傾平直不等厚表層模型;(c)單傾凸起等厚表層模型;(d)單傾凹陷等厚表層模型

4 結語

通過對四種雙層均勻介質模型的數值模擬，討論了傾斜地表條件下利用多道瞬態瑞雷面波法探測巖土體深度的應用效果，得到以下結論：

1)傾斜地表相對水平地表對瑞雷波的傳播造成影響，在傾斜條件下激發產生的地震記錄更易產生基階-高階波的分離。

2)起伏地形影響瑞雷波的傳播，凸起地形下的勘探結果優于凹陷地形。

3)在傾斜地表模型的條件下，震源置于排列系統下方的深度勘探效果要優于上部激震方式，可能與上部產生第二高階波的分離有關。

綜上所述，在傾斜地表條件下，通過多道瞬態瑞雷面波采集的數據利用遺傳算法基階-高階聯合反演，可以獲得反映實際深度的結果，反演深度相對模型表層存在一定的誤差。

[1] HASKELL N A. The dispersion of surface waves on multi-layered media.[J]. Bulletin of the Seismological Society of America, 1953, 43(1):86-103.

[2] ABO-ZENA A. Dispersion function computations for unlimited frequency values[J]. Geophysical Journal International, 2007, 58(1):91-105.

[3] 張碧星, 蘭從慶. 分層介質中面波的能量分布[J]. 聲學學報:中文版, 1998(2):97-106.

ZHANG B X, LAN C Q. Energy distribution of surface waves in stratified media[J]. Acta Acustica,1998(2):97-106.(In Chinese)

[4] 凡友華, 劉家琦. 層狀介質中瑞雷面波的頻散研究[J]. 哈爾濱工業大學學報, 2001, 33(5):577-581.

FAN Y H, LIU J Q. Research on the dispersion of rayleigh waves in multilayered media[J]. Journal of Harbin Institute of Technology, 2001, 33(5):577-581. (In Chinese)

[5] 潘佳鐵, 吳慶舉, 李永華,等. 基于快速廣義反射透射系數方法的面波群速度計算[J]. 地球物理學進展, 2009, 24(6):2030-2035.

PAN J T, WU Q J, LI Y H, et al.Group velocities computation of surface waves based on the fast generalized R/T coefficient method[J].Progress in Geophysics，2009, 24(6):2030-2035. (In Chinese)

[6] 袁臘梅, 凡友華, 孫書榮. 廣義反射-透射系數算法的無量綱化[J]. 地震學報, 2009, 31(4):377-384.

YUAN L M, FAN Y H, SUN S R. Making the calculation of generalized reflection transmission coefficients dimensionless [J]. Acta Seismologica Sinica, 2009, 31(4):377-384. (In Chinese)

[7] HESTHOLM S, RUUD B. 2D finite-difference elastic wave modelling including surface topography 1[J]. Geophysical Prospecting, 1996, 42(5):371-390.

[8] ROBERTSSON J O A. A numerical free-surface condition for elastic finite-difference modeling in the presence of topography[J]. Geophysics, 1996, 61(6):1921-1934.

[9] 熊章強, 南坤, 唐圣松. 起伏地表二維瑞雷面波正演模擬及波場分析[J]. 物探化探計算技術, 2010, 32(5):470-475.

XIONG Z Q, NAN K, TANG S S.Forward modeling and wave field analysis of 2D Rayleigh wave based on irregular topography [J].Computing Techniques For Geophysical And Geochemical Exploration,2010, 32(5):470-475. (In Chinese)

[10] 江凡, 林智勇. 瑞雷波在滑動面探測中的應用效果討論[J]. 福州大學學報:自然科學版, 2014, 42(2):299-304.

JIANG F, LIN Z Y. Discussion of the application effect of Rayleigh wave in the detection of sliding surface[J]. Journal of Fuzhou University( Natural Science Edition),2014, 42(2):299-304. (In Chinese)

[11] 張大洲, 熊章強, 顧漢明. 高精度瑞雷波有限差分數值模擬及波場分析[J]. 地球物理學進展, 2009, 24(4):1313-1319.

ZHANG D Z, XIONG Z Q, GU H M. Namerical modeling of Rayleigh wave using high-accuracy finite-difference method and wave field analysis[J].Progress In Geophysics, 2009, 24(4):1313-1319.(In Chinese)

[12] 琚長輝, 王偉, 王祥春,等. 基于起伏地形的瑞利面波有限差分數值模擬[J]. 東華理工大學學報:自然科學版, 2015, 38(3):319-324.

JU C H,WANG W,WANG X C, et al. Finite-difference numerical modeling of Rayleigh wave based on relief surface [J]. Journal of East China Institute of Technology, 2015, 38(3):319-324. (In Chinese)

[13] 劉雪明, 凡友華, 翟佳羽,等. 四種典型地層的瑞雷波“之”字型頻散曲線數值模擬研究[J]. 地球物理學報, 2009, 52(12):3042-3050.

LIU X M,FAN Y H,ZHAI J Y, et al. Numerical simulation of Rayleigh wave zigzag dispersion curves of four typical strata[J].Chinese Journal Of Geophysics., 2009, 52(12):3042-3050. (In Chinese)

StratigraphicdivisioneffectofRayleighwaveexplorationintheinclinedundulation

XI Chaoqiang, ZHANG Pingsong, LI Jianling, DING Meiqin

(Anhui University of Science and Technology,Huainan 232001,China)

In the complex surface conditions, the use of Rayleigh wave exploration, the need to consider the terrain on Rayleigh wave propagation. In this paper, the Rayleigh wave signal acquisition is designed by arranging the geophone array in different directions, and the FV frequency spectrum of the seismic record after fast Fourier transform is discussed by using the high precision staggered-grid finite difference method. By comparing the signal frequency spectrum with the theoretical multi-order dispersion curve, the variation characteristics of Rayleigh wave propagation energy under different terrain conditions are analyzed, indicating that the superposition dispersion curve is extracted. The genetic algorithm is used to calculate the ground-wave inversion and the thickness of the initial model. And then the error curves of the dispersion curves of the initial model and FK method are used to calculate the global fitting error. The results show that the higher frequency dispersion of the upper excitation is observed in the multi - channel transient Rayleigh wave exploration. Therefore, the inversion of the seismic records of the lower excitation results in a relatively small thickness error.

Rayleigh wave; tilt; stratigraphic division; joint inversion

P 631.4

A

10.3969/j.issn.1001-1749.2017.05.14

2016-09-02 改回日期： 2016-11-21

席超強(1993-)，男，碩士，主要從事瑞雷面波數值模擬研究， E-mail：cichaoqiang626@163.com。

1001-1749(2017)05-0669-08