蘇悅 楊振濤 楊博 趙亮

基于矢量波數(shù)變換法的主動源瑞雷波多模式提取方法在近地表地層結(jié)構(gòu)探測中的應(yīng)用研究

蘇悅1楊振濤2,?楊博3趙亮4

1.北京大學(xué)地球與空間科學(xué)學(xué)院, 北京 100871; 2.南方科技大學(xué)地球與空間科學(xué)系, 深圳 518055; 3.中南大學(xué)地球科學(xué)與信息物理學(xué)院, 長沙 410083; 4.東方地球物理公司物探技術(shù)研究中心, 涿州 072751; ?通信作者, E-mail: yangzt@sustech.edu.cn

基于瞬態(tài)多道面波分析法(MASW)的數(shù)據(jù)采集方式, 對矢量波數(shù)變換法(VWTM)、相移法和高分辨率拉東變換法 3 種瑞雷波頻散成像方法進(jìn)行對比分析。首先, 采用合成地震數(shù)據(jù)對比 3 種方法的成像效果, 發(fā)現(xiàn) VWTM 在成像精度和成像質(zhì)量方面都有很大的優(yōu)勢, 成像分辨率高, 與理論頻散曲線擬合準(zhǔn)確, 高階模式成像效果突出, 且抗噪能力也優(yōu)于另兩種方法。然后, 采用實際多道瞬態(tài)面波數(shù)據(jù)進(jìn)行對比分析, 發(fā)現(xiàn)與相移法和拉東變換法相比, VWTM 在瑞雷波基階模式成像精度和高階模式成像質(zhì)量方面仍具有優(yōu)勢。采用遺傳算法對比基階與高階頻散曲線反演結(jié)果, 發(fā)現(xiàn)含有高階模式的多模式頻散曲線聯(lián)合反演的多解性明顯降低, 反演結(jié)果精度更高, 也更穩(wěn)定。在實際瞬態(tài)瑞雷波勘探實踐中, VWTM 能夠高質(zhì)量地提取多模式瑞雷波頻散特征, 可為多模式面波頻散聯(lián)合反演提供可靠的基礎(chǔ)數(shù)據(jù)。研究結(jié)果表明, 將 VWTM 法與有效的反演方法相結(jié)合, 可以極大地提高瑞雷波淺地層探測能力和探測精度。

瑞雷波; 矢量波數(shù)變換法(VWTM); 多道面波分析法(MASW); 高階頻散曲線; 遺傳算法

利用瑞雷波頻散特性的勘探方法是一種重要的近地表地球物理勘探方法, 其中的瞬態(tài)多道面波分析(multichannel analysis of surface waves, MASW) 方法[1]在過去的 20 年中廣泛地應(yīng)用于工程地球物理勘探領(lǐng)域。目前, 從多道面波記錄中獲取頻散曲線的常見方法有-法[2]、-變換法[3]、相移法(shift phase)[1]和拉東變換法(Radon transform)[4?6]等。-法方法簡單, 易于實現(xiàn), 但需要滿足空間域與時間域采樣間隔相等的條件, 且對總道數(shù)要求較高, 如果波場記錄中存在壞道, 會對結(jié)果產(chǎn)生較大的影響, 使橫向分辨率大大降低, 從而降低提取頻散曲線的精度[7]。-變換法和線性拉東變換法均將數(shù)據(jù)沿一系列直線進(jìn)行疊加, 高階模式成像質(zhì)量較好, 但基階模式頻散曲線在低頻段成像效果較差, 容易出現(xiàn)假頻現(xiàn)象和端點效應(yīng)[7?8]。相移法是先將各道地震數(shù)據(jù)進(jìn)行傅里葉變換, 然后在空間域進(jìn)行積分, 在不同頻率下進(jìn)行速度掃描, 對基階模式的成像效果較好, 但對高階模式的成像分辨率不高[8?9]。

綜上所述, 現(xiàn)有利用瑞雷波頻散特性的勘探方法對低頻部分和高階模式的成像均具有局限性。隨著工程勘查對探測精度要求的不斷提高, 對面波勘探也提出更多的要求。研究發(fā)現(xiàn), 高階面波在地層結(jié)構(gòu)反演分析中起著關(guān)鍵作用[9?14]。若能在頻散能量圖中有效地提高高階模式的成像質(zhì)量, 將極大地提高面波反演的精度和穩(wěn)定性。因此, 本研究分析空間格林函數(shù)的矢量波數(shù)變換法(vector wavenumber transformation method, VWTM)(在噪聲成像中也稱為頻率?貝塞爾變換法(frequency-Bessel transform method, F-J)[15]的成像能力, 并與相移法及高分辨率線性拉東變換法進(jìn)行對比。在采集道數(shù)不變的情況下, VWTM能提高瑞雷波高階模式的成像質(zhì)量和抗噪能力; 通過 VWTM 提取瑞雷波頻散曲線進(jìn)行多模式聯(lián)合反演, 可明顯提高瑞雷波的探測精度。

1 矢量波數(shù)變換(VWTM)分析法原理

在水平半無限空間模型中, 震源為原點, 在柱坐標(biāo)中, 地表觀測點觀測到的瑞雷波可表示為

其中,表示頻率。根據(jù)水平層狀介質(zhì)中地震波傳播理論[15?16], 地表觀測點與震源間的格林函數(shù)在頻率域的計算公式可以表示為

其中,是格林核函數(shù),0是第一類零階貝塞爾函數(shù)。對 G(,)進(jìn)行矢量波數(shù)變換[9,15], 得到

其中,表示波數(shù)。當(dāng)空間為水平層狀介質(zhì)時,=||, 其矢量波數(shù)變換結(jié)果僅依賴波數(shù)矢量的模(), 與波數(shù)矢量的方向無關(guān), 此時式(4)可表示為

將式(3)代入式(5), 可得式(6)[9]:

根據(jù)貝塞爾函數(shù)的正交性:

將式(7)代入式(6), 則式(6)可以簡化為

將式(2)代入式(8), 可得

式(8)中, 格林函數(shù)的核函數(shù)值與確定面波頻散特性的久期函數(shù)值成反比[16]:

2 合成數(shù)據(jù)分析

2.1 成像對比分析

為了測試 VWTM 的成像能力, 首先采用離散波數(shù)法[17?18], 模擬垂直于地面單力點源產(chǎn)生的地震波場, 然后采用 3 種方法(VWTM、相移法和高分辨率線性拉東變換法)進(jìn)行成像, 并與廣義反射?透射系數(shù)方法[19?20]得到的理論頻散曲線進(jìn)行疊加對比。

本研究采用含低速層的 4 層水平層狀模型, 地層模型參數(shù)見表 1。合成地震數(shù)據(jù)的震源時間函數(shù)采用雷克子波近似, 中心頻率為 10Hz, 延遲為 0.2 s。震源和檢波器呈線性排列, 道間距為 1m, 最小偏移距為 10m, 數(shù)據(jù)長度為 2s, 采樣間隔為 1ms, 合成 84 道地震數(shù)據(jù)(圖 1)。分別采用 24 和 84 道地震數(shù)據(jù), 應(yīng)用 VWTM、相移法和高分辨率線性拉東變換法[5]進(jìn)行計算, 得到頻散能量圖, 再與理論頻散曲線疊加, 對比其成像效果(圖2)。

表1 水平層狀地層模型參數(shù)[9]

通過對比成像結(jié)果發(fā)現(xiàn), 當(dāng)參與處理的地震道數(shù)增加(即排列增長)時, 3 種方法的頻散能量圖在低頻區(qū)域的成像質(zhì)量都明顯提高, 表明瑞雷波勘探深度和精度會隨著臺站數(shù)量和排列長度的增加而提高。圖 2 顯示, 相移法的整體分辨率較低, 尤其是高階面波的成像分辨率明顯弱于其他兩種方法; 當(dāng)只有 24 道地震數(shù)據(jù)參與成像時效果非常差, 10Hz以下的部分幾乎完全不能識別, 隨著地震道數(shù)增加, 基階模式成像質(zhì)量明顯提高, 但高階模式成像質(zhì)量仍較差。高分辨率線性拉東變換法成像分辨率較高, 尤其在地震道數(shù)較少時, 也能得到分辨率較高的頻散圖像, 但與理論頻散曲線對比出現(xiàn)一定的偏差(圖 2); 隨著地震道數(shù)增加, 拉東變換法整體成像效果有所提高, 但在小于 5Hz 區(qū)域基階模式成像與理論頻散曲線出現(xiàn)較大偏差。與前兩種方法相比, VWTM 的成像質(zhì)量具有明顯優(yōu)勢, 尤其在 4~30Hz頻率范圍內(nèi)基階和高階模式的成像質(zhì)量和精度都明顯優(yōu)于其他兩種方法; 雖然地震道數(shù)較少時不如高分辨率線性拉東變換法的分辨率高, 但地震道數(shù)增加后分辨率明顯提高, 成像精度明顯優(yōu)于其他兩種方法。因此, VWTM 在高階模式的成像質(zhì)量和基階模式的成像精度方面具有巨大的優(yōu)勢。

2.2 抗噪能力分析

為分析成像方法的抗噪能力, 我們模擬隨機(jī)壞道, 每道加入均值為 0, 標(biāo)準(zhǔn)差為 1 的隨機(jī)干擾, 噪聲水平為原始數(shù)據(jù)中每道振幅最大值平均值的 20% (圖 3), 并隨機(jī)將各道的噪聲水平放大 1 倍、縮小50%或不變。分別采用 24 道和 84 道地震數(shù)據(jù), 應(yīng)用 VWTM、相移法和高分辨率線性拉東變換法進(jìn)行計算, 得到頻散能量圖(圖4)。

從圖 4 看出, 在處理含有噪聲的數(shù)據(jù)時, 相移法成像質(zhì)量明顯降低, 基階模式成像質(zhì)量相對穩(wěn)定, 但高頻部分(>21Hz)很難有效地成像, 且高階模式無法成像; 高分辨率拉東變換法受干擾影響較大, 尤其當(dāng)參與計算的地震道數(shù)較少時, 在高頻頻率區(qū)域(>23Hz)和低頻區(qū)域(<5Hz), 基階模式成像與理論頻散曲線存在較大的偏差, 且高階模式成像精度大大降低; VWTM 仍然是受干擾影響最小的方法, 即使在 24 道數(shù)據(jù)的情況下, 基階和高階模式仍能高質(zhì)量地成像, 且精度高, 僅在高頻部分(>25Hz)的成像受到干擾。

3 橫波速度結(jié)構(gòu)反演

3.1 遺傳算法簡介

遺傳算法是一種借鑒生物界遺傳機(jī)制和自然選擇的高效并行、隨機(jī)全局優(yōu)化搜索算法, 通過模擬自然進(jìn)化的過程來搜索最優(yōu)解, 只計算正演和適應(yīng)性函數(shù), 可以避免矩陣求逆和偏導(dǎo)數(shù)計算而可能產(chǎn)生的數(shù)值不穩(wěn)定或精度丟失等問題。該方法 1975 年由 Holland[21]提出, 現(xiàn)已廣泛推廣并成功地應(yīng)用于多個領(lǐng)域, 也是常用的瑞雷波頻散曲線反演方法。

作為一種搜索尋優(yōu)技術(shù), 遺傳算法從問題參數(shù)編碼出發(fā), 并按照一定的操作規(guī)則(如選擇、交叉和變異等)初始化種群, 然后按照一定的標(biāo)準(zhǔn)更新群體, 實現(xiàn)一代遺傳; 如此反復(fù)迭代, 逐步逼近問題的解。因為面波頻散曲線反演是一種高度非線性、多維度的復(fù)雜問題, 因此本文的反演研究采用遺傳算法的實數(shù)編碼遺傳算法來進(jìn)行問題的求解, 相對于二進(jìn)制編碼, 這種編碼更適合求解大規(guī)模復(fù)雜函數(shù)優(yōu)化問題, 且求解精度高, 算法運(yùn)行快, 結(jié)果更穩(wěn)定[22]。早期的遺傳算法中使用的交叉和變異算子是固定不變的, 本研究采用自適應(yīng)的交叉和變異算子, 以便靈活地應(yīng)對面波頻散曲線反演問題。自適應(yīng)的交叉和變異算子的概率計算公式為

3.2 模型反演

選取 84 道數(shù)據(jù), 采用 VWTM 得到的頻散能量圖(圖 2(a2)), 提取頻散曲線, 分別用不同的模式進(jìn)行反演, 并與 4 層理論模型進(jìn)行對比。

3.2.1 僅用基階模式進(jìn)行反演

從圖 6(a)可以看出, 僅采用基階模式反演得到的基階頻散曲線與實測基階頻散點在高頻部分能夠很好地吻合, 但在低頻區(qū)域(<5Hz)與理論頻散曲線和實測頻散點之間仍有一定程度的誤差。

圖 6(b)為 10 次反演的結(jié)果(每次反演迭代 50次, 目標(biāo)函數(shù)均趨于收斂), 得到 1~3 層的層厚標(biāo)準(zhǔn)差分別為 0.18, 1.43 和 2.92 m, 1~4 層的橫波速度標(biāo)準(zhǔn)差分別為 0.50, 15.01, 30.09 和 47.72 m/s?？梢钥吹? 僅采用基階模式反演, 雖然頻散曲線趨于擬合, 但多次反演的結(jié)果仍具有較大的不穩(wěn)定性。

3.2.2 多模式聯(lián)合反演

從圖 7(a)看出, 對所有實測頻散點進(jìn)行多模式聯(lián)合反演得到的理論頻散曲線與實測頻散點能夠很好地吻合。圖 7(b)為采用 10 次反演的結(jié)果(每次反演迭代 50 次, 目標(biāo)函數(shù)均趨于收斂), 得到 1~3 層的層厚標(biāo)準(zhǔn)差分別為 0.07, 0.91, 1.72 m, 1~4 層的橫波速度標(biāo)準(zhǔn)差分別為 0.18, 6.19, 10.71 和 29.28m/s。可以看到, 加入高階模式的約束后, 反演得到的頻散曲線不僅能很好地擬合, 而且反演精度和穩(wěn)定性都比僅采用基階模式反演有很大程度的提高。

4 實例研究

如前所述, 采用合成數(shù)據(jù) VWTM 可以精確地提取瑞雷波高階模式頻散曲線, 隨后采用遺傳算法多模式聯(lián)合可以高精度地反演地層結(jié)構(gòu)。為了研究VWTM 對實際數(shù)據(jù)的適用性, 我們在常州北郊長江灘涂地區(qū)進(jìn)行線性排列瞬態(tài)瑞雷波的采集。采集數(shù)據(jù)的道間距為 1m, 最小偏移距為 10m, 數(shù)據(jù)長度為 1s, 采樣頻率為 1000Hz, 24 道采集, 如圖 8 所示。分別應(yīng)用 VWTM、相移法和高分辨率線性拉東變換法進(jìn)行處理, 對比 3 種方法對實際數(shù)據(jù)的成像效果, 如圖 9 所示。從圖中可以看到, VWTM 的成像質(zhì)量明顯優(yōu)于其他兩種方法, 尤其高階模式的成像優(yōu)勢明顯, 分辨率高; 相移法的整體分辨率都很低, 低頻區(qū)域(<10Hz)的邊界效應(yīng)嚴(yán)重, 高階面波成像不理想; 高分辨率線性拉東變換法保持很高的分辨率, 但低頻區(qū)域(<5Hz)的效果很差, 誤差明顯, 且難以分辨不同模式的高階面波。

選取圖 9(a)中 VWTM 提取多模式頻散曲線的結(jié)果, 根據(jù)搜集到的地質(zhì)資料建立一個 7 層初始模型進(jìn)行反演。通過提取頻散能量圖中的能量極值點獲得頻散點, 并對頻散點進(jìn)行遺傳算法多模式反演, 結(jié)果如圖 10 所示。反演得到的頻散曲線與實測頻散點基本上吻合, 地層結(jié)構(gòu)為速度遞增正常沉積地層。該結(jié)果與實際勘察資料相符。

5 結(jié)論

本文將一種新的多道面波分析方法——矢量波數(shù)變換法(VWTM)與傳統(tǒng)相移法以及高分辨率線性拉東變換法對合成數(shù)據(jù)和實際數(shù)據(jù)進(jìn)行成像對比分析, 發(fā)現(xiàn) VWTM 在成像精度和成像質(zhì)量方面都有巨大的優(yōu)勢, 尤其是 VWTM 的高階模式成像能力和抗噪能力明顯優(yōu)于其他兩種方法; 在實際瑞雷波探測中, 采集道數(shù)不是很多的情況下, VWTM 仍能保證高階模式成像的分辨率。因此, VWTM 能為實際地層多模式聯(lián)合反演提供更為豐富的頻散信息。

在得到多模式頻散曲線的基礎(chǔ)上, 本文采用改進(jìn)的實數(shù)編碼自適應(yīng)遺傳算法對基階、高階頻散曲線進(jìn)行反演對比, 發(fā)現(xiàn)對含有高階模式的頻散信息進(jìn)行聯(lián)合反演時多解性明顯降低, 反演結(jié)果精度更高, 也更穩(wěn)定。在實際瞬態(tài)瑞雷波勘探中, VWTM能夠有效地高質(zhì)量地提取多模式瑞雷波頻散特征, 為基階與高階面波頻散聯(lián)合反演提供可靠的數(shù)據(jù)基礎(chǔ)。VWTM 結(jié)合有效的反演方法, 可以極大地提高近地表地層結(jié)構(gòu)的探測精度, 尤其在城市工程勘察中具有巨大的應(yīng)用潛力。

VWTM 是基于檢波器與震源之間的格林函數(shù)進(jìn)行頻率?波數(shù)域貝塞爾變換, 所以檢波器可以任意排列, 為了與其他方法進(jìn)行比較, 本文僅采用線性觀測進(jìn)行對比分析。本研究基于水平層狀介質(zhì)模型, 將震源近似為雷克子波, 對實際情況做了近似處理。然而, 實際地下結(jié)構(gòu)非常復(fù)雜，震源時間函數(shù)也不是簡單的雷克子波。因此, 本文的 VWTM方法仍有較大的改善空間。

[1] Park C B, Miller R D, Xia J. Multichannel analysis of surface waves. Geophysics, 1999, 64(3): 800?808

[2] Capon J. High-resolution frequency-wavenumber spec-trum analysis. Proceedings of the IEEE, 2005, 57(8): 1408?1418

[3] Mcmechan G A, Yedlin M J. Analysis of dispersive by wave field transformation. Geophyscis, 1981, 6(46): 869?874

[4] 余欽, 范沈操, 牛濱華. Radon 變換的 MATLAB 實現(xiàn). 物探化探計算技術(shù), 2000, 22(4): 345?350

[5] Luo Y, Xia J, Miller R D, et al. Rayleigh-wave mode separation by high-resolution linear Radon transform. Geophysical Journal International, 2009, 179(1): 254? 264

[6] 潘東明, 胡明順, 崔若飛, 等. 基于拉東變換的瑞雷面波頻散分析與應(yīng)用. 地球物理學(xué)報, 2010, 53 (11): 2760?2766

[7] 盧建旗. 多道面波分析方法及其應(yīng)用研究[D]. 北京: 中國地震局工程力學(xué)研究所, 2013

[8] 邵廣周, 李慶春. 聯(lián)合應(yīng)用-變換法和相移法提取面波頻散曲線. 石油地球物理勘探, 2010, 45(6): 836?840

[9] 楊振濤, 陳曉非, 潘磊, 等. 基于矢量波數(shù)變換法(VWTM)的多道 Rayleigh 波分析方法. 地球物理學(xué)報, 2019, 62(1): 298?305

[10] Zhang B X, Lu L Y, Bao G S. A study on zigzag dis-persion curves in Rayleigh wave exploration. Chinese Journal of Geophysics, 2002, 45(2): 265?276

[11] 張碧星, 魯來玉, 鮑光淑. 瑞利波勘探中“之”字形頻散曲線研究. 地球物理學(xué)報, 2002, 45(2): 263? 274

[12] 凡友華. 考慮高階模的 Rayleigh 波勘探應(yīng)用研究[D]. 北京: 北京大學(xué), 2003

[13] Pan L, Chen X, Wang J, et al. Sensitivity analiysis of dispersion curves of Rayleigh waves with fundamen-tal and higher modes. Geophyscial Journal Interna-tional, 2019, 216(2): 1276?1303

[14] Xia J, Miller R D, Park C B, et al. Inversion of high frequency surface waves with fundamental and higher modes. Journal of Applied Geophysics, 2003, 52(1): 45?57

[15] Wang J, Wu G, Chen X. Frequency-Bessel transform method for effective imaging of higher-mode Ray-leigh dispersion curves from ambient seismic noise data. Journal of Geophysical Research Solid Earth, 2019, 124(4): 3708?3723

[16] Chen X. Seismogram synthesis in multi-layered half-space partⅠ. theoretical formulations. Earthquake Research in China, 1999, 13(2): 53?78

[17] Zhang H M, Chen X F, Chang S. An efficient numeri-cal method for computing synthetic seismograms for a layered half-space with sources and receivers at close or same depths. Pure & Applied Geophysics, 2003, 160(3/4): 467?486

[18] 王建楠. 背景噪音提取高階頻散曲線的矢量波數(shù)變換方法[D]. 合肥: 中國科學(xué)技術(shù)大學(xué), 2017

[19] Chen X. A systematic and efficient method of compu-ting normal modes for multilayered half-space. Geo-physical Journal International, 1993, 115(2): 391?409

[20] 何耀鋒, 陳蔚天, 陳曉非. 利用廣義反射?透射系數(shù)方法求解含低速層水平層狀介質(zhì)模型中面波頻散曲線問題. 地球物理學(xué)報, 2006, 49(4): 1074?1081

[21] Holland J H. Adaptation in natural and artificial sys-tems. Ann Arbor, 1975, 6(2): 126?137

[22] 周永華. 實數(shù)編碼遺傳算法雜交算子組合研究[D]. 廣州: 華南理工大學(xué), 2003

[23] 毛承英. 基于改進(jìn)遺傳算法的瑞雷波頻散曲線反演[D]. 長沙: 中南大學(xué), 2010

Application Research of Active Source Rayleigh Wave Multi-Mode Extraction Method Based on Vector Wavenumber Transformation Method in Near Surface Stratigraphic Structure Detection

SU Yue1, YANG Zhentao2,?, YANG Bo3, ZHAO Liang4

1. School of Earth and Space Sciences, Peking University, Beijing 100871; 2. Department of Earth and Space Sciences, Southern University of Science and Technology, Shenzhen 518055; 3. School of Geosciences and Info-physics, Central South University, Changsha 410083; 4. Bureau of Geophysical Prospecting Inc., China National Petroleum Corporation, Zhuozhou 072751; ? Corresponding author, E-mail: yangzt@sustech.edu.cn

Based on the array of the transient multi-channel analysis of surface waves (MASW), Rayleigh wave dispersion imaging methods including vector wavenumber transformation method (VWTM), phase shift, and high-resolution linear Radon transform are used to carry out comparative analysis. First, the dispersion imaging results of synthetic data are compared to analyze the quality of the three methods. The VWTM is found superior to the other two in terms of resolution, accuracy of higher modes, fitting degree with theoretical dispersion curves and anti-noise property. In the actual Rayleigh wave detection, the VWTM can still has advantages in foundamental mode imaging accuracy and higher order modes imaging quality. Then the genetic algorithm is employed to invert the dispersion curves of fundamental mode and higher modes. The result is more accurate and stable when higher modes are combined in the inversion, it also performs much better in reducing non-uniqueness of inversion. The research shows that the VWTM can extract multi-mode characteristics of Rayleigh wave dispersion effectively and with high quality in field prospecting, which provides a reliable data foundation for joint inversion of multi-mode of Rayleigh wave dispersion curves. In conclusion, the VWTM combined with effective inversion method can remarkably improve the accuracy of near surface stratigraphic structure detection and has great potential in enhancing Rayleigh wave exploration capability.

Rayleigh wave; vector wavenumber transformation method (VWTM); multi-channel analysis of surface waves (MASW); higher-mode dispersion curves; genetic algorithm

10.13209/j.0479-8023.2020.021

國家重點研發(fā)計劃(2018YFC0603600)和國家自然科學(xué)基金(41974047)資助

2019?05?12;

2019?06?29