董智開(kāi) 段文勝 肖承文 胡天躍,? 張獻(xiàn)兵,?

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

董智開(kāi)1段文勝2肖承文2胡天躍1,?張獻(xiàn)兵1,?

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

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

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

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

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

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

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

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

1.1 正演模型

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

1.2 快速標(biāo)量傳遞法正演瑞雷波頻散曲線

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

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

重新定義Menke65B9法的矢量:

其中,

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

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

2 非線性反演算法

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

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

2.1 MC 輔助識(shí)別地層類型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3)選擇。將選擇概率(m)作為標(biāo)準(zhǔn), 選出優(yōu)良個(gè)體作為父本繁殖下一代。(m)的計(jì)算公式為

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

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

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

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

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

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

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

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

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

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

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

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

3.1 無(wú)噪聲理論模型計(jì)算及不同方法反演誤差對(duì)比

在加入線性約束條件的 GA 反演結(jié)果基礎(chǔ)上, 將搜索范圍適當(dāng)?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é)果附近進(jìn)行SA搜索。

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

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

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

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

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

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

利用 MC 判定出地層類型后, 將約束范圍設(shè)定為適用于對(duì)應(yīng)地層類型的數(shù)值進(jìn)行反演。計(jì)算得到含高速硬夾層模型和含低速軟夾層模型頻散曲線的最終擬合誤差分別為 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ù)的相對(duì)誤差在 1%以內(nèi), 表明該方法對(duì)不同類型地層的模型均有較好的反演效果。

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

3.2 含噪聲理論模型計(jì)算

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

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

4 實(shí)際模型計(jì)算實(shí)例

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

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

圖13 垂直 V 分量記錄

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

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

5 結(jié)論

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

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

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

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

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

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

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

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

[3] 黃忠賢, 李紅誼, 胥頤. 中國(guó)西部及鄰區(qū)巖石圈 S波速度結(jié)構(gòu)面波層析成像. 地球物理學(xué)報(bào), 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] 金聰, 楊文海, 羅登貴, 等. 面波頻散曲線提取方法對(duì)比分析. 地球物理學(xué)進(jìn)展, 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年中國(guó)地球科學(xué)聯(lián)合學(xué)術(shù)年會(huì)——專題 14: 地下介質(zhì)結(jié)構(gòu)及其變化的地震面波、背景噪聲及尾波研究論文集. 北京: 中國(guó)地球物理學(xué)會(huì), 2014: 762

[13] 羅銀河, 夏江海, 劉江平, 等. 基階與高階瑞利波聯(lián)合反演研究. 地球物理學(xué)報(bào), 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] 李欣欣, 李慶春. 利用改進(jìn)的 F-K 變換法提取瑞雷波的頻散曲線. 地球物理學(xué)進(jìn)展, 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)用于瑞雷波頻散曲線提取研究. 地震工程與工程振動(dòng), 2014, 34(5): 21?27

[20] 楊振濤, 陳曉非, 潘磊, 等. 基于矢量波數(shù)變換法(VWTM)的多道 Rayleigh 波分析方法. 地球物理學(xué)報(bào), 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ì)中軸對(duì)稱柱面瑞利面波頻散函數(shù)的計(jì)算. 地震工程與工程振動(dòng), 2001, 21(3): 1?5

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

[27] 凡友華, 劉家琦, 肖柏勛. 計(jì)算瑞利波頻散曲線的快速矢量傳遞算法. 湖南大學(xué)學(xué)報(bào)(自然科學(xué)版), 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)造的遺傳算法. 地球物理學(xué)報(bào), 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] 蔡偉, 宋先海, 袁士川, 等. 基于螢火蟲(chóng)和蝙蝠群智能算法的瑞雷波頻散曲線反演. 地球物理學(xué)報(bào), 2018, 61(6): 2409?2420

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

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

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

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

[49] 陳華根, 李麗華, 許惠平, 等. 改進(jìn)的非常快速模擬退火算法. 同濟(jì)大學(xué)學(xué)報(bào)(自然科學(xué)版), 2006, 34(8): 1121?1125

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

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

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

[53] 崔建文. 一種改進(jìn)的全局優(yōu)化算法及其在面波頻散曲線反演中的應(yīng)用. 地球物理學(xué)報(bào), 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]. 西安: 西安石油大學(xué), 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] 張曉陽(yáng), 郭珍, 唐拓, 等. 基于快速標(biāo)量傳遞算法的頻散曲線正演模擬. 價(jià)值工程, 2015, 34(33): 130?132

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

[61] 劉建宙. 有限差分法在瑞雷波波場(chǎng)正演模擬中的應(yīng)用[D]. 北京: 中國(guó)地質(zhì)大學(xué)(北京), 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

國(guó)家科技重大專項(xiàng)(2016ZX05004003)和國(guó)家自然科學(xué)基金(41674122)資助

2019?05?10;

2019?06?30