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

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

0 引言

瑞利面波是一種由縱波與橫波干涉形成沿自由表面?zhèn)鞑サ牡卣鸩ā?887年，英國(guó)學(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è)含低速層的擬二維橫波速度剖面；魯來(lái)玉等[9]和羅銀河等[10]分別采用遺傳算法和最小二乘方法實(shí)現(xiàn)了不同尺度多階面波頻散聯(lián)合反演，可有效提高單純只依靠基階面波的反演精度。

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

針對(duì)上述算法中存在的問(wèn)題，結(jié)合前人研究成果，本文提出一種基于貝葉斯馬爾科夫蒙特卡洛(MCMC)理論的隨機(jī)反演方法，通過(guò)先驗(yàn)信息建立解空間，采用馬爾科夫鏈擾動(dòng)模擬算法，計(jì)算模型頻散曲線與實(shí)際曲線的似然函數(shù)，再利用轉(zhuǎn)移概率規(guī)則引導(dǎo)搜索過(guò)程，可得到與實(shí)際頻散曲線匹配最佳的最優(yōu)解，獲得大概率的橫波速度后驗(yàn)概率分布，從而減少反演的多解性，提高分辨率。本文通過(guò)典型理論模型和實(shí)測(cè)面波數(shù)據(jù)驗(yàn)證基于貝葉斯理論的隨機(jī)反演方法的可行性。實(shí)測(cè)數(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í)測(cè)數(shù)據(jù)的概率)，p(d)是實(shí)際數(shù)據(jù)的準(zhǔn)確性，與m無(wú)關(guān)可以視為常數(shù)。

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

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

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

(3)

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

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

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

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

(4)

式中，k是浮動(dòng)核的數(shù)量，I為約束核的數(shù)量，dpk為浮動(dòng)核的深度，vsk是浮動(dòng)核的橫波速度，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—無(wú)先驗(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的情況下觀測(cè)到實(shí)際數(shù)據(jù)的可能性，實(shí)際上是通過(guò)對(duì)模型m進(jìn)行正演模擬得到頻散曲線G(m)并計(jì)算其與實(shí)際頻散曲線的相似度來(lái)實(shí)現(xiàn)的。假設(shè)數(shù)據(jù)di(i=1,2,…,Ndata)的每對(duì)頻率fi與相速度PVi相互獨(dú)立，則似然函數(shù)：

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

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

a—改變核的vs;b—改變核的深度;c—產(chǎn)出一個(gè)浮動(dòng)核;d—移除一個(gè)浮動(dòng)核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è)擾動(dòng)方式Fig.2 Illustration of four possible perturbations to a current model

(14)

(15)

所以這兩種擾動(dòng)：

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

(16)

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

(17)

(18)

對(duì)于移除浮動(dòng)核：

(19)

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

(20)

(21)

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

(22)

因此：

(23)

對(duì)于移除Voronoi核，|J|death=|J-1|birth=1，所以對(duì)于上述4種擾動(dòng)，雅克比項(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)對(duì)解空間內(nèi)的隨機(jī)模型進(jìn)行預(yù)反演，其結(jié)果作為反演的初始模型m；

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

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

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

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

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

為驗(yàn)證本文所述隨機(jī)反演方法，分別采用4層vs遞增模型和4層含低速夾層模型開展反演測(cè)試。表1給出了4層vs遞增模型的基本參數(shù)。圖4為4層vs遞增模型基階頻散曲線的反演結(jié)果。圖4a為無(wú)約束反演結(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，相比于無(wú)約束反演結(jié)果，帶約束反演最佳模型和vs的后驗(yàn)概率密度分布與真實(shí)模型更好地吻合。圖4c為無(wú)約束反演的反演結(jié)果(最佳模型)與實(shí)際模型頻散曲線，兩條曲線在所有頻率范圍內(nèi)基本重合。圖4d為有約束反演的最佳模型與實(shí)際模型頻散曲線。圖4e為兩種模式下的迭代誤差曲線，有約束反演(紅線)迭代1 000次左右就能將誤差降低到無(wú)約束反演(藍(lán)線)迭代8 000次的誤差大小，且目標(biāo)函數(shù)最終的迭代誤差也更小。該模型測(cè)試驗(yàn)證了對(duì)于速度遞增模型，本文提出的隨機(jī)反演策略即使在無(wú)約束條件下也能得到較好的反演結(jié)果。

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

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

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

a—無(wú)約束反演結(jié)果;b—有先驗(yàn)約束反演結(jié)果;c—無(wú)約束反演的最佳模型與實(shí)際模型頻散曲線;d—無(wú)約束反演頻散曲線誤差;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í)測(cè)地震數(shù)據(jù)測(cè)試

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

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

將上述一維反演結(jié)果按照空間位置分布組合成擬二維剖面，如圖8所示。圖8a為無(wú)約束反演結(jié)果組合成的擬二維剖面，黑色實(shí)線為GPR分層界面，從圖中可以看出：面波頻散反演的橫波速度結(jié)果前兩層分界面與GPR數(shù)據(jù)基本匹配，但后3層分界面誤差較大。圖8b為GPR深度約束反演結(jié)果組合成的擬二維剖面，前3層分界面與GPR數(shù)據(jù)基本匹配，雖然底層界面仍然存在誤差，但相比無(wú)約束反演，誤差已經(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—無(wú)約束反演結(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í)測(cè)數(shù)據(jù)反演Fig.7 Inversion results of real surface wave data

a—無(wú)約束反演結(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é)果對(duì)比Fig.8 Comparison of quasi-two-dimensional profile inversion results

4 結(jié)論

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

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

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

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