基于快速標量傳遞算法的瑞雷波頻散曲線反演研究

董智開 段文勝 肖承文 胡天躍,? 張獻兵,?

基于快速標量傳遞算法的瑞雷波頻散曲線反演研究

董智開1段文勝2肖承文2胡天躍1,?張獻兵1,?

1.北京大學地球與空間科學學院, 北京 100871; 2.中國石油塔里木油田公司, 庫爾勒 841000; ?通信作者, E-mail: tianyue@pku.edu.cn (胡天躍), zxb@pku.edu.cn (張獻兵)

為了提高近地表瑞雷波頻散曲線反演的效率和精度, 引入快速標量傳遞算法來計算瑞雷波頻散曲線正演理論值。通過對比加入線性約束條件前后遺傳算法(GA)與模擬退火法(SA)在反演瑞雷波頻散曲線中的表現(xiàn), 提出將計算速度快的蒙特卡洛法(MC)作為輔助手段來快速識別地層類型, 然后在 GA 和 SA 中加入線性約束條件來提高收斂速度, 并將 GA 得到的反演結(jié)果作為 SA 的初始狀態(tài), 同時適當?shù)乜s小搜索范圍, 通過聯(lián)合反演來克服 GA 的早熟問題。用上述方法計算和驗證三層地層模型、含噪聲數(shù)據(jù)以及工區(qū)實際模型, 結(jié)果表明該方法高效、準確、穩(wěn)定性強, 有很強的全局尋優(yōu)能力, 并具有一定的抗噪能力。

淺地表; 瑞雷波; 頻散曲線; 快速標量傳遞法; 非線性反演

瑞雷面波具有能量強、振幅衰減弱、高信噪比及在層狀介質(zhì)中有頻散特性等特點[1]。近年來, 隨著高精度地震數(shù)據(jù)采集技術(shù)的發(fā)展, 因瑞雷波勘探在地球物理探測和工程勘探中的優(yōu)越性而被用于獲取淺地表、巖石圈或地幔的橫波速度結(jié)構(gòu)信 息[2?6]。瑞雷波勘探主要包括野外數(shù)據(jù)采集[7]、頻散曲線提取[8?10]和頻散曲線反演[11?14]3 個步驟。其中, 瑞雷波頻散曲線的提取是獲取觀測數(shù)據(jù)的重要一步, 在變換域發(fā)展出多種方法, 主要有 F-K 變換法[15?17]、相移法[18]、拉東變換法[8]、S 變換[1]、小波變換[19]以及在多模態(tài)面波頻散曲線提取中效果顯著的矢量波數(shù)變換法(vector wavenumber transform method, VWTM)[20]。不同頻率的一系列面波可以有相同的相速度, 相同頻率下相速度最小者為基階面波, 其次為一階、二階等高階面波[21]。在面波勘探研究領(lǐng)域, 無論是基階面波反演, 還是多模態(tài)高階面波的反演, 對于高效的頻散曲線正演計算方法的探究以及對反演優(yōu)化算法的選擇和改進都具有重要的意義。各向同性水平層狀均勻介質(zhì)的瑞雷波頻散曲線正演方法主要有 Haskell 算法[22]、Knopoff-Schwab 算法[23]、矩陣與 Abo-Zena 算法[24]、快速標量傳遞算法與快速矢量傳遞算法[25?27]以及廣義反射?透射系數(shù)法[26]等, 其中快速標量傳遞算法通過將傳遞矩陣中的元素表示成沒有量綱的實數(shù), 可以提高計算速度, 廣義反射?透射系數(shù)法則在高頻厚層介質(zhì)情況下的計算結(jié)果更加準確。

反演方法分為局部線性化反演方法和非線性反演方法[28], 局部線性化反演方法包括最小二乘法[29]、自適應(yīng)阻尼最小二乘法[30]和 Occam 算法[31]。由于此類方法采用泰勒級數(shù)展開, 近似地將瑞雷波方程線性化, 并用擾動法計算雅克比矩陣, 因此難免在計算過程中帶來誤差。瑞雷波頻散曲線反演具有高度非線性、多參數(shù)和多極值的特點[32], 要求線性化反演方法有較好的初始模型。傳統(tǒng)的非線性反演方法包括蒙特卡洛法(Monte Carlo, MC)[33?34]、遺傳算法(genetic algorithm, GA)[35?37]和模擬退火法(simulated annealing, SA)[38?40], 具有全局搜索能力, 但收斂速度慢、精度較低, 容易出現(xiàn)早熟收斂的現(xiàn)象。此外, 模式識別算法[41]、粒子群算法[42]、差分化算法[43]、螢火蟲和蝙蝠智能算法[44]及相關(guān)的改進算法[45]也應(yīng)用到面波頻散曲線的反演中。

針對傳統(tǒng)反演方法中的問題, 解決方法分為 3個方向: 1)對瑞雷波頻散曲線反演問題中目標函數(shù)的改進[46]; 2)對算法的優(yōu)化(如快速退火算法[47?48]、杰出個體保護[49]以及自適應(yīng)交叉變異運算[50])和在遺傳算法的基礎(chǔ)上加入線性約束[51?52]都可以在一定程度上解決早熟收斂的問題, 加快收斂速度, 但方法的單一性在一定程度上帶來結(jié)果的不準確性; 3)采用聯(lián)合反演方法, 充分利用各個方法的特點, 取長補短(比如 SA 與 GA 相結(jié)合的 SAGA 聯(lián)合反演方法[53?54]、線性反演方法與非線性反演方法相結(jié)合的自適應(yīng) GA 與阻尼最小二乘法(damped least squares, DLS)聯(lián)合反演算法[55]以及瑞雷波與垂直電測深數(shù)據(jù)聯(lián)合反演[56]), 但多種反演方法的同時應(yīng)用難免犧牲計算速度, 影響數(shù)據(jù)處理的效率。

為了提高反演效率和精度, 本文引入計算效率較高的快速標量傳遞算法來正演理論頻散曲線, 考察 MC 識別地層類型的可能性, 為后續(xù)在啟發(fā)式蒙特卡洛法中加入線性約束條件提供依據(jù), 并基于三層遞增型地層模型, 分析線性約束條件在 SA 和GA 中的作用。利用先驗信息[57]或者計算速度快的MC 多次反演輔助識別出地層類型后, 在 GA 和 SA中加入相應(yīng)的線性約束條件, 提高收斂速度, 之后將 GA 得到的反演結(jié)果作為 SA 的初始狀態(tài), 以線性約束下 GA 的反演結(jié)果為參考, 適當縮小搜索范圍, 進行聯(lián)合反演, 克服 GA 的早熟和 SA 搜索效率低的缺陷, 壓制多解性, 達到同時提高正、反演計算效率和結(jié)果準確性的目的, 并通過對三層地層模型、含噪聲數(shù)據(jù)和工區(qū)實際模型的計算, 證實方法的有效性。

1 瑞雷波頻散曲線正演算法

1.1 正演模型

近地表一般為風化層、砂礫巖、沙漠、戈壁和黃土等, 其共同特點是橫波速度較慢, 地震記錄中面波十分發(fā)育。本研究基于前人的研究結(jié)果[31?33], 并根據(jù)近地表實際情況, 設(shè)計較為常見的遞增型、含高速硬夾層以及含低速軟夾層 3 種三層水平均勻各向同性地層模型(表 1), 設(shè)定地層密度均為 1.9 g/cm3, 第一層和第二層的厚度均為 5 m, 第三層為無限厚度介質(zhì)。

1.2 快速標量傳遞法正演瑞雷波頻散曲線

快速標量傳遞算法[25?27]是在 Abo-Zena 算法[24]和 Menke 算法[58]的基礎(chǔ)上改進的一種快速正演瑞雷波頻散曲線的算法, 其基本思想是重新定義Menke算法中的矢量, 將原本形式復(fù)雜、帶有復(fù)數(shù)運算及有量綱運算的傳遞矩陣分解為 3 個形式簡單的五階矩陣的乘積形式, 各矩陣中有些元素為零, 某些元素在矩陣乘法運算中可以消去, 從而改進為標量運算和無量綱運算, 可以大大地提高計算速度[59], 具體方法如下。

表1 理論模型參數(shù)

重新定義Menke65B9法的矢量:

其中,

基于以上快速標量傳遞算法, 采用二分法計算正演理論頻散曲線: 1)固定頻率的數(shù)值, 對面波相速度給定某一掃描區(qū)間和搜索步長, 使面波相速度從掃描區(qū)間自下而上取值; 2)搜索有根區(qū)間, 利用快速標量傳遞算法計算相鄰相速度值對應(yīng)的頻散方程值, 只要找到一個有根區(qū)間, 就調(diào)用二分法的求根子函數(shù)求得零點值, 其對應(yīng)的相速度值就是該頻率下的瑞雷波頻散曲線值; 3)最后, 改變頻率的數(shù)值, 重復(fù)第 2 步, 得到不同頻率數(shù)值下對應(yīng)的使得頻散方程成立的面波相速度值。

圖 1 為采用快速標量傳遞算法正演的表 1 中 3種地層的理論頻散曲線, 在 100Hz 以內(nèi)可見 3 種模式的瑞雷波頻散曲線, 且不同類型地層模型的頻散曲線形態(tài)有較明顯的區(qū)別。正演計算是反演工作的基礎(chǔ), 由于高階模式瑞雷波的反演方法與基階模式的瑞雷波反演方法基本上相同, 因此本文只反演基階模式的瑞雷波頻散曲線(灰色曲線)。

2 非線性反演算法

瑞雷波頻散方程的非線性隱式形式可以用式(8)[55]表示:

常規(guī)非線性反演方法包括蒙特卡洛法(MC)、遺傳算法(GA)和模擬退火法(SA)等, 其中 MC 是所有非線性反演算法的鼻祖, GA 和 SA 則是一種啟發(fā)式的蒙特卡洛法。瑞雷波頻散曲線反演的目標函數(shù)定義為觀測數(shù)據(jù)與相速度值的均方差(RMS), 即

2.1 MC 輔助識別地層類型

傳統(tǒng)的蒙特卡洛法又稱嘗試法或窮舉法, 該方法根據(jù)給出的參數(shù)約束范圍, 在范圍內(nèi)隨機產(chǎn)生一個初始模型0, 通過正演計算目標函數(shù)值, 對初始模型0, 加入一個隨機擾動得到1, 保證1的模型空間在給出的約束范圍之內(nèi); 再次通過正演計算目標函數(shù)值, 如果1對應(yīng)的目標函數(shù)值小于0對應(yīng)的目標函數(shù)值, 則新模型被接收, 否則被排斥, 進行下一輪的迭代, 直至達到搜索次數(shù)上限或目標函數(shù)值小于預(yù)設(shè)值。即迭代中止條件為

由于沒有其他判定法則和復(fù)雜的搜索方式, 因此計算機產(chǎn)生隨機數(shù)的速度很快, 從而使 MC 的計算效率很高, 使用起來更加便捷。本文測試 MC 在寬泛的搜索范圍下識別地層類型的能力, 設(shè)定厚度的搜索范圍為 4.00~6.25m, 遞增型和含低速軟夾層地層橫波速度的搜索范圍為 200~500m/s, 含高速硬夾層地層橫波速度的搜索范圍為 200~600m/s。

從 MC 的收斂過程、反演橫波速度剖面和頻散曲線的擬合情況(圖 2)可以看出, 遞增型和含高速軟夾層地層模型在 510 次左右的時候停止更新, 最終頻散曲線擬合的均方差分別為 3.3781 和 5.3584, 含低速軟夾層地層模型的最終擬合方差為 1.9972, 反演結(jié)果總體上能夠反映地層模型的速度變化特征, 識別出地層類型。在先驗信息不充分的情況下, 可以利用 MC 多次反演, 得到具有統(tǒng)計意義的地層類型判定結(jié)果, 以此作為高效識別地層類型的輔助方法, 為后續(xù)在啟發(fā)式蒙特卡洛法中加入線性約束條件提供依據(jù)。

2.2 線性約束在 SA 與 GA 中的作用

完成地層類型識別后, 對橫波速度加入相應(yīng)的線性約束條件。以遞增型地層模型為例, 線性約束條件應(yīng)為

其中,V為第層的橫波速度。下面比較在 SA 與GA 中加入線性約束條件對收斂過程及反演結(jié)果的影響。

2.2.1 SA 的基本原理及線性約束條件的加入

SA 是一種啟發(fā)式蒙特卡洛算法, 源于統(tǒng)計熱力學。它模擬物質(zhì)退火時內(nèi)部結(jié)構(gòu)的變化過程: 物質(zhì)先在高溫下融化, 然后逐漸冷卻?；舅枷胧? 生成一系列參數(shù)向量來模擬粒子的熱運動, 通過緩慢地減小一個模擬溫度的控制參數(shù), 使模擬的熱系統(tǒng)最終冷卻結(jié)晶, 達到系統(tǒng)能量最小值。

利用 SA 反演瑞雷波頻散曲線的過程分為以下步驟。

1)在給定的約束范圍內(nèi), 隨機產(chǎn)生一個初始模型0, 或事先給定一個初始模型。設(shè)定系統(tǒng)的初始溫度0、溫度變化率、最大迭代次數(shù)max以及擬合誤差閥值等參數(shù), 初次迭代開始時取=1。

3)根據(jù) Metropolis 接受準則, 采用下式進行判斷:

其中,為狀態(tài)躍遷的概率,b為 Boltzmann 常數(shù)。根據(jù)式(12)計算得到的概率值進行判定, 然后進行下一步的計算。如果=1 或>, 則接受新模型, 否則保持原來的模型。

4)采用雙曲下降型退火機制T=0a進行迭代, 每迭代一次,值加 1。迭代中止條件同 MC, 重復(fù)步驟 2 和步驟 3。

計算過程中, 設(shè)置初始溫度為 5°, 溫度變化率為 0.96。因為已判定為遞增型地層, 為了測試 SA的全局尋優(yōu)能力, 給定一個遠離全局最優(yōu)的遞增型地層初始模型: 橫波速度從上至下依次為 300, 310和 420m/s, 初始厚度值從上至下依次為 4.5 和 5.5m, 第三層為半無限空間的介質(zhì)。由于在實際工作中, 通常會結(jié)合先驗信息給出合理的約束范圍, 在利用 MC 識別出地層類型后, 結(jié)合測井信息, 搜索范圍可按類給出, 因此本試驗中將所有待反演參數(shù)的約束范圍設(shè)定為理論值的 0.8~1.25 倍。加入線性約束條件指在每次退火之前, 即在進行 Metropolis接受準則判斷之前, 先對加入隨機擾動的新模型進行篩選, 符合線性約束條件的新模型進入下一步的操作。

在相同搜索次數(shù)下, 未加入線性約束條件的SA 最終擬合均方差為 9.8018, 加入線性約束條件的 SA 最終擬合均方差為 3.4640。從線性約束下反演所得參數(shù)值與理論值對比(表 2)來看, SA 的結(jié)果整體上好于常規(guī) SA 方法, 收斂過程(圖 3)也表明, 加入線性約束的 SA 收斂過程更加連續(xù)穩(wěn)定, 在 200次搜索后依然保持較快的向全局最優(yōu)解收斂速度, 能夠獲得更加精確的結(jié)果。反演結(jié)果也說明, SA 雖然具有全局搜索能力, 但當初始狀態(tài)不在全局最優(yōu)解附近時, 收斂至全局最優(yōu)解將變得非常困難, 顯示其搜索效率的低下。除第一層的反演速度與理論值基本上吻合外, 其他兩層反演的速度和厚度都相差較遠, 目標函數(shù)收斂于距初始模型較近的局部最優(yōu)解附近, 意味著跳出至距初始模型值更遠的全局最優(yōu)解將花費大量的時間。

2.2.2 GA 的基本原理及線性約束條件的加入

GA 是一種模擬生物進化的自然選擇和遺傳過程, 依據(jù)適者生存原則建立的仿生式蒙特卡洛法, 利用 GA 反演瑞雷波頻散曲線的過程包括編碼、初始種群的生成、選擇、交叉和變異 5 個部分。

1)編碼。采用二進制編碼, 將反演參數(shù)組轉(zhuǎn)換生成 0 和 1 字符, 每一個參數(shù)對應(yīng)一串預(yù)設(shè)固定長度的二進制編碼, 一個參數(shù)組編碼連接起來形成一條染色體。

2)初始種群的生成。按照二進制編碼格式, 通過隨機函數(shù)產(chǎn)生具有一定群體規(guī)模的初始種群, 其中每一條染色體為一個個體。利用正演方法, 計算每個個體對應(yīng)的目標函數(shù)值(m)(m指第個個體,=1, 2, …,,為群體規(guī)模), 并將目標函數(shù)轉(zhuǎn)換為適應(yīng)度函數(shù):

其中,(m)為個體m的適應(yīng)度,max是一個大于目標函數(shù)最大值的較大整數(shù)。

3)選擇。將選擇概率(m)作為標準, 選出優(yōu)良個體作為父本繁殖下一代。(m)的計算公式為

然后, 利用計算得到的選擇概率, 采用經(jīng)典輪盤賭的方式, 循環(huán)選擇父本, 生成新的群體。在 0 與 1之間生成隨機數(shù), 若(1)+(2)+…+(m–1)<<(1)+(2)+…+(m), 則個體m被保留。個體的適應(yīng)值越大, 目標函數(shù)值越小, 個體被選中的幾率也越大。通過這樣選擇, 優(yōu)良個體被保留, 并遺傳至下一代, 劣質(zhì)個體被淘汰。

表2 SA 搜索范圍及加入線性約束條件前后的反演結(jié)果

4)交叉。對選擇產(chǎn)生的新群體進行隨機性的兩兩配對, 對部分配對的個體, 利用一定的概率, 在參數(shù)編碼的低位進行基因互換(染色體對應(yīng)編碼互換), 目的是在優(yōu)良個體周圍進行隨機優(yōu)化。

5)變異。為了維持種群的多樣性, 讓群體進化到后期依然保持獲得優(yōu)良基因的能力, 利用一定的變異概率, 對個體參數(shù)編碼進行反轉(zhuǎn)。交叉和變異過程見圖 4, 其中虛線部分為交叉后從其他個體獲得的序列。

6)循環(huán)選擇、交叉、變異過程, 直至達到中止條件。

根據(jù)文獻[60], 需要對參數(shù)給予適當?shù)脑O(shè)置。計算過程中, 設(shè)置二進制編碼長度為 7, 群體規(guī)模為 32, 交叉概率為 0.6, 變異概率為 0.02, 搜索次數(shù)為 200, 搜索范圍和中止條件同 SA。在 GA 中加入線性約束條件的操作方法: 在每一次迭代過程中, 完成選擇、交叉、變異的過程之后, 對群體中的優(yōu)良個體進行篩選, 滿足線性約束條件的個體進入下一代。

經(jīng)典 GA 最終擬合均方差為 1.3312, 加入線性約束條件的最終擬合均方差為 0.5492。整體上, GA的反演結(jié)果已經(jīng)較為接近理論值(表 3), 加入線性約束后的 GA 得到更加精準的結(jié)果, 尤其第三層介質(zhì)橫波速度的反演結(jié)果明顯優(yōu)于經(jīng)典 GA, 但第二層介質(zhì)的橫波速度反演結(jié)果仍不太理想。從收斂過程(圖 5)來看, 加入線性約束的 GA 有明顯的改善, 其收斂過程更加穩(wěn)定連續(xù), 在迭代搜索 60 次后進入早熟狀態(tài), 而經(jīng)典 GA 僅在反演初期有較好的收斂速度, 后續(xù)收斂速度明顯變慢, 提前進入早熟收斂狀態(tài)。

3 線性約束聯(lián)合反演

第 2 節(jié)的計算結(jié)果表明, 線性約束條件可以使SA 和 GA 具有更快的收斂速度和更精確的收斂結(jié)果。GA 搜索效率很高, 但后期容易因群體多樣性的降低而進入早熟狀態(tài)。SA 搜索面廣, 既可以向目標函數(shù)增大的方向進行搜索, 也可以向目標函數(shù)減小的方向進行尋找, 但搜索效率低下。

在針對瑞雷波頻散曲線的前期反演中, 加入線性約束的 GA 反演結(jié)果已經(jīng)收斂于靠近全局極小值的某一個局部最小值, 因此可以用加入線性約束條件的 GA 反演結(jié)果作為 SA 的初始狀態(tài)進行聯(lián)合反演。其中, SA 也加入線性約束條件, 相當于讓 GA幫助 SA 找到一種能量較低的平衡態(tài), 再利用 SA 在這種狀態(tài)下跳出局部最小值的突出能力, 達到收斂至全局最優(yōu)的目的。

表3 GA搜索范圍及加入線性約束條件前后的反演結(jié)果

該方法的流程如圖 6 所示: 先利用 MC 或者先驗信息, 快速地識別地層類型后, 在 SA 和 GA 中加入相應(yīng)的線性約束條件, 給予合理的約束范圍, 先進行 GA 搜索, 進入早熟狀態(tài)后, 用 GA 反演的結(jié)果作為 SA 的初始狀態(tài)。由于 GA 結(jié)果已經(jīng)靠近全局極小值, 因此可以以 GA 反演結(jié)果為參考, 適當?shù)乜s小搜索范圍, 進行 SA 搜索, 直至達到中止條件。

3.1 無噪聲理論模型計算及不同方法反演誤差對比

在加入線性約束條件的 GA 反演結(jié)果基礎(chǔ)上, 將搜索范圍適當?shù)乜s小, 從理論值的 0.8~1.25 倍調(diào)整為 GA 反演結(jié)果的 0.9~1.1 倍, 以 GA 反演結(jié)果作為 SA 的初至狀態(tài), 設(shè)置初始溫度為 5°, 溫度變化率為 0.96, 在 GA 的反演結(jié)果附近進行SA搜索。

經(jīng)過聯(lián)合反演, 頻散曲線的最終擬合誤差為0.0929, 計算得到第一層至第三層的橫波速度依次為 249.95, 250.92 和 448.85 m/s, 第一層和第二層的厚度分別為 5.01 和 5.00 m。收斂過程(圖 7)顯示, 在GA 進入早熟狀態(tài)的情況下, SA 能夠以此為起點, 在反演初期迅速跳出局部最小值的限制, 最終收斂至靠近全局最優(yōu)解的位置, 頻散曲線擬合度很高。反演得到的數(shù)值結(jié)果和橫波速度剖面圖(圖 8)表明, 聯(lián)合反演得到的橫波速度和地層厚度與理論值幾乎完全吻合, 證明該方法對于瑞雷波頻散曲線非線性反演問題具有很強的全局尋優(yōu)能力。

圖 9 對比線性約束下的聯(lián)合反演、未加入線性約束的 SA、加入線性約束的 SA、未入線性約束的GA 和加入線性約束的 GA 擬合的頻散曲線, 可以看出, SA 方法受初始狀態(tài)影響, 擬合效果不太理想, 而 GA 和聯(lián)合反演方法擬合的頻散曲線均能與觀測數(shù)據(jù)較好地吻合, 從側(cè)面反映出瑞雷波反演問題的多解性。

定義各參數(shù)反演結(jié)果的相對誤差:

以遞增型地層模型為例, 加入線性約束的反演方法能得到比經(jīng)典方法更貼近理論值的速度剖面, 且相對誤差(表 4)減小, 說明線性約束可在一定程度上壓制多解性, 提高結(jié)果的準確性。GA 對面波頻散曲線有較好的反演效果, 最小誤差為 0.29%, 但最大誤差仍然達到 4.53%。因此, 單純使用 GA進行反演, 結(jié)果中可能出現(xiàn)個別參數(shù)反演誤差過大的現(xiàn)象, 影響整體精度, 不能滿足實際工作的要求, 無法作為最終結(jié)果使用, 但可以作為基礎(chǔ)供二次反演參考。聯(lián)合反演結(jié)果的最大誤差僅 0.2%, 最小誤差約為 0, 且速度剖面圖(圖 10)中顯示出極高的吻合度, 進一步證明該方法的有效性。

采用加入線性約束的聯(lián)合反演算法, 計算另外兩類地層模型的線性約束條件。含高速硬夾層地層模型的線性約束條件為

含低速軟夾層地層模型的線性約束條件為

利用 MC 判定出地層類型后, 將約束范圍設(shè)定為適用于對應(yīng)地層類型的數(shù)值進行反演。計算得到含高速硬夾層模型和含低速軟夾層模型頻散曲線的最終擬合誤差分別為 0.3236 和 0.0685; 含高速硬夾層地層第一層至第三層的橫波速度依次為 249.68, 552.01 和 296.91m/s, 第一層和第二層的厚度分別為 4.98 和 4.97m; 含低速軟夾層地層第一層至第三層的橫波速度依次為 350.08, 315.12 和 450.2m/s, 第一層和第二層的厚度分別為 4.98 和 5.04m。圖11 顯示, 頻散曲線和反演得到的速度剖面幾乎都與理論值完全吻合, 各參數(shù)的相對誤差在 1%以內(nèi), 表明該方法對不同類型地層的模型均有較好的反演效果。

表4 不同方法反演結(jié)果的相對誤差(%)

3.2 含噪聲理論模型計算

在瑞雷波勘探的反演過程中, 將頻散曲線提取的瑞雷波相速度值作為觀測數(shù)據(jù), 其主要噪聲來源于提取方法的局限性和地震記錄的信噪比。選擇直接在理論頻散曲線(觀測數(shù)據(jù))上加噪聲, 疊加了提取方法和原始地震記錄帶來的噪聲。為了檢驗本文方法的抗噪能力, 在 3 種地層模型的瑞雷波頻散曲線觀測數(shù)據(jù)中加入 10%的隨機噪聲, 搜索范圍與對應(yīng)的無噪聲模型相同。設(shè)置 MC 搜索次數(shù)為 2000, SA 初始溫度為 5°, 溫度變化率為 0.96, GA 編碼長度為 7, 群體規(guī)模為 32, 交叉概率為 0.6, 變異概率為 0.02, GA 搜索次數(shù)為 100。早熟收斂后進入 SA操作, 設(shè)置 SA 搜索次數(shù)為 2000, 用線性約束 GA 與SA聯(lián)合反演方法進行計算。

采用上述加線性約束聯(lián)合反演方法, 含噪聲遞增型、含噪聲含高速硬夾層和含噪聲低速軟夾層地層模型頻散曲線的最終擬合均方差分別為 3.6935, 3.6384 和 3.1181。計算結(jié)果表明, 對含噪聲數(shù)據(jù), MC 均能很好地識別地層類型。含噪聲遞增型地層模型反演結(jié)果的最大誤差為 2.60%(第二層厚度), 最小誤差為 0.22%(第一層橫波速度); 含噪聲含高速硬夾層地層模型反演結(jié)果的最大誤差為 3.38% (第三層橫波速度), 最小誤差為 0.33%(第一層橫波速度); 含噪聲含低速軟夾層地層模型反演結(jié)果最大誤差為 7.10%(第一層厚度), 最小誤差為 0.21%(第一層橫波速度)。結(jié)合速度剖面、頻散曲線的吻合程度(圖 12), 說明橫波速度參數(shù)反演結(jié)果整體上好于厚度參數(shù)反演結(jié)果, 淺層的反演結(jié)果略好于深層的反演結(jié)果, 遞增型地層反演效果最佳, 雖然與無噪聲數(shù)據(jù)的反演結(jié)果相比精度均有所下降, 但速度曲線吻合得很好, 表明反演結(jié)果與理論值相差很小, 精度可以滿足實際工作要求, 可見該方法具有一定的抗噪能力。

4 實際模型計算實例

某工區(qū)某井的調(diào)查結(jié)果如下: 淺地表低速帶厚度為80 m, 從上至下依次為松散黃土層、黃土層和黃土夾沙層, 對應(yīng)的縱波速度分別為 419, 765 和976m/s, 橫波速度約為 201, 368 和 469m/s, 第一層和第二層的厚度分別為 8.1 和 36.4m, 各層的密度約為 1.6, 1.79 和 2.22g/cm3。

將該地區(qū)的實際模型近似為一維層狀模型, 采用高精度有限差分法[61?62], 對該模型進行全波場正演模擬。雷克子波主頻為 20Hz, 子波寬度為 0.1s, 差分網(wǎng)格大小d=d=2m, 網(wǎng)格數(shù)為 202×202, 吸收層厚度為 100m, 采樣間隔為 0.2ms。上表面設(shè)置為自由邊界條件, 設(shè)置兩側(cè)為 PML (perfect match layer)吸收邊界條件, 在自由表面進行激發(fā)和接收,接收時間為 1.2s。用最小偏移距 30m, 道間距 1m抽取 96 道垂直 V 分量記錄。圖 13 顯示, 模擬地震記錄中面波能量最強, 直達波次之, 同相軸清晰可見。與體波勘探不同, 面波記錄在面波勘探中不再作為干擾信號被去除, 而是作為有效信息被保留。

圖13 垂直 V 分量記錄

獲取垂直V分量地震記錄后, 選取適當?shù)念l率和相速度搜索范圍, 采用 F-K 法計算其頻率相速度譜。歸一化后的 F-K 域能量譜和 F-V 域能量譜見圖 14, 可見求得的頻率波數(shù)譜和頻率相速度譜均連續(xù)且清晰, 以基階頻散曲線能量為最強, 高階頻散曲線能量相對較弱。將每一頻率值對應(yīng)能量最強的相速度值連成曲線, 獲得頻散曲線觀測數(shù)據(jù)。

設(shè)置第一層至第三層的橫波速度搜索范圍分別為161~251, 294~460和375~586 m/s, 第一層和第二層的厚度搜索范圍分別為 6~10 和 29~46m。利用上述聯(lián)合反演方法進行反演, 瑞雷波頻散曲線最終擬合均方差為 0.1479。計算得到第一層至第三層的橫波速度依次為 200.92, 364.96 和 455.84m/s, 第一層和第二層的厚度分別為 8.07 和 34.92m。反演結(jié)果的最小相對誤差為 0.04% (第一層橫波速度), 最大相對誤差為 4.07% (第二層厚度), 其中第一層的橫波速度和厚度以及第二層的橫波速度反演結(jié)果相對誤差均在 1%以內(nèi)。頻散曲線(圖 15)的吻合度非常高, 數(shù)值計算結(jié)果和速度剖面圖(圖 16)都表明, 該方法對較淺層參數(shù)的反演結(jié)果非常好, 但由于面波勘探問題的多解性較強, 較深層的反演結(jié)果略差于淺層。雖然無法完全避免多解性的出現(xiàn), 但整體上可以較好地還原該地區(qū)的近地表橫波速度結(jié)構(gòu), 進一步證明該方法的有效性和實際應(yīng)用價值。

5 結(jié)論

本文構(gòu)建常見的近地表三層地層模型, 引入快速標量傳遞法進行頻散曲線的正演計算, 考察 MC在寬泛約束條件下識別地層類型的能力。基于三層遞增型地層模型, 對比加入線性約束條件前后 SA與 GA 瑞雷波頻散曲線的反演效果。識別出地層類型之后, 在 SA 和 GA 中加入線性約束條件, 將 GA的反演結(jié)果作為 SA 搜索的初始狀態(tài), 進行聯(lián)合反演, 并通過對三層地層模型、含噪聲模型和工區(qū)實際模型的正演模擬與反演計算, 得到以下結(jié)論。

1)在先驗信息不足的情況下, 可以利用 MC 計算速度快的特點, 在寬泛約束條件下, 輔助識別地層類型, 為后續(xù)加入線性約束條件提供依據(jù)。

2)線性約束條件的加入可以使得 SA 和 GA 的收斂過程更加穩(wěn)定和快速, 并在一定程度上提高反演的精度, 使 GA 的早熟收斂滯后發(fā)生。

3)雖然 SA 具有較強的全局搜索能力, 但搜索效率低下, 如果給定的初始狀態(tài)較差, SA 收斂至全局最優(yōu)解的速度將十分緩慢。GA 的搜索效率很高, 對瑞雷波頻散曲線的反演能夠取得較好的效果, 但早熟收斂的情況依然存在。

4)在快速標量傳遞法正演瑞雷波頻散曲線的基礎(chǔ)上, 采用加入線性約束條件的 GA 與 SA 進行聯(lián)合反演, 使得反演速度從正演和反演兩方面得到提升, 并且計算結(jié)果顯示該方法既能克服 GA 的早熟收斂, 也能利用 SA 搜索面廣的特點而跳出局部極小值, 壓制多解性, 具有很強的全局尋優(yōu)能力, 能夠得到與理論值相對誤差較小的解, 并具有一定的抗噪能力。

5)本文提出的方法屬于概率性算法的一種, 通常情況下均能收斂到理論值附近, 但在多次計算過程中會有一定的幾率出現(xiàn)異常解。如有此類情況發(fā)生, 需要進行辨別和剔除。

[1] 馬見青, 李慶春, 樊金生, 等. 基于 S 變換的瑞利面波頻散分析. 地球科學與環(huán)境學報, 2010, 32(3): 319?323

[2] 胡家富, 段永康, 胡毅力, 等. 利用 Rayleigh 波反演淺土層的剪切波速度結(jié)構(gòu). 地球物理學報, 1999, 42(3): 393?400

[3] 黃忠賢, 李紅誼, 胥頤. 中國西部及鄰區(qū)巖石圈 S波速度結(jié)構(gòu)面波層析成像. 地球物理學報, 2014, 57(12): 3994?4004

[4] Strobbia C L, Emam A E, Al-Genai J, et al. Rayleigh wave inversion for the near-surface characterization of shallow targets in a heavy oil field in Kuwait. First Break, 2010, 28(5): 103?109

[5] Xia Jianghai. Estimation of near-surface shear-wave velocities and quality factors using multichannel ana-lysis of surface-wave methods. Journal of Applied Geophysics, 2014, 103: 140?151

[6] Socco L V, Foti S, Boiero D. Surface-wave analysis for building near-surface velocity models — estab-lished approaches and new perspectives. Geophysics, 2010, 75(5): 75A83?75A102

[7] Zhang S X, Chan L S, Xia J. The selection of field acquisition parameters for dispersion images from multichannel surface wave data. Pure and Applied Geophysics, 2004, 161(1): 185?201

[8] Luo Y, Xia J, Miller R D, et al. Rayleigh-wave dis-persive energy imaging using a high-resolution linear radon transform. Pure and Applied Geophysics, 2008, 165(5): 903?922

[9] 王季. 基于 F-K 變換的井下多道瑞利波頻散曲線提取. 煤田地質(zhì)與勘探, 2012, 40(2): 75?77

[10] 金聰, 楊文海, 羅登貴, 等. 面波頻散曲線提取方法對比分析. 地球物理學進展, 2016, 31(6): 2735?2742

[11] Lu L, Wang C, Zhang B. Inversion of multimode Rayleigh waves in the presence of a low-velocity layer: numerical and laboratory study. Geophysical Journal of the Royal Astronomical Society, 2007, 168(3): 1235?1246

[12] 楊雅慧, 曲杉, 陳裕華, 等. 基于面波頻散曲線反演方法綜論// 2014年中國地球科學聯(lián)合學術(shù)年會——專題 14: 地下介質(zhì)結(jié)構(gòu)及其變化的地震面波、背景噪聲及尾波研究論文集. 北京: 中國地球物理學會, 2014: 762

[13] 羅銀河, 夏江海, 劉江平, 等. 基階與高階瑞利波聯(lián)合反演研究. 地球物理學報, 2008, 51(1): 242?249

[14] Zhang S, Chan L. Possible effects of misidentified mode number on Rayleigh wave inversion. Journal of Applied Geophysics, 2003, 53(1): 17?29

[15] 李欣欣, 李慶春. 利用改進的 F-K 變換法提取瑞雷波的頻散曲線. 地球物理學進展, 2017, 32(1): 197?203

[16] Shen H, Chen C, Yan Y, et al. Multiple-transient surface wave phase velocity analysis in expanded f-k domain and its application. Journal of Seismic Exp-loration, 2016, 25: 299?319

[17] Shen C, Wang A, Wang L, et al. Resolution equiva-lence of dispersion-imaging methods for noise-free high-frequency surface-wave data. Journal of Applied Geophysics, 2015, 122: 167?171

[18] Zhang S X, Chan L S. Possible effects of misi-dentified mode number on Rayleigh wave inversion. Journal of Applied Geophysics, 2003, 53(1): 17?29

[19] 房曰榮, 沈斐敏. 小波分析應(yīng)用于瑞雷波頻散曲線提取研究. 地震工程與工程振動, 2014, 34(5): 21?27

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

[21] Gao L, Xia J, Pan Y, et al. Reason and condition for mode kissing in MASW method. Pure and Applied Geophysics, 2016, 173(5): 1627?1638

[22] Haskell N A. The dispersion of surface waves on multilayered media. Bulletin of the seismological So-ciety of America, 1953, 43(1): 17?34

[23] Schwab F A, Knopoff L. Fast surface wave and free mode computations. Methods in Computational Phy-sics Advances in Research & Applications, 1972, 11: 87?180

[24] Anas A Z. Dispersion function computations for un-limited frequency values. Geophysical Journal Inter-national, 1979, 58(1): 91?105

[25] 凡友華, 肖柏勛, 劉家琦. 層狀介質(zhì)中軸對稱柱面瑞利面波頻散函數(shù)的計算. 地震工程與工程振動, 2001, 21(3): 1?5

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

[27] 凡友華, 劉家琦, 肖柏勛. 計算瑞利波頻散曲線的快速矢量傳遞算法. 湖南大學學報(自然科學版), 2002, 29(5): 25?30

[28] Xia J, Miller R D, Park C B. Estimation of near-surface shear-wave velocity by inversion of Rayleigh waves. Geophysics, 1999, 64(3): 691?700

[29] Lu L, Wang C, Zhang B. Inversion of multimode Rayleigh waves in the presence of a low-velocity layer: numerical and laboratory study. Geophysical Journal of the Royal Astronomical Society, 2010, 168(3): 1235?1246

[30] 宋先海, 肖柏勛, 黃榮榮, 等. 用等厚薄層權(quán)重自適應(yīng)迭代阻尼最小二乘法反演瑞雷波頻散曲線. 物探與化探, 2003, 16(4): 212?216

[31] Constable S C, Parker R L, Constable C G. Occam’s inversion: a practical algorithm for generating smooth models from electromagnetic sounding data. Geophy-sics, 1987, 52(3): 289?300

[32] Moro G D, Pipan M, Gabrielli P. Rayleigh wave dis-persion curve inversion via genetic algorithms and marginal posterior probability density estimation. Journal of Applied Geophysics, 2007, 61(1): 39?55

[33] Maraschini M, Foti S. A Monte Carlo multimodal in-version of surface waves. Geophysical Journal Inter-national, 2010, 182(3): 1557?1566

[34] Yamanaka H, Chimoto K. Variability of shallow soil amplification from surface-wave inversion using the Markov-chain Monte Carlo method. Soil Dynamics and Earthquake Engineering, 2018, 107: 141?151

[35] 石耀霖, 金文. 面波頻散反演地球內(nèi)部構(gòu)造的遺傳算法. 地球物理學報, 1995, 38(2): 189?198

[36] Tokeshi K, Harutoonian P, Leo C J, et al. Use of sur-face waves for geotechnical engineering applications in Western Sydney. Advances in Geosciences, 2013, 35: 37?44

[37] Yamanaka H, Ishida H. Application of genetic algori-thm to an inversion of surface-wave dispersion data. Bulletin of the Seismological Society of America, 2003, 86(2): 436?444

[38] Martínez M D, Lana X, Badal J, et al. Simulated an-nealing in 3-D seismic modelling and elastic structure of the mediterranean basin. Physics & Chemistry of the Earth Part A Solid Earth & Geodesy, 1999, 24 (3): 253?260

[39] Martínez M D, Lana X, Olarte J, et al. Inversion of Rayleigh wave phase and group velocities by simu-lated annealing. Physics of the Earth & Planetary In-teriors, 2000, 122(1): 3?17

[40] Beaty K S, Schmitt D R, Sacchi M. Simulated an-nealing inversion of multimode Rayleigh wave dis-persion curves for geological structure. Geophysical Journal International, 2002, 151(2): 622?631

[41] Song X, Gu H, Zhang X, et al. Pattern search al-gorithms for nonlinear inversion of high-frequency Rayleigh-wave dispersion curves. Computers and Geosciences, 2008, 34(6): 611?624

[42] Song X, Tang L I, Xiaochun L V, et al. Application of particle swarm optimization to interpret Rayleigh wave dispersion curves. Journal of Applied Geophy-sics, 2012, 84(9): 1?13

[43] Kranjcic D, Stumberger G. Differential evolution-based identification of the nonlinear Kaplan turbine model. IEEE Transactions on Energy Conversion, 2014, 29(1): 178?187

[44] 蔡偉, 宋先海, 袁士川, 等. 基于螢火蟲和蝙蝠群智能算法的瑞雷波頻散曲線反演. 地球物理學報, 2018, 61(6): 2409?2420

[45] 于東凱, 宋先海, 江東威, 等. 改進蜂群算法及其在面波頻散曲線反演中的應(yīng)用. 地球物理學報, 2018, 61(4): 1482?1495

[46] 蔡偉, 宋先海, 袁士川, 等. 新的瑞雷波多模式頻散曲線反演目標函數(shù). 地球科學, 2017, 42(9): 1608?1622

[47] 姚姚. 地球物理非線性反演模擬退火法的改進. 地球物理學報, 1995, 38(5): 643?650

[48] 王山山, 李燦平, 李青仁, 等. 快速模擬退火地震反演. 地球物理學報, 1995, 38(增刊1): 123?124

[49] 陳華根, 李麗華, 許惠平, 等. 改進的非?？焖倌M退火算法. 同濟大學學報(自然科學版), 2006, 34(8): 1121?1125

[50] 毛承英. 基于改進遺傳算法的瑞雷波頻散曲線反演[D]. 長沙: 中南大學, 2010: 30?44

[51] 閆月鋒, 孫成禹. 瑞雷波頻散曲線反演的線性約束問題研究//2014 年中國地球科學聯(lián)合學術(shù)年會——專題 19: 地震波傳播與成像論文集. 北京: 中國地球物理學會, 2014: 1167

[52] 閆月鋒. 石油勘探中地震面波的非線性反演方法研究[D]. 青島: 中國石油大學(華東), 2016: 57?61

[53] 崔建文. 一種改進的全局優(yōu)化算法及其在面波頻散曲線反演中的應(yīng)用. 地球物理學報, 2004, 47(3): 521?527

[54] Iglesias A, Atienza V M C, Shapiro N M, et al. Crustal structure of south-central Mexico estimated from the inversion of surface-wave dispersion curves using genetic and simulated annealing algorithms. Geophysical International, 2013, 40(3): 181?190

[55] 雷宇航. 自適應(yīng)GA與DLS聯(lián)合反演瑞雷波頻散曲線方法研究[D]. 西安: 西安石油大學, 2018: 30?62

[56] Senkaya M, Karsl? H. Joint inversion of Rayleigh-wave dispersion data and vertical electric sounding data: synthetic tests on characteristic sub-surface mo-dels. Geophysical Prospecting, 2016, 64(1): 228?246

[57] Calderon-Macias C, Luke B. Improved parameteriza-tion to invert Rayleigh-wave data for shallow profiles containing stiff inclusions. Geophysics, 2007, 72(1): U1?U10

[58] Menke W. Comment on ‘Dispersion function compu-tations for unlimited frequency values’ by Anas Abo-Zena. Geophysical Journal of the Royal Astronomi-cal Society, 1979, 59(2): 315?323

[59] 張曉陽, 郭珍, 唐拓, 等. 基于快速標量傳遞算法的頻散曲線正演模擬. 價值工程, 2015, 34(33): 130?132

[60] 李杰, 楊婧, 陳宣華. 面波頻散曲線遺傳算法反演的程序設(shè)計. 地球物理學進展, 2013, 28(5): 2693?2700

[61] 劉建宙. 有限差分法在瑞雷波波場正演模擬中的應(yīng)用[D]. 北京: 中國地質(zhì)大學(北京), 2014: 19?45

[62] Connolly D P, Giannopoulos A, Forde M C. A higher order perfectly matched layer formulation for finite-difference time-domain seismic wave modeling. Geo-physics, 2015, 80(1): T1?T16

Inversion Research of Rayleigh Wave Dispersion Curve Based on Fast Scalar Transfer Algorithm

DONG Zhikai1, DUAN Wensheng2, XIAO Chengwen2, HU Tianyue1,?, ZHANG Xianbing1,?

1. School of Earth and Space Sciences, Peking University, Beijing 100871; 2. Research Institute of Exploration and Development, Tarim Oilfield, PetroChina, Korla 841000; ? Corresponding authors, E-mail: tianyue@pku.edu.cn (HU Tianyue), zxb@pku.edu.cn (ZHANG Xianbing)

In order to improve the efficiency and accuracy of the inversion of Rayleigh wave dispersion curves near the surface, fast scalar transfer algorithm which has the characteristics of high computational efficiency is introduced to calculate the forward theoretical value of Rayleigh wave dispersion curve. The performances of genetic algorithm (GA), simulated annealing algorithm (SA) in the inversion of Rayleigh wave dispersion curves before and after adding linear constraints are compared. On this basis, linear constraints are added to GA and SA to improve the speed of convergence, and Monte Carlo method (MC) with fast computing speed is used to identify the types of formation as a supplementary means. Then the inversion results obtained by GA are taken as the initial state of SA as well as narrowing search scope appropriately, and this kind of joint inversion is carried out to overcome the premature problem of GA. Using the above method to calculate the three-layer model, noise-containing data and actual model of the work area. The results show that the method above is efficient, accurate and stable, and it has strong ability of global optimization and anti-noise ability to a certain extent.

shallow surface; Rayleigh wave; dispersion curve; fast scalar transfer algorithm; nonlinear inversion

10.13209/j.0479-8023.2020.020

國家科技重大專項(2016ZX05004003)和國家自然科學基金(41674122)資助

2019?05?10;

2019?06?30