劉輝，李靜,2，曾昭發(fā)，王天琪

(1.吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，吉林 長春 130021;2.地球信息探測儀器教育部重點(diǎn)實(shí)驗(yàn)室 (吉林大學(xué))，吉林 長春 130026)

0 引言

瑞利面波是一種由縱波與橫波干涉形成沿自由表面?zhèn)鞑サ牡卣鸩ā?887年，英國學(xué)者Rayleigh首先發(fā)現(xiàn)并證明均勻半空間中瑞利面波的存在[1-2]。面波頻散曲線反演是獲取淺地表地下空間、深部巖石圈和地球構(gòu)造橫波速度的重要方法技術(shù)[3-5]。Park等[6]提出的多道面波分析方法(MASW)廣泛應(yīng)用于淺地表面波數(shù)據(jù)采集處理；Xia等[7]采取奇異值分解和L-M方法，實(shí)現(xiàn)了瑞利波基階頻散曲線反演；Song等[8]使用Occam算法進(jìn)行一維頻散曲線反演后獲得了一個(gè)含低速層的擬二維橫波速度剖面；魯來玉等[9]和羅銀河等[10]分別采用遺傳算法和最小二乘方法實(shí)現(xiàn)了不同尺度多階面波頻散聯(lián)合反演，可有效提高單純只依靠基階面波的反演精度。

常規(guī)基于迭代的最小二乘線性反演方法依賴于初始模型且存在多極值、容易陷入局部最小、反演精度低等問題?；谌肿顑?yōu)的非線性反演方法由于不依賴于初始模型，已廣泛應(yīng)用于地球物理反演。Beaty等[11]利用模擬退火法在一定概率下可以接受劣解的搜索方式實(shí)現(xiàn)了多階瑞利波聯(lián)合反演；宋先海等將模式識別算法[12]和粒子群算法[13]應(yīng)用到瑞利波頻散曲線反演中；蔡偉等[14]將螢火蟲和蝙蝠群智能算法應(yīng)用于瑞利波頻散曲線反演；于東凱等分別利用改進(jìn)蜂群算法[15]和蚱蜢算法[16]實(shí)現(xiàn)了瑞利波頻散曲線反演。在天然地震被動源面波成像研究中，Aki[17]提出利用背景噪聲研究地下結(jié)構(gòu)，從臺陣信號中提取瑞利波頻散信息的空間自相關(guān)法(SPAC)；徐佩芬等[18]利用SPAC法從微動信號中提取頻散曲線并利用遺傳算法進(jìn)行反演獲得S波速度結(jié)構(gòu)；張寶龍等[19]利用擴(kuò)展空間自相關(guān)法(ESPAC)提取頻散曲線結(jié)合面波層析成像技術(shù)獲得了地殼淺層三維S波速度結(jié)構(gòu)。雖然傳統(tǒng)算法如遺傳算法具有對初始模型要求寬松、無需求梯度、易于實(shí)現(xiàn)并行化等優(yōu)點(diǎn)，較適合于解決非線性、多參數(shù)的瑞利波頻散曲線反演問題，但也存在計(jì)算量大、局部搜索能力差和收斂穩(wěn)定性差等缺陷，使其在頻散曲線反演中受到一定的限制。

針對上述算法中存在的問題，結(jié)合前人研究成果，本文提出一種基于貝葉斯馬爾科夫蒙特卡洛(MCMC)理論的隨機(jī)反演方法，通過先驗(yàn)信息建立解空間，采用馬爾科夫鏈擾動模擬算法，計(jì)算模型頻散曲線與實(shí)際曲線的似然函數(shù)，再利用轉(zhuǎn)移概率規(guī)則引導(dǎo)搜索過程，可得到與實(shí)際頻散曲線匹配最佳的最優(yōu)解，獲得大概率的橫波速度后驗(yàn)概率分布，從而減少反演的多解性，提高分辨率。本文通過典型理論模型和實(shí)測面波數(shù)據(jù)驗(yàn)證基于貝葉斯理論的隨機(jī)反演方法的可行性。實(shí)測數(shù)據(jù)應(yīng)用中，利用GPR解釋剖面進(jìn)行深度約束作為先驗(yàn)?zāi)Ｐ托畔?，在此基礎(chǔ)上開展面波隨機(jī)反演，從而獲得近地表橫波速度結(jié)構(gòu)。

1 貝葉斯隨機(jī)反演原理

根據(jù)貝葉斯理論可以用d表示實(shí)際數(shù)據(jù)，用m表示橫波(Vs)速度模型：

p(m│d)=p(d│m)p(m)/p(d),

∝p(d│m)p(m),

(1)

其中，p(m│d)是模型在給定數(shù)據(jù)的情況下的后驗(yàn)概率，p(m)是引入數(shù)據(jù)之前有關(guān)模型的先驗(yàn)信息，p(d│m)是似然函數(shù)(在給定特定模型(m)的情況下觀察到實(shí)測數(shù)據(jù)的概率)，p(d)是實(shí)際數(shù)據(jù)的準(zhǔn)確性，與m無關(guān)可以視為常數(shù)。

實(shí)際數(shù)據(jù)由N對頻率(f)和瑞利波相速度(PV)組成的頻散曲線以及偏差(σ)組成：

d=[f1f2…fN,PV1PV2…PVN] , (2)

σ=[σ1σ1…σN] 。

(3)

相速度偏差(σ)是表示所拾取的頻散曲線的不確定性的量，其大小可取頻散能量圖的各頻率的相速度能量集中的區(qū)間長度。

使用Voronoi核來表示模型參數(shù)[20]。如圖1所示?；疑珵関s模型空間，每一層都包含一個(gè)約束核和多個(gè)浮動核，在每個(gè)深度的vs值由離其最近的核決定，Voronoi核在模型空間中移動，其中浮動核可以在模型空間的深度范圍內(nèi)移動，而約束核只能在其固定的層移動，于是模型參數(shù)為：

m=[k,dp1,dp2,…,dpk,vs1,vs2,…,vsk,I,

dpc1,dpc2,…,dpcI,vsc2,…,vscI],

(4)

式中，k是浮動核的數(shù)量，I為約束核的數(shù)量，dpk為浮動核的深度，vsk是浮動核的橫波速度，dpci是約束核的深度，vsci是約束核的橫波速度。

先驗(yàn)信息p(m)可以分為兩部分：

p(m)=p(m|n)p(n) 。

(5)

其中,p(n)是Voronoi核數(shù)的先驗(yàn)，設(shè)區(qū)間I={n∈N|nmin

(6)

其中Δn=(nmax-nmin)，由于Voronoi核的深度c和速度v相互獨(dú)立，所以p(m|n)可拆分為：

p(m|n)=p(c|n)p(v|n) 。

(7)

速度的先驗(yàn)區(qū)間為J={vi∈R|vmin

(8)

a—無先驗(yàn)信息的vs均勻半空間模型;b—有先驗(yàn)信息的三層vs空間模型a—1 layer model without constraints;b—3 layer model with GPR constraints圖1 使用Voronoi核進(jìn)行模型參數(shù)化Fig.1 Model parameterization using Voronoi nuclei

其中,Δv=(vmax-vmin)。由于每個(gè)Voronoi核的速度相互獨(dú)立，則:

(9)

(10)

結(jié)合式(5)～(10)先驗(yàn)概率密度函數(shù)：

(11)

似然函數(shù)p(d│m)表示在給定模型m的情況下觀測到實(shí)際數(shù)據(jù)的可能性，實(shí)際上是通過對模型m進(jìn)行正演模擬得到頻散曲線G(m)并計(jì)算其與實(shí)際頻散曲線的相似度來實(shí)現(xiàn)的。假設(shè)數(shù)據(jù)di(i=1,2,…,Ndata)的每對頻率fi與相速度PVi相互獨(dú)立，則似然函數(shù)：

對每個(gè)新模型進(jìn)行正演計(jì)算以確定其是否達(dá)到可接受的標(biāo)準(zhǔn)：

式中：q(m′|m)是從模型m變化到模型m′的概率，q(m|m′)是從模型m′變化到模型m的概率，J為雅可比矩陣。對于不改變模型尺度，改變核的vs(圖2a)和改變核的深度(圖2b)，其擾動是對稱的，即從m變化到m′的概率等于從m′變化到m的概率:

a—改變核的vs;b—改變核的深度;c—產(chǎn)出一個(gè)浮動核;d—移除一個(gè)浮動核a—change vs of a nuclcus;b—move a nuclcus to a different depth;c—give birth to a new floating nucleus;d—remove a floating nucleus圖2 模型的4個(gè)擾動方式Fig.2 Illustration of four possible perturbations to a current model

(14)

(15)

所以這兩種擾動：

q(m′|m)=q(m|m′) 。

(16)

對于產(chǎn)生一個(gè)新的浮動核(圖2c)和移除一個(gè)浮動核(圖2d)，模型尺寸發(fā)生了變化，且擾動不對稱。由于Voronoi核的深度和速度相互獨(dú)立，對于產(chǎn)生新核：

(17)

(18)

對于移除浮動核：

(19)

式(13)中，J是從模型m變化到模型m′的雅可比矩陣，表示m與m′的尺度變化，即只有在兩個(gè)不同尺度的模型之間擾動(產(chǎn)生或移除Voronoi核)時(shí)才需計(jì)算雅克比矩陣，如果當(dāng)前模型與擾動模型具有相同尺度，則|J|=1，可以忽略。對于產(chǎn)生Voronoi核，從m到m′的雙射映射h：

(20)

(21)

vi是新核位置擾動前的速度值。模型空間分為離散空間(核深度)和連續(xù)空間(速度)，uc是用于離散空間轉(zhuǎn)換的離散變量，而uv是用于連續(xù)空間轉(zhuǎn)換的連續(xù)變量。Denison等[22]研究表明對于離散變換，雅克比項(xiàng)|J|=1，可以忽略，所以雅克比項(xiàng)僅考慮以下變量：

(22)

因此：

(23)

對于移除Voronoi核，|J|death=|J-1|birth=1，所以對于上述4種擾動，雅克比項(xiàng)都等于1，可以忽略。

在求得新模型的a值后，將其與隨機(jī)數(shù)u～U[0,1]進(jìn)行比較：如果a>u，則接受新模型m′并將其添加到鏈中，如果a

面波隨機(jī)反演的基本流程如圖3所示，其基本實(shí)現(xiàn)流程如下：

圖3 面波隨機(jī)反演基本流程Fig.3 Basic flowchart of surface wave stochastic inversion

1)根據(jù)地震面波資料提取頻散曲線；

2)根據(jù)先驗(yàn)信息建立解空間；

3)對解空間內(nèi)的隨機(jī)模型進(jìn)行預(yù)反演，其結(jié)果作為反演的初始模型m；

4)計(jì)算初始模型m的理論頻散曲線與似然函數(shù)；

5)隨機(jī)擾動產(chǎn)生新模型m′，并計(jì)算理論頻散曲線及似然函數(shù)；

6)計(jì)算a，若a>u(0,1),則將新模型m′加入到馬爾科夫鏈中，否則保留原模型m；

7)重復(fù)步驟5)、6)，直至滿足迭代停止條件。

2 面波頻散隨機(jī)反演模型測試

為驗(yàn)證本文所述隨機(jī)反演方法，分別采用4層vs遞增模型和4層含低速夾層模型開展反演測試。表1給出了4層vs遞增模型的基本參數(shù)。圖4為4層vs遞增模型基階頻散曲線的反演結(jié)果。圖4a為無約束反演結(jié)果，反演結(jié)果(最佳模型)與真實(shí)模型在深度7 m的第2、3層分界面處有所差異，vs的后驗(yàn)概率分布與真實(shí)模型基本擬合；圖4b為帶約束反演結(jié)果，vs的先驗(yàn)約束根據(jù)頻散曲線設(shè)為300～1 300m/s(vs的下邊界為0.9倍的相速度最小值；vs的上邊界為1.1倍的相速度最大值)。深度約束為3層，分別為0～2 m、3～15 m、16～30 m，相比于無約束反演結(jié)果，帶約束反演最佳模型和vs的后驗(yàn)概率密度分布與真實(shí)模型更好地吻合。圖4c為無約束反演的反演結(jié)果(最佳模型)與實(shí)際模型頻散曲線，兩條曲線在所有頻率范圍內(nèi)基本重合。圖4d為有約束反演的最佳模型與實(shí)際模型頻散曲線。圖4e為兩種模式下的迭代誤差曲線，有約束反演(紅線)迭代1 000次左右就能將誤差降低到無約束反演(藍(lán)線)迭代8 000次的誤差大小，且目標(biāo)函數(shù)最終的迭代誤差也更小。該模型測試驗(yàn)證了對于速度遞增模型，本文提出的隨機(jī)反演策略即使在無約束條件下也能得到較好的反演結(jié)果。

表1 4層vs遞增模型參數(shù)

a—無約束反演結(jié)果;b—有先驗(yàn)約束反演結(jié)果;c—無約束反演的最佳模型與實(shí)際模型頻散曲線;d—有先驗(yàn)約束反演的最佳模型與實(shí)際模型頻散曲線;e—迭代誤差曲線a—inversion results without constraints;b—inversion results with prior constraints;c—dispersion curves of the best fitting model and true model for unconstrained inversion;d—dispersion curves of the best fitting model and the true model with prior constraint inversion;e—misfit iteration curve圖4 4層vs遞增模型基階頻散曲線反演Fig.4 Fundamental dispersion curve inversion of four-layer vs increasing model

表2給出了4層含低速夾層測試模型的基本參數(shù)。圖5為4層含低速夾層模型基階頻散曲線的反演結(jié)果。圖5a為無約束反演結(jié)果，最佳模型(反演結(jié)果)與真實(shí)模型誤差較大；圖5b為有約束反演結(jié)果，此處vs的先驗(yàn)約束范圍同樣設(shè)為300～1 300 m/s，深度約束為3層，分別為0～2 m、3～17 m、18～30 m，最佳模型和vs的后驗(yàn)概率密度分布與真實(shí)模型基本吻合。圖5c為無約束反演的最佳模型與實(shí)際模型頻散曲線，兩曲線基本重合，最佳模型與真實(shí)模型的反演誤差主要來自于基階波對低速層的不敏感以及無約束隨機(jī)反演局部搜索能力差。圖5d為無約束反演頻散曲線誤差，可以看出在不同頻率范圍內(nèi)存在一定的頻散曲線擬合誤差。圖5e和圖5f為帶約束反演的最佳模型與實(shí)際模型頻散曲線和頻散曲線誤差，相比無約束反演，有約束反演誤差較小且整體幾乎等于零，表明了約束隨機(jī)反演具有較高的反演精度。圖5g為兩種模式下的迭代誤差曲線，無約束反演(藍(lán)線)在迭代8 000次左右時(shí)誤差達(dá)到最小且趨于穩(wěn)定，而有約束反演(紅線)迭代3 000次左右就能達(dá)到同樣的誤差大小且最終的迭代誤差也更小。上述模型測試表明，對于含低速夾層的速度結(jié)構(gòu)，由于面波頻散信息復(fù)雜，存在模態(tài)混疊現(xiàn)象，需要在一定約束條件下才能得到滿意的反演結(jié)果。

表2 4層含低速夾層模型參數(shù)

a—無約束反演結(jié)果;b—有先驗(yàn)約束反演結(jié)果;c—無約束反演的最佳模型與實(shí)際模型頻散曲線;d—無約束反演頻散曲線誤差;e—有先驗(yàn)約束反演的最佳模型與實(shí)際模型頻散曲線;f—帶約束反演頻散曲線誤差;g—迭代誤差曲線a—stochastic inversion results without constraints;b—stochastic inversion results with prior constraints;c—dispersion curve of the best fitting model and true model for unconstrained inversion;d—error of dispersion curve for unconstrained inversion;e—dispersion curves of the best fitting model and the true model with prior constraint inversion;f—error of dispersion curve with prior constraint inversion;g—misfit iteration curve圖5 4層含低速夾層模型基階頻散曲線反演Fig.5 Fundamental dispersion curve inversion of low-velocity layer model

3 實(shí)測地震數(shù)據(jù)測試

為進(jìn)一步驗(yàn)證隨機(jī)反演方法的可行性與可靠性，采用淺層主動源面波實(shí)測數(shù)據(jù)開展橫波速度反演測試。實(shí)測數(shù)據(jù)來自于中國安徽省淮南市某校園操場[23]，開展的地球物理方法包括淺層主動源面波與探地雷達(dá)(GPR)方法。面波數(shù)據(jù)沿測線共24炮，采用24個(gè)4 Hz檢波器，道間距為1 m，最小偏移距10 m，震源采用錘擊震源；GPR數(shù)據(jù)包含以第13炮集記錄測量點(diǎn)為中心的共中心點(diǎn)道集(CMP)記錄以及覆蓋整個(gè)地震測線的共偏移記錄。圖6a為處理探地雷達(dá)數(shù)據(jù)得到的GPR剖面。圖6b為基于卷積神經(jīng)網(wǎng)絡(luò)(CNN)[24]進(jìn)行GPR層界面自動識別結(jié)果，根據(jù)GPR剖面與CNN識別結(jié)果將地下劃分為4層，黑色實(shí)線為層界面。圖7a為經(jīng)過處理僅保留面波成分的第13炮地震記錄。圖7b為由炮集記錄提取的頻散能量，根據(jù)能量的最大值提取出基階頻散曲線。圖7c無約束反演結(jié)果，根據(jù)vs的后驗(yàn)概率分布可以看出，在10 m深度范圍內(nèi)地下共4層，其中第3層為低速夾層，深度為4 ～8 m，反演的最佳模型在淺層與vs后驗(yàn)概率分布基本匹配，但深層誤差較大。圖7d為GPR結(jié)果深度約束反演結(jié)果，vs約束根據(jù)頻散曲線設(shè)為100～500 m/s，與無約束反演相比，第3層的vs后驗(yàn)概率分布更加集中，且與反演的最佳模型更加匹配。圖 7e為迭代誤差曲線，可見在有約束條件下，隨著迭代次數(shù)的增加，GPR約束反演的收斂速度更快。

圖6 探地雷達(dá)處理剖面(a)與基于卷積神經(jīng)網(wǎng)絡(luò)的GPR層界面識別結(jié)果(b)Fig.6 Processed GPR profile(a)and GPR layer recognition results based on Convolutional Neural Networks(b)

將上述一維反演結(jié)果按照空間位置分布組合成擬二維剖面，如圖8所示。圖8a為無約束反演結(jié)果組合成的擬二維剖面，黑色實(shí)線為GPR分層界面，從圖中可以看出：面波頻散反演的橫波速度結(jié)果前兩層分界面與GPR數(shù)據(jù)基本匹配，但后3層分界面誤差較大。圖8b為GPR深度約束反演結(jié)果組合成的擬二維剖面，前3層分界面與GPR數(shù)據(jù)基本匹配，雖然底層界面仍然存在誤差，但相比無約束反演，誤差已經(jīng)有所降低。

結(jié)合前人鉆孔研究成果[25]，地下結(jié)構(gòu)分層解釋結(jié)論為：第1層由松軟的沙子、礫石、土壤組成，深度為1～2.2 m，橫波速度為100～150 m/s；第2層由緊湊且致密的沙子、礫石、土壤組成，深度為2.2～4.2 m，橫波速度為350～500 m/s；第3層為黏土，深度為4.2～8 m，橫波速度為120～180 m/s；基巖層由風(fēng)化程度不同的砂巖組成，深度大于8 m，橫波速度為400～500 m/s。前人鉆孔深度解釋與本文約束隨機(jī)反演橫波速度結(jié)構(gòu)有較好的吻合性。

a—第13炮地震記錄;b—由炮集記錄提取的頻散能量;c—無約束反演結(jié)果;d—GPR深度約束反演結(jié)果;e—迭代誤差曲線a—one shot gather of 13th shot;b—extracted dispersion curve;c—unconstrained stochastic inversion result;d—stochastic inversion with GPR constrained;e—iterations圖7 實(shí)測數(shù)據(jù)反演Fig.7 Inversion results of real surface wave data

a—無約束反演結(jié)果組合成的擬2D剖面;b—GPR深度約束反演結(jié)果組合成的擬2D剖面a—2D velocity profile of unconstrained stochastic inversion;b—2D velocity profile of GPR constrained stochastic inversion圖8 擬二維剖面反演結(jié)果對比Fig.8 Comparison of quasi-two-dimensional profile inversion results

4 結(jié)論

本文提出了基于貝葉斯反演理論的隨機(jī)反演方法并應(yīng)用于面波頻散曲線反演。利用先驗(yàn)信息深度約束分別對理論和實(shí)際面波數(shù)據(jù)進(jìn)行無約束、有約束隨機(jī)頻散反演對比測試，結(jié)果表明：

1)隨機(jī)反演方法是一種高效、穩(wěn)健的反演方法，具有較高的全局尋優(yōu)能力，且不依賴于初始模型。

2)帶約束的隨機(jī)反演可以有效且可靠地應(yīng)用于面波頻散曲線反演中，與傳統(tǒng)無約束反演相比，隨機(jī)反演可以很好地融合先驗(yàn)信息，減少多解性，提高反演效率與反演精度。

本文所使用的先驗(yàn)約束為GPR深度約束，而vs速度范圍約束是由頻散曲線確定的，如果實(shí)際測區(qū)有地質(zhì)資料可以確定地下各層的vs速度范圍，將其加入到先驗(yàn)約束中可以進(jìn)一步減少多解性。