基于局部采樣MCMC方法的時移探地雷達反演

王升超，韓立國，鞏向博，張盼

吉林大學地球探測科學與技術學院，長春 130026

0 引言

馬爾科夫鏈蒙特卡洛方法(MCMC)是基于數學概率分布模擬的一種方法，目前已經廣泛應用于理論物理、信號通信、醫學等領域，在地球物理領域的應用也在不斷的發展中，特別是應用于地球物理反演問題(Grandis et al.,1999;Malinverno and Leaney,2000).相比于其他反演方法，MCMC方法的優點在于全局尋優，不會陷入局部極小，并且不依賴于準確的先驗模型，可以引入更為復雜的先驗信息.

對于地球物理反演問題，基于最小二乘的方法得到了廣泛的應用，Menke(1989)提出了確定性最小二乘法，而Tarantola和Valette(1982)提出了最小二乘的概率反演方法，在貝葉斯框架下，反演問題的解是一個概率密度函數，稱之為后驗概率密度函數.Hansen等(2006)提出了一種基于連續模擬的線性走時層析成像的概率反演方法，他們利用了經典最小二乘反演和克里格法(Journel and Huijbregts，1978)的等價性.Nielsen等(2010)給出了這種方法在井間探地雷達數據中的應用.Gloaguen等(2005a，b)提出了一種基于克里格法的誤差模擬的相關方法，相當于概率最小二乘法，并將其應用于井間探地雷達層析成像.Giroux和Gloaguen(2012)將這種方法用于各向異性速度場的反演.這些方法僅對嚴格線性反演問題有效，并且依賴于描述噪聲模型和先驗模型的高斯統計的固有假設，先驗模型必須以高斯形式給出，先驗模型由均值和協方差模型定義.國內，張廣智等(2011)將MCMC方法應用于疊前地震數據反演，在波阻抗反演方面取得了一定效果.殷長春等(2014)將概率反演的方法應用到航空電磁領域，克服了初始模型的影響，成功反演了深度低阻層.王朋巖等(2015)成功使用MCMC方法反演了巖性參數.孫月成(2018)則將基于MCMC 算法的地質統計學反演應用在油藏模擬當中，取得了良好的效果.

但在實際反演問題當中，先驗信息往往比高斯模型描述的更為復雜，而且先驗信息也不能用準確的數學公式表達.Hansen等(2008b)展示了拓展的Metropolis算法(Mosegaard et al.，1995)在非線性井間層析成像問題中的應用，其中先驗模型是非高斯的，并可以由任何地質統計學方法定義.擴展的Metropolis算法不需要先驗信息的準確數學表達式，可以用來采樣后驗概率密度函數，包括高度非線性的反演問題.在拓展的Metropolis算法中，只要有一個能夠對先驗概率密度函數進行采樣的“黑箱”算法，就可以完成對先驗信息的采樣任務.Hansen等(2008a,2012)提出了連續吉布斯采樣的方法，連續吉布斯算法在采樣先驗信息時，可以有效的控制擾動的步長，具有很高的采樣效率(Gómez-Hernández and Journel，1993).這樣，將連續吉布斯采樣作為一種“黑箱”算法應用于拓展的Metropolis算法中，使得在進行概率反演時，對先驗模型的選擇變得非常靈活(Journel and Zhang，2006)，可以對任何地質統計學算法定義的先驗信息進行采樣，既可用于相對簡單的兩點統計模型(例如基于高斯的先驗模型)，也可以用于更復雜的多點的統計模型，通過多點算法可以靈活地模擬高熵、低熵結構或其組合的實現，使得在求解概率反演時，可以引入更復雜的先驗模型信息(Hansen et al.,2013).

后驗概率分布的求解是通過模型提供的先驗信息，數據的不確定性和模型參數等綜合求解.其中所有可用的先驗信息都由先驗概率密度函數描述.多種信息綜合反映了后驗樣本模型的變化，為后驗采樣提供了依據.概率反演方法除了能提供后驗模型參數的協方差外，也可以解決復雜的問題，如地質連通性的概率或流體的停留時間(Mosegaard,1998).

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

時移反演有兩種反演策略，稱為連續反演法和雙差法(Waldhauser and Ellsworth，2000)，Huang等(2020)在雙差法的基礎上實現了基于目標導向的全波形時移反演，提高了時移反演的計算效率.兩種方法的主要區別在于雙差法減小了非目標區的影響，提高了目標區的反演精度.在時移反演中，作者使用雙差法結合MCMC方法對探地雷數據進行反演.對于時移探地雷達數據，要解決目標位置的變化，需要對兩組數據在不同觀測點進行至少兩次反演計算.兩次反演計算量大，非目標區域的反演也會影響目標區域的精度.為了解決這個問題，作者采用了局部采樣的MCMC反演方法，在雙差法的第二次反演過程中，采用局部采樣即只對目標區域采樣的方法代替了全局采樣法.這樣，減少了非采樣區域對時移位置的影響，在提高計算效率的同時，增加了目標區域的反演精度.

本文中，作者演示了如何通過使用拓展的Metropolis算法結合使用連續吉布斯采樣的地質統計學算法定義的先驗信息對時移的探底雷達數據實現MCMC的反演.將MCMC方法應用于探地雷達時移反演中，首先利用連續吉布斯采樣和Metropolis算法對反演問題進行求解.然后，結合雙差法得到目標區域的變化.作者將該方法應用到GPR模擬數據中，測試該方法的效果，分析對比了局部采樣的MCMC反演和全局采樣的MCMC的反演誤差，證明了局部采樣MCMC反演的有效性.

1 MCMC反演

1.1 概率反演公式

對于地球物理反演問題，地下構造可以用一組模型參數m來表示，觀測數據可以用一組數據d來表示.正演問題是指通過一個函數映射f來獲取觀測數據d(Tarantola and Valette，1982)：

d=f(m)，

(1)

函數f通常基于相應的物理關系，對應的反演問題的公式可以表示為

m=f-1(d)，

(2)

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

反演問題是通過先驗信息來獲取模型參數，基于概率方法的反演公式可以表示為

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

(3)

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

(4)

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

1.2 基于射線的正演方法

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

(5)

快速推進法是一種具有高精度、高效率、無條件穩定等特性的走時計算方法(孫章慶，2008)，利用多階快速推進法可以有效的求解程函方程.這種正演模型是非線性的，因為程函方程對應于波動方程的高頻近似，通常被稱為高頻射線近似.信號源和接收器之間的走時d可由式(6)給出：

(6)

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

1.3 拓展的Metropolis 算法

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

(1)提出了一個候選模型mpro，它是在當前模型mcur中給出一個擾動后產生的新模型，同時也是先驗概率密度函數的一次實現.

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

(7)

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

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

(2)計算當前模型和候選模型的似然函數，L(mcur)和L(mpro).

(4)如果建議的模型被接受，那么用mpro代替當前模型mcur的轉換被接受，即mcur=mpro，否則，仍然在mcur位置，再次產生一個隨機擾動進行下一次迭代.

2 雙差法時移反演

本文主要采用雙差法來解決時移反演問題，在雙差法的第二次反演時，采用局部采樣的MCMC方法進行反演，該雙差法的主要流程如圖1所示.

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

3 合成的GPR數據測試

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

圖2 數據觀測系統分布，共702組發射源和接收點組合Fig.2 Recording geometry,ND=702 pairs of sources (red crosses)and receivers (black dots)are represented by a connecting black line

圖3 起始速度模型mT1，監測速度模型mT2和速度擾動Fig.3 Initial model mT1,monitor model mT2 and the perturbation

使用速度模型mT1和mT2正演模擬獲得相應的初至時間數據data1 和 data2作為觀測數據，圖4顯示了相應的正演數據分布，其中黑色曲線代表data1，紅色曲線代表data2.由圖4可以看出，兩組正演數據的變化趨勢大體相同，只在局部位置有差異，可以認為是速度擾動區域在正演數據中產生的影響.

圖4 正演模擬的初至走時數據Fig.4 The simulated traveltime data1 and data2

在反演之前，首先要根據先驗信息對先驗概率密度函數進行采樣，先驗采樣是一個隨機過程，任何滿足先驗條件的模型都有可能出現.圖5顯示了對先驗概率密度函數進行連續吉布斯采樣后，得到的五個獨立實現的先驗速度模型，由圖5可以看出，每個先驗模型的速度分布各不相同.

圖5 先驗概率密度函數的五個獨立隨機采樣結果Fig.5 Five statistically independent realizations of the priori probability density

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

圖6 后驗概率密度函數的五個采樣樣本Fig.6 Five statistically independent realizations of the posterior probability density

為進一步檢驗反演的準確性，我們對這些后驗采樣得到的模型進行正演模擬，正演的走時數據如圖7所示，其中紅色曲線為觀測數據data1，黑色曲線為多個后驗采樣模型的正演數據.可以發現，通過拓展的Metropolis算法獲得的后驗樣本，其正演的結果與觀測數據data1有很高的一致性，證明了反演的方法的有效性.

圖7 后驗樣本正演得到的走時數據Fig.7 The traveltime data of data and posterior

雙差法時移反演中，第二次反演求解速度模型mT2時，需要選擇第一次反演的結果作為初始模型，我們選擇后驗樣本的一個結果作為起始模型進行第二次MCMC反演.不同于第一次反演時的全采樣方法，第二次反演時我們只對目標區域使用連續的吉布斯算法進行采樣.為對比局部采樣的反演效果，我們同時設置了一組模型全采樣的數據進行對比，具體結果如圖8所示.

圖8a顯示了使用全采樣的擴展Metropolis算法采樣后驗概率密度函數得到的五個速度模型，用圖8a的后驗模型與第一次反演得到的初始模型相減，得到擾動變化如圖8b所示.圖8c是用局部采樣的MCMC反演得到的后驗分布，圖8d是用圖8c與初始反演模型相減的結果.對比圖8a、c和模型mT2，可以看出圖8a、c在淺層位置都反演出了一定的高速分布，但是圖8a淺層高速分布的位置范圍變化較大，高速位置分布隨機并不集中在目標區域，圖8c則主要集中在目標區域位置.用圖8b、d對比更加明顯，全采樣的反演相對初始模型的變化位置，除淺層外，深層也有變化.

圖8 (a)T2時刻全采樣反演結果；(b)全采樣擾動變化；(c)T2時刻局部采樣反演結果；(d)局部采樣擾動變化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)全采樣法的目標區域反演結果；(b)局部采樣法的目標區域反演結果Fig.9 (a)Target area of all sampling method;(b)Target area of local sampling method

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

4 結論

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

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

(2)使用局部采樣的MCMC方法只對目標區域進行采樣，即連續吉布斯采樣法范圍限定為時移目標區域，減小了非目標區域的影響，使目標區域反演結果更加準確.具體在時移反演時，對比了全采樣的MCMC方法和局部采樣的MCMC方法，分析了兩者時移反演結果的誤差，證明了局部采樣的MCMC反演對目標區域的反演誤差更小，準確度更高.

(3)將MCMC方法應用到GPR數據時移反演中，通過對GPR的模擬數據進行反演，理論模型和反演結果基本符合，有效的反演出探地雷達數據的時移擾動，準確的反映了地下介質不同時間內的速度變化，說明了MCMC反演方法的有效性和可靠性.