王升超，韓立國(guó)，鞏向博，張盼

吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

馬爾科夫鏈蒙特卡洛方法(MCMC)是基于數(shù)學(xué)概率分布模擬的一種方法，目前已經(jīng)廣泛應(yīng)用于理論物理、信號(hào)通信、醫(yī)學(xué)等領(lǐng)域，在地球物理領(lǐng)域的應(yīng)用也在不斷的發(fā)展中，特別是應(yīng)用于地球物理反演問(wèn)題(Grandis et al.,1999;Malinverno and Leaney,2000).相比于其他反演方法，MCMC方法的優(yōu)點(diǎn)在于全局尋優(yōu)，不會(huì)陷入局部極小，并且不依賴于準(zhǔn)確的先驗(yàn)?zāi)Ｐ?，可以引入更為?fù)雜的先驗(yàn)信息.

對(duì)于地球物理反演問(wèn)題，基于最小二乘的方法得到了廣泛的應(yīng)用，Menke(1989)提出了確定性最小二乘法，而Tarantola和Valette(1982)提出了最小二乘的概率反演方法，在貝葉斯框架下，反演問(wèn)題的解是一個(gè)概率密度函數(shù)，稱之為后驗(yàn)概率密度函數(shù).Hansen等(2006)提出了一種基于連續(xù)模擬的線性走時(shí)層析成像的概率反演方法，他們利用了經(jīng)典最小二乘反演和克里格法(Journel and Huijbregts，1978)的等價(jià)性.Nielsen等(2010)給出了這種方法在井間探地雷達(dá)數(shù)據(jù)中的應(yīng)用.Gloaguen等(2005a，b)提出了一種基于克里格法的誤差模擬的相關(guān)方法，相當(dāng)于概率最小二乘法，并將其應(yīng)用于井間探地雷達(dá)層析成像.Giroux和Gloaguen(2012)將這種方法用于各向異性速度場(chǎng)的反演.這些方法僅對(duì)嚴(yán)格線性反演問(wèn)題有效，并且依賴于描述噪聲模型和先驗(yàn)?zāi)Ｐ偷母咚菇y(tǒng)計(jì)的固有假設(shè)，先驗(yàn)?zāi)Ｐ捅仨氁愿咚剐问浇o出，先驗(yàn)?zāi)Ｐ陀删岛蛥f(xié)方差模型定義.國(guó)內(nèi)，張廣智等(2011)將MCMC方法應(yīng)用于疊前地震數(shù)據(jù)反演，在波阻抗反演方面取得了一定效果.殷長(zhǎng)春等(2014)將概率反演的方法應(yīng)用到航空電磁領(lǐng)域，克服了初始模型的影響，成功反演了深度低阻層.王朋巖等(2015)成功使用MCMC方法反演了巖性參數(shù).孫月成(2018)則將基于MCMC 算法的地質(zhì)統(tǒng)計(jì)學(xué)反演應(yīng)用在油藏模擬當(dāng)中，取得了良好的效果.

但在實(shí)際反演問(wèn)題當(dāng)中，先驗(yàn)信息往往比高斯模型描述的更為復(fù)雜，而且先驗(yàn)信息也不能用準(zhǔn)確的數(shù)學(xué)公式表達(dá).Hansen等(2008b)展示了拓展的Metropolis算法(Mosegaard et al.，1995)在非線性井間層析成像問(wèn)題中的應(yīng)用，其中先驗(yàn)?zāi)Ｐ褪欠歉咚沟模⒖梢杂扇魏蔚刭|(zhì)統(tǒng)計(jì)學(xué)方法定義.擴(kuò)展的Metropolis算法不需要先驗(yàn)信息的準(zhǔn)確數(shù)學(xué)表達(dá)式，可以用來(lái)采樣后驗(yàn)概率密度函數(shù)，包括高度非線性的反演問(wèn)題.在拓展的Metropolis算法中，只要有一個(gè)能夠?qū)ο闰?yàn)概率密度函數(shù)進(jìn)行采樣的“黑箱”算法，就可以完成對(duì)先驗(yàn)信息的采樣任務(wù).Hansen等(2008a,2012)提出了連續(xù)吉布斯采樣的方法，連續(xù)吉布斯算法在采樣先驗(yàn)信息時(shí)，可以有效的控制擾動(dòng)的步長(zhǎng)，具有很高的采樣效率(Gómez-Hernández and Journel，1993).這樣，將連續(xù)吉布斯采樣作為一種“黑箱”算法應(yīng)用于拓展的Metropolis算法中，使得在進(jìn)行概率反演時(shí)，對(duì)先驗(yàn)?zāi)Ｐ偷倪x擇變得非常靈活(Journel and Zhang，2006)，可以對(duì)任何地質(zhì)統(tǒng)計(jì)學(xué)算法定義的先驗(yàn)信息進(jìn)行采樣，既可用于相對(duì)簡(jiǎn)單的兩點(diǎn)統(tǒng)計(jì)模型(例如基于高斯的先驗(yàn)?zāi)Ｐ?，也可以用于更復(fù)雜的多點(diǎn)的統(tǒng)計(jì)模型，通過(guò)多點(diǎn)算法可以靈活地模擬高熵、低熵結(jié)構(gòu)或其組合的實(shí)現(xiàn)，使得在求解概率反演時(shí)，可以引入更復(fù)雜的先驗(yàn)?zāi)Ｐ托畔?Hansen et al.,2013).

后驗(yàn)概率分布的求解是通過(guò)模型提供的先驗(yàn)信息，數(shù)據(jù)的不確定性和模型參數(shù)等綜合求解.其中所有可用的先驗(yàn)信息都由先驗(yàn)概率密度函數(shù)描述.多種信息綜合反映了后驗(yàn)樣本模型的變化，為后驗(yàn)采樣提供了依據(jù).概率反演方法除了能提供后驗(yàn)?zāi)Ｐ蛥?shù)的協(xié)方差外，也可以解決復(fù)雜的問(wèn)題，如地質(zhì)連通性的概率或流體的停留時(shí)間(Mosegaard,1998).

作者將討論使用MCMC方法對(duì)井間探地雷達(dá)(GPR)數(shù)據(jù)進(jìn)行初至走時(shí)的時(shí)移反演問(wèn)題，探地雷達(dá)(GPR)井間層析成像也是近地表地質(zhì)構(gòu)造和地球物理參數(shù)層析成像的常用方法.這種走時(shí)數(shù)據(jù)對(duì)電磁波速的地下變化很敏感，這與介電常數(shù)有關(guān)，而介電常數(shù)受地下含水量的強(qiáng)烈影響(Topp et al.，1980).因此，時(shí)移探地雷達(dá)反演結(jié)果可以反映地下含水量在不同時(shí)間的變化.

時(shí)移反演有兩種反演策略，稱為連續(xù)反演法和雙差法(Waldhauser and Ellsworth，2000)，Huang等(2020)在雙差法的基礎(chǔ)上實(shí)現(xiàn)了基于目標(biāo)導(dǎo)向的全波形時(shí)移反演，提高了時(shí)移反演的計(jì)算效率.兩種方法的主要區(qū)別在于雙差法減小了非目標(biāo)區(qū)的影響，提高了目標(biāo)區(qū)的反演精度.在時(shí)移反演中，作者使用雙差法結(jié)合MCMC方法對(duì)探地雷數(shù)據(jù)進(jìn)行反演.對(duì)于時(shí)移探地雷達(dá)數(shù)據(jù)，要解決目標(biāo)位置的變化，需要對(duì)兩組數(shù)據(jù)在不同觀測(cè)點(diǎn)進(jìn)行至少兩次反演計(jì)算.兩次反演計(jì)算量大，非目標(biāo)區(qū)域的反演也會(huì)影響目標(biāo)區(qū)域的精度.為了解決這個(gè)問(wèn)題，作者采用了局部采樣的MCMC反演方法，在雙差法的第二次反演過(guò)程中，采用局部采樣即只對(duì)目標(biāo)區(qū)域采樣的方法代替了全局采樣法.這樣，減少了非采樣區(qū)域?qū)r(shí)移位置的影響，在提高計(jì)算效率的同時(shí)，增加了目標(biāo)區(qū)域的反演精度.

本文中，作者演示了如何通過(guò)使用拓展的Metropolis算法結(jié)合使用連續(xù)吉布斯采樣的地質(zhì)統(tǒng)計(jì)學(xué)算法定義的先驗(yàn)信息對(duì)時(shí)移的探底雷達(dá)數(shù)據(jù)實(shí)現(xiàn)MCMC的反演.將MCMC方法應(yīng)用于探地雷達(dá)時(shí)移反演中，首先利用連續(xù)吉布斯采樣和Metropolis算法對(duì)反演問(wèn)題進(jìn)行求解.然后，結(jié)合雙差法得到目標(biāo)區(qū)域的變化.作者將該方法應(yīng)用到GPR模擬數(shù)據(jù)中，測(cè)試該方法的效果，分析對(duì)比了局部采樣的MCMC反演和全局采樣的MCMC的反演誤差，證明了局部采樣MCMC反演的有效性.

1 MCMC反演

1.1 概率反演公式

對(duì)于地球物理反演問(wèn)題，地下構(gòu)造可以用一組模型參數(shù)m來(lái)表示，觀測(cè)數(shù)據(jù)可以用一組數(shù)據(jù)d來(lái)表示.正演問(wèn)題是指通過(guò)一個(gè)函數(shù)映射f來(lái)獲取觀測(cè)數(shù)據(jù)d(Tarantola and Valette，1982)：

d=f(m)，

(1)

函數(shù)f通常基于相應(yīng)的物理關(guān)系，對(duì)應(yīng)的反演問(wèn)題的公式可以表示為

m=f-1(d)，

(2)

求解反演問(wèn)題的主要困難是如何而求解反演算子f-1，實(shí)際操作中會(huì)出現(xiàn)反演算子難以求解或者不存在的情況，此外，正演算子f是基于對(duì)正確物理關(guān)系的某種近似，有一定的誤差.在GPR反演的研究中，模型m代表了地下速度的層析成像結(jié)果，觀測(cè)數(shù)據(jù)是波在發(fā)射源和接收器之間的走時(shí).

反演問(wèn)題是通過(guò)先驗(yàn)信息來(lái)獲取模型參數(shù)，基于概率方法的反演公式可以表示為

σM(m)=kρM(m)L(m)，

(3)

反演問(wèn)題的解σM(m)是后驗(yàn)概率分布，其中k是歸一化因子.先驗(yàn)概率密度ρM(m)描述了與數(shù)據(jù)無(wú)關(guān)的模型參數(shù)先驗(yàn)信息.L(m)是似然函數(shù)，它是一種概率的度量，用于衡量給定模型相關(guān)參數(shù)與給定觀測(cè)數(shù)據(jù)不確定性模型的匹配程度，具體公式為

(4)

ρD(g(m))描述測(cè)量的不確定性，通常與記錄數(shù)據(jù)的儀器中的不確定性有關(guān).θ(d|m)表示由于使用不完善的正演方法或不完善的參數(shù)化而引起的建模誤差.μD(d)描述信息的統(tǒng)一性，以確保參數(shù)化對(duì)坐標(biāo)系中的更改保持不變.在大多數(shù)情況下，我們可以假設(shè)μD(d)是一個(gè)常數(shù).

1.2 基于射線的正演方法

求解探地雷達(dá)(GPR)數(shù)據(jù)初至走時(shí)的反演問(wèn)題時(shí)，正演的觀測(cè)數(shù)據(jù)是走時(shí)，也就是波在發(fā)射源和接收器之間的傳播時(shí)延.關(guān)于走時(shí)的計(jì)算，有幾種類型的正演方法，我們使用基于程函方程的方法(Hassouna and Farag,2007).在程函方程中，沿曲線u(x)的到達(dá)時(shí)間，以速度場(chǎng)m(x)定義傳播速度(Sethian and Popovici，1999)：

(5)

快速推進(jìn)法是一種具有高精度、高效率、無(wú)條件穩(wěn)定等特性的走時(shí)計(jì)算方法(孫章慶，2008)，利用多階快速推進(jìn)法可以有效的求解程函方程.這種正演模型是非線性的，因?yàn)槌毯匠虒?duì)應(yīng)于波動(dòng)方程的高頻近似，通常被稱為高頻射線近似.信號(hào)源和接收器之間的走時(shí)d可由式(6)給出：

(6)

其中G(x)是靈敏度核，它描述了每個(gè)模型參數(shù)(在Fresnell區(qū)內(nèi))對(duì)走時(shí)的靈敏度.G(x)可以在廣泛的假設(shè)下計(jì)算，我們利用波動(dòng)方程的高頻近似來(lái)計(jì)算靈敏度核G(x)，靈敏度核可以用連接源和接收器的射線來(lái)描述，這種類型的正演模型稱為基于射線的正演模型.

1.3 拓展的Metropolis 算法

在現(xiàn)實(shí)中，大多數(shù)地球物理反演問(wèn)題都是非線性的，而且用非高斯統(tǒng)計(jì)來(lái)描述.我們需要一種不需要先驗(yàn)概率密度顯式表達(dá)式的算法.拓展的Metropolis算法可以通過(guò)使用連續(xù)的吉布斯采樣作為黑箱算法來(lái)解決這個(gè)問(wèn)題，黑箱算法能夠在先驗(yàn)概率密度足夠的情況下執(zhí)行隨機(jī)游走(Mosegaard and Tarantola，1995).拓展的Metropolis算法包含兩個(gè)主要的步驟：

(1)提出了一個(gè)候選模型mpro，它是在當(dāng)前模型mcur中給出一個(gè)擾動(dòng)后產(chǎn)生的新模型，同時(shí)也是先驗(yàn)概率密度函數(shù)的一次實(shí)現(xiàn).

(2)決定接受或是拒絕當(dāng)前的模型，模型的接受依據(jù)Metropolis算法的接受概率：

(7)

式中，L(mpro)/L(mcur)為候選模型與當(dāng)前模型之間的似然函數(shù)之比.如果接受，則候選模型取代當(dāng)前模型，即達(dá)成了一個(gè)后驗(yàn)概率密度采樣的實(shí)現(xiàn).否則，候選模型將被拒絕，當(dāng)前模型將再次循環(huán)，所以在每次迭代中，后驗(yàn)概率密度的樣本量都會(huì)增加.然而，基于MCMC的采樣方法計(jì)算量大，對(duì)于時(shí)移反演問(wèn)題，需要對(duì)兩個(gè)時(shí)刻分別進(jìn)行反演，進(jìn)一步增加了計(jì)算量，并且非標(biāo)區(qū)域的反演會(huì)影響目標(biāo)區(qū)域的精度.針對(duì)這一問(wèn)題，作者使用了局部采樣的MCMC反演方法，在時(shí)移反演的第二次反演時(shí)，只對(duì)目標(biāo)區(qū)域進(jìn)行采樣，減少了計(jì)算量，提高了目標(biāo)區(qū)域的反演精度.

局部采樣的MCMC反演將依靠作為擴(kuò)展Metropolis算法一部分的連續(xù)的吉布斯采樣器完成.我們使用連續(xù)的吉布斯采樣器對(duì)ρM(m)直接采樣.流程如下：

(2)計(jì)算當(dāng)前模型和候選模型的似然函數(shù)，L(mcur)和L(mpro).

(4)如果建議的模型被接受，那么用mpro代替當(dāng)前模型mcur的轉(zhuǎn)換被接受，即mcur=mpro，否則，仍然在mcur位置，再次產(chǎn)生一個(gè)隨機(jī)擾動(dòng)進(jìn)行下一次迭代.

2 雙差法時(shí)移反演

本文主要采用雙差法來(lái)解決時(shí)移反演問(wèn)題，在雙差法的第二次反演時(shí)，采用局部采樣的MCMC方法進(jìn)行反演，該雙差法的主要流程如圖1所示.

圖1 雙差法時(shí)移反演流程圖Fig.1 Double difference strategy

3 合成的GPR數(shù)據(jù)測(cè)試

為了驗(yàn)證該方法的有效性，我們模擬了一個(gè)合成的探地雷達(dá)數(shù)據(jù)，由一個(gè)平均速度為1.45×10-1m·ns-1的多變量高斯概率分布生成的7 m×13 m(87×47像素，0.15 m×0.15 m)的參考模型mT1，方差為3×10-4m2ns-2，其球形協(xié)方差模型具有6 m的各向同性范圍，這與Looms等(2010)和Hansen等(2013)的觀測(cè)系統(tǒng)和參數(shù)設(shè)置是相似的.圖2顯示了探地雷達(dá)(GPR)鉆孔記錄的觀測(cè)系統(tǒng)，其中紅色的點(diǎn)為發(fā)射源，黑色圓點(diǎn)為接收器，總共ND=702對(duì).

圖2 數(shù)據(jù)觀測(cè)系統(tǒng)分布，共702組發(fā)射源和接收點(diǎn)組合Fig.2 Recording geometry,ND=702 pairs of sources (red crosses)and receivers (black dots)are represented by a connecting black line

圖3 起始速度模型mT1，監(jiān)測(cè)速度模型mT2和速度擾動(dòng)Fig.3 Initial model mT1,monitor model mT2 and the perturbation

使用速度模型mT1和mT2正演模擬獲得相應(yīng)的初至?xí)r間數(shù)據(jù)data1 和 data2作為觀測(cè)數(shù)據(jù)，圖4顯示了相應(yīng)的正演數(shù)據(jù)分布，其中黑色曲線代表data1，紅色曲線代表data2.由圖4可以看出，兩組正演數(shù)據(jù)的變化趨勢(shì)大體相同，只在局部位置有差異，可以認(rèn)為是速度擾動(dòng)區(qū)域在正演數(shù)據(jù)中產(chǎn)生的影響.

圖4 正演模擬的初至走時(shí)數(shù)據(jù)Fig.4 The simulated traveltime data1 and data2

在反演之前，首先要根據(jù)先驗(yàn)信息對(duì)先驗(yàn)概率密度函數(shù)進(jìn)行采樣，先驗(yàn)采樣是一個(gè)隨機(jī)過(guò)程，任何滿足先驗(yàn)條件的模型都有可能出現(xiàn).圖5顯示了對(duì)先驗(yàn)概率密度函數(shù)進(jìn)行連續(xù)吉布斯采樣后，得到的五個(gè)獨(dú)立實(shí)現(xiàn)的先驗(yàn)速度模型，由圖5可以看出，每個(gè)先驗(yàn)?zāi)Ｐ偷乃俣确植几鞑幌嗤?

圖5 先驗(yàn)概率密度函數(shù)的五個(gè)獨(dú)立隨機(jī)采樣結(jié)果Fig.5 Five statistically independent realizations of the priori probability density

在時(shí)移反演中，為了求得時(shí)移位置的變化，我們首先對(duì)起始時(shí)刻T1的觀測(cè)數(shù)據(jù)data1進(jìn)行MCMC反演，用拓展的Metropolis算法采樣獲得后驗(yàn)概率密度函數(shù)的樣本.圖6展示了五個(gè)后驗(yàn)采樣樣本，樣本中可以明顯看出速度分布淺層主要為低速的藍(lán)色，深層為高速的黃色，與給出的速度模型mT1特征一致.

圖6 后驗(yàn)概率密度函數(shù)的五個(gè)采樣樣本Fig.6 Five statistically independent realizations of the posterior probability density

為進(jìn)一步檢驗(yàn)反演的準(zhǔn)確性，我們對(duì)這些后驗(yàn)采樣得到的模型進(jìn)行正演模擬，正演的走時(shí)數(shù)據(jù)如圖7所示，其中紅色曲線為觀測(cè)數(shù)據(jù)data1，黑色曲線為多個(gè)后驗(yàn)采樣模型的正演數(shù)據(jù).可以發(fā)現(xiàn)，通過(guò)拓展的Metropolis算法獲得的后驗(yàn)樣本，其正演的結(jié)果與觀測(cè)數(shù)據(jù)data1有很高的一致性，證明了反演的方法的有效性.

圖7 后驗(yàn)樣本正演得到的走時(shí)數(shù)據(jù)Fig.7 The traveltime data of data and posterior

雙差法時(shí)移反演中，第二次反演求解速度模型mT2時(shí)，需要選擇第一次反演的結(jié)果作為初始模型，我們選擇后驗(yàn)樣本的一個(gè)結(jié)果作為起始模型進(jìn)行第二次MCMC反演.不同于第一次反演時(shí)的全采樣方法，第二次反演時(shí)我們只對(duì)目標(biāo)區(qū)域使用連續(xù)的吉布斯算法進(jìn)行采樣.為對(duì)比局部采樣的反演效果，我們同時(shí)設(shè)置了一組模型全采樣的數(shù)據(jù)進(jìn)行對(duì)比，具體結(jié)果如圖8所示.

圖8a顯示了使用全采樣的擴(kuò)展Metropolis算法采樣后驗(yàn)概率密度函數(shù)得到的五個(gè)速度模型，用圖8a的后驗(yàn)?zāi)Ｐ团c第一次反演得到的初始模型相減，得到擾動(dòng)變化如圖8b所示.圖8c是用局部采樣的MCMC反演得到的后驗(yàn)分布，圖8d是用圖8c與初始反演模型相減的結(jié)果.對(duì)比圖8a、c和模型mT2，可以看出圖8a、c在淺層位置都反演出了一定的高速分布，但是圖8a淺層高速分布的位置范圍變化較大，高速位置分布隨機(jī)并不集中在目標(biāo)區(qū)域，圖8c則主要集中在目標(biāo)區(qū)域位置.用圖8b、d對(duì)比更加明顯，全采樣的反演相對(duì)初始模型的變化位置，除淺層外，深層也有變化.

圖8 (a)T2時(shí)刻全采樣反演結(jié)果；(b)全采樣擾動(dòng)變化；(c)T2時(shí)刻局部采樣反演結(jié)果；(d)局部采樣擾動(dòng)變化Fig.8 (a)Full sampling inversion model of T2;(b)The perturbation of the inversion;(c)Local sampling inversion model of T2;(d)The perturbation of the local sampling inversion

圖9 (a)全采樣法的目標(biāo)區(qū)域反演結(jié)果；(b)局部采樣法的目標(biāo)區(qū)域反演結(jié)果Fig.9 (a)Target area of all sampling method;(b)Target area of local sampling method

圖10a為設(shè)計(jì)的擾動(dòng)的數(shù)值分布，圖10b為相應(yīng)位置全采樣方法得到的數(shù)值分布，圖10c為局部采樣數(shù)值分布.明顯看出圖10c的數(shù)值分布與圖10a更為接近.經(jīng)過(guò)計(jì)算圖10b的均值為1.8×10-2m·ns-1，圖10c的均值為3.1×10-2m·ns-1，更接近擾動(dòng)位置設(shè)計(jì)的均值4.5×10-2m·ns-1，證明了局部采樣MCMC方法的有效性.

4 結(jié)論

(1)本文針對(duì)時(shí)移反演問(wèn)題，提出了一種基于MCMC方法求解的通用框架，將MCMC算法與雙差法相結(jié)合，對(duì)時(shí)移的目標(biāo)位置進(jìn)行有效的反演.MCMC反演采用擴(kuò)展Metropolis算法結(jié)合連續(xù)吉布斯采樣求解，雙差法反演時(shí)將第一次反演的結(jié)果作為初始模型，并使用局部采樣的方式，減小了計(jì)算量，提高了目標(biāo)區(qū)域的求解效率.

圖直方圖分布；(b)全采樣目標(biāo)位置直方圖分布；(c)局部采樣目標(biāo)位置直方圖分布Fig.10 (a) histogram;(b)All sampling histogram;(c)Local sampling histogram

(2)使用局部采樣的MCMC方法只對(duì)目標(biāo)區(qū)域進(jìn)行采樣，即連續(xù)吉布斯采樣法范圍限定為時(shí)移目標(biāo)區(qū)域，減小了非目標(biāo)區(qū)域的影響，使目標(biāo)區(qū)域反演結(jié)果更加準(zhǔn)確.具體在時(shí)移反演時(shí)，對(duì)比了全采樣的MCMC方法和局部采樣的MCMC方法，分析了兩者時(shí)移反演結(jié)果的誤差，證明了局部采樣的MCMC反演對(duì)目標(biāo)區(qū)域的反演誤差更小，準(zhǔn)確度更高.

(3)將MCMC方法應(yīng)用到GPR數(shù)據(jù)時(shí)移反演中，通過(guò)對(duì)GPR的模擬數(shù)據(jù)進(jìn)行反演，理論模型和反演結(jié)果基本符合，有效的反演出探地雷達(dá)數(shù)據(jù)的時(shí)移擾動(dòng)，準(zhǔn)確的反映了地下介質(zhì)不同時(shí)間內(nèi)的速度變化，說(shuō)明了MCMC反演方法的有效性和可靠性.