重力數(shù)據(jù)探測(cè)地下目標(biāo)的貝葉斯方法研究

黃佳喜 邊少鋒 紀(jì) 兵

1 海軍工程大學(xué)電氣工程學(xué)院,武漢市解放大道717號(hào),430033 2 中國(guó)地質(zhì)大學(xué)(武漢)地質(zhì)探測(cè)與評(píng)估教育部重點(diǎn)實(shí)驗(yàn)室,武漢市魯磨路388號(hào),430074

準(zhǔn)確掌握地下目標(biāo)分布信息可為合理安全規(guī)劃利用地下空間、監(jiān)測(cè)預(yù)警地質(zhì)災(zāi)害、強(qiáng)化地下國(guó)防工程探測(cè)與防護(hù)能力、消除邊境安全隱患等提供技術(shù)支撐[1-2]。確定目標(biāo)的空間位置和形狀是地下目標(biāo)探測(cè)的核心任務(wù),基于線性反演方法容易陷入局部最優(yōu)解的問題[3-5],而屬于非線性概率反演方法的貝葉斯推斷理論上能得到全局最優(yōu)解[6-8],近年來在各領(lǐng)域都得到了廣泛應(yīng)用[9]。

基于此,本文針對(duì)地下目標(biāo)形狀和位置信息均未知的情況,區(qū)分面測(cè)量和少量測(cè)線測(cè)量?jī)煞N典型應(yīng)用場(chǎng)景,通過模擬實(shí)驗(yàn)對(duì)貝葉斯推斷探測(cè)地下目標(biāo)的可行性和有效性進(jìn)行分析。

1 研究方法

1.1 規(guī)則形體重力場(chǎng)正演

空洞、隧道、地下建筑等一般具有規(guī)則形體,可用球體、長(zhǎng)方體、水平圓柱體等近似,這樣能有效簡(jiǎn)化問題,且能通過簡(jiǎn)單規(guī)則形體的組合疊加研究更復(fù)雜的地下目標(biāo)[10]。因此,本文以水平長(zhǎng)方體為研究對(duì)象,模擬分析地下目標(biāo)引起的重力和重力梯度異常。

假設(shè)在空間直角坐標(biāo)系下觀測(cè)點(diǎn)的坐標(biāo)為(xp,yp,zp),長(zhǎng)方體的范圍為x1≤x≤x2,y1≤y≤y2,z1≤z≤z2,z軸向下為正,則相應(yīng)的計(jì)算公式[11]為:

gz=GΔρ[xln(y+r)+yln(x+r)-

(1)

(2)

令ta、tb(t=x、y、z)代表積分上限和下限,則式(1)和式(2)中各參量的積分上、下限可定義為ta=t1-tp、tb=t2-tp。當(dāng)目標(biāo)體的走向長(zhǎng)度比截面的尺度和埋藏深度大得多,且觀測(cè)剖面位于目標(biāo)長(zhǎng)軸中部時(shí),可將其視為二度體,無需考慮長(zhǎng)軸方向的重力場(chǎng)變化,此時(shí)可對(duì)式(1)和式(2)進(jìn)行進(jìn)一步簡(jiǎn)化[10]。

如果局部測(cè)量坐標(biāo)系xoy與長(zhǎng)方體坐標(biāo)系uov不重合,即目標(biāo)體的走向與測(cè)線存在如圖1所示的銳角α(逆時(shí)針為正)時(shí),應(yīng)先進(jìn)行坐標(biāo)旋轉(zhuǎn),然后再進(jìn)行重力和重力梯度計(jì)算。水平旋轉(zhuǎn)矩陣定義如下:

(3)

圖1 測(cè)量坐標(biāo)系與長(zhǎng)方體坐標(biāo)系示意圖

(4)

1.2 貝葉斯推斷基礎(chǔ)

1.2.1 貝葉斯公式

貝葉斯公式的密度函數(shù)形式定義為[12]:

(5)

式中,π(θ)為未知參數(shù)θ的先驗(yàn)分布,即在抽取樣本Q之前對(duì)θ的認(rèn)識(shí),π(θ|q)為θ的后驗(yàn)分布,即給定樣本Q=q條件下θ的條件分布,分母稱為證據(jù)或Q的邊緣密度。

為了顯式地描述推斷過程中依賴的模型M,也可將上式寫為:

(6)

為計(jì)算出證據(jù)f(q|M),需要邊緣化(通過離散求和或積分)所有可能的f(θ|M),即根據(jù)給定模型邊緣化所有θ的先驗(yàn)[13]。

1.2.2 似然函數(shù)

測(cè)線上第i個(gè)點(diǎn)的觀測(cè)值qi可表示為:

qi=si+nb+ninst

(7)

式中,si為目標(biāo)產(chǎn)生的信號(hào),可由式(1)或式(2)得到,nb為背景場(chǎng)噪聲,ninst為測(cè)量設(shè)備產(chǎn)生的噪聲。觀測(cè)數(shù)據(jù)與正演模型之間的差異(噪聲)序列可視為服從均值為0,協(xié)方差為Σ的正態(tài)分布,其中協(xié)方差Σ=Φb+Φinst由背景場(chǎng)噪聲協(xié)方差和設(shè)備噪聲協(xié)方差組成。在此基礎(chǔ)上,將觀測(cè)序列的似然函數(shù)表示為:

(8)

有了總體信息模型和樣本數(shù)據(jù),只要獲取先驗(yàn)信息,就可根據(jù)式(8)的貝葉斯公式推斷參數(shù)的后驗(yàn)分布。

1.2.3 先驗(yàn)信息

如果僅知道地下目標(biāo)位于測(cè)量區(qū)域內(nèi),且目標(biāo)呈水平長(zhǎng)方體形狀,則長(zhǎng)方體的形狀參數(shù)(長(zhǎng)、寬、高)、位置信息(中心坐標(biāo)(x0,y0,z0)和目標(biāo)深度d、走向α及引起信號(hào)異常的剩余密度Δρ均為未知參數(shù)。在進(jìn)行貝葉斯推斷之前需要確定這些未知參數(shù)的先驗(yàn)信息,具體內(nèi)容將在實(shí)驗(yàn)部分展開介紹。至此,確定了總體模型、似然函數(shù)和先驗(yàn)信息,之后可以進(jìn)行貝葉斯推斷,推斷過程見圖2。

圖2 貝葉斯推斷基本流程

2 實(shí)驗(yàn)結(jié)果分析

2.1 基于區(qū)域測(cè)量數(shù)據(jù)的目標(biāo)探測(cè)實(shí)驗(yàn)

為驗(yàn)證本文方法的可行性,首先用1個(gè)20 m×15 m×5 m、頂部位于地面以下15 m、剩余密度為-2 300 kg/m3的長(zhǎng)方體模擬地下隧道,同時(shí)以長(zhǎng)方體為中心在60 m×60 m范圍內(nèi)分別按3 m、6 m、9 m間距進(jìn)行地面重力、垂直重力梯度測(cè)量。假設(shè)該區(qū)域內(nèi)無明顯的地質(zhì)干擾和地形起伏,對(duì)觀測(cè)數(shù)據(jù)分別加入不同程度的高斯白噪聲,然后驗(yàn)證探測(cè)效果。圖3(a)為加入6 μGal白噪聲的重力異常分布,圖3(b)為加入6 E白噪聲的垂直重力梯度異常分布。

圖3 模擬的區(qū)域重力異常和垂直重力梯度異常

根據(jù)上節(jié)內(nèi)容,進(jìn)行貝葉斯推斷前需要確定所有未知參數(shù)的先驗(yàn)信息。考慮到長(zhǎng)方體的長(zhǎng)(l)、寬(w)、高(h)和埋深(d)必須為正值,因此可假設(shè)其先驗(yàn)服從伽馬分布;長(zhǎng)方體中心的水平坐標(biāo)(x0,y0)不可能位于測(cè)區(qū)以外且取值可正可負(fù),因此設(shè)定其服從正態(tài)分布;地下空洞一般因質(zhì)量虧空才產(chǎn)生異常信號(hào),因此剩余密度的先驗(yàn)可設(shè)定為服從均值為負(fù)、具有一定標(biāo)準(zhǔn)差的正態(tài)分布。表1詳細(xì)列舉了各參數(shù)的定義及先驗(yàn)信息,相應(yīng)的概率密度分布見圖4。

表1 未知參數(shù)的定義及先驗(yàn)分布

圖4 未知參數(shù)的先驗(yàn)概率密度分布

結(jié)合表1的先驗(yàn)信息及圖3的觀測(cè)數(shù)據(jù),在Python中調(diào)用PYMC3庫(kù)函數(shù)即可實(shí)現(xiàn)貝葉斯推斷過程。對(duì)采樣過程進(jìn)行收斂性和相關(guān)性分析后,最終得到如表2所示的后驗(yàn)統(tǒng)計(jì)信息。

以長(zhǎng)度參數(shù)為例,從表2中可以看出,即使以3 m間距、噪聲標(biāo)準(zhǔn)差僅為2 μGal的精度進(jìn)行微重力測(cè)量,反演結(jié)果的94% HDI(相應(yīng)的上下界分別對(duì)應(yīng)表中3% HDI和97% HDI)位于均值為4.1的概率區(qū)間(0.1,9.4)內(nèi),概率取值區(qū)間不僅分散,且與真值的偏差很大。如果使用9 m間距、噪聲標(biāo)準(zhǔn)差為6 E的垂直重力梯度測(cè)量,其反演結(jié)果的94%概率區(qū)間為(18.5,20.9),94% HDI的均值為19.7 m,概率取值區(qū)間不僅相對(duì)集中,結(jié)果也更接近真值。結(jié)合其他參數(shù)的反演結(jié)果可以看出,微重力反演淺源目標(biāo)的效果遠(yuǎn)不如垂直重力梯度。這是因?yàn)樵撃繕?biāo)引起的重力異常信號(hào)最大量級(jí)也只有百μGal級(jí),接近商用重力儀的精度極限,很容易受外界噪聲的干擾;而相應(yīng)的地面垂直重力梯度異常信號(hào)接近100 E,即使加入標(biāo)準(zhǔn)差為6 E的噪聲,仍然具有很高的信噪比。另外,從表2的反演結(jié)果可以看出,減小測(cè)點(diǎn)間距進(jìn)行加密測(cè)量并沒有使反演結(jié)果得到明顯改善,但測(cè)點(diǎn)間距過大時(shí),如采用12 m的測(cè)量間距,就會(huì)出現(xiàn)反演結(jié)果發(fā)散的現(xiàn)象(文中未列出實(shí)驗(yàn)結(jié)果)。

總體來看,水平位置、剩余密度及長(zhǎng)、寬等參數(shù)的反演效果要明顯優(yōu)于深度、高度等垂向參數(shù)。圖5以表2中9 m測(cè)量間距、6 E噪聲的結(jié)果為例直觀地反映了這種差異,其中水平位置、長(zhǎng)、寬等參數(shù)的75% HDI與真值相差不足1 m;而深度、高度等參數(shù)的75% HDI與真值最大相差5 m,這是因?yàn)橹亓?重力梯度數(shù)據(jù)本身就沒有垂向分辨能力,不管真實(shí)深度如何,反演結(jié)果總是趨于近地表。熊盛青[14]將這一現(xiàn)象稱為“趨膚效應(yīng)”或模型“上漂”,下一步將嘗試結(jié)合深度加權(quán)技術(shù)克服該現(xiàn)象。

圖5 模型參數(shù)的貝葉斯后驗(yàn)概率分布

為驗(yàn)證先驗(yàn)信息對(duì)推斷結(jié)果的影響,將表1中γ分布函數(shù)的參數(shù)調(diào)整為(3,1)或(4,0.5),并重復(fù)上述實(shí)驗(yàn)。反演結(jié)果與表2無明顯差異,說明先驗(yàn)信息取值范圍具有合理性。因此,貝葉斯推斷能合理估計(jì)探測(cè)目標(biāo)中的未知參數(shù),地面重力梯度觀測(cè)比重力異常探測(cè)更具探測(cè)優(yōu)勢(shì)。

2.2 基于少量測(cè)線重力梯度數(shù)據(jù)的目標(biāo)探測(cè)實(shí)驗(yàn)

在一些特殊應(yīng)用領(lǐng)域,受地形地物特征或觀測(cè)條件的限制,往往無法進(jìn)行大面積區(qū)域測(cè)量,甚至無法進(jìn)行地面測(cè)量。因此,有必要研究低空平臺(tái)少量測(cè)線數(shù)據(jù)探測(cè)地下目標(biāo)的可行性。針對(duì)這一應(yīng)用場(chǎng)景,顧及上一實(shí)驗(yàn)中垂直重力梯度探測(cè)效果明顯優(yōu)于微重力的結(jié)果,本文僅利用重力梯度數(shù)據(jù)分析貝葉斯推斷的應(yīng)用效果。

Romaides等[15]曾對(duì)一處地下地鐵停車場(chǎng)進(jìn)行地面微重力測(cè)量,該地下設(shè)施長(zhǎng)150 m、寬18.3 m、高8.5 m,頂部距地表3.7 m。因剩余密度未知,本文取值為-1 900 kg/m3,并在目標(biāo)中部生成2條與目標(biāo)體長(zhǎng)軸呈35°方向的測(cè)線(圖1)。測(cè)線長(zhǎng)90 m,以長(zhǎng)方體長(zhǎng)軸為中心,測(cè)點(diǎn)間距分別設(shè)為2 m、4 m、8 m和12 m,并加入不同程度的白噪聲(表3)。圖6(a)為地面(虛線)及30 m高度上(實(shí)線)重力和垂直重力梯度正演結(jié)果。盡管地面重力異常(藍(lán)色虛線)量級(jí)明顯低于文獻(xiàn)[15]中的實(shí)測(cè)值,但較小的剩余密度更有助于驗(yàn)證本文方法的可行性。圖6(b)為30 m高度上的重力梯度正演結(jié)果(藍(lán)色曲線)及加入6 E白噪聲的結(jié)果(紅色曲線)。

表3 少量測(cè)線垂直重力梯度數(shù)據(jù)的貝葉斯推斷結(jié)果

圖6 不同高度上的重力異常和重力梯度測(cè)線數(shù)據(jù)

若已知該目標(biāo)呈長(zhǎng)方體狀,且可視為水平無限長(zhǎng)二度體,那么在目標(biāo)探測(cè)時(shí)就無需考慮長(zhǎng)軸方向的中心坐標(biāo)及長(zhǎng)度,其余未知參數(shù)的先驗(yàn)信息同表2。此外,因目標(biāo)體的走向未知,需要增加水平走向參數(shù)α,根據(jù)定義,其取值在(-90°～90°),可將其先驗(yàn)分布設(shè)定為均勻分布U(-90°,90°)。

推斷結(jié)果如表3所示,可以看出,只有1條測(cè)線時(shí),盡管噪聲的標(biāo)準(zhǔn)差達(dá)到正演信號(hào)最大值的20%左右,但除水平走向α之外,其余參數(shù)的真值基本都位于后驗(yàn)的94% HDI內(nèi),但參數(shù)概率取值區(qū)間較上節(jié)較為分散,結(jié)果不一定真實(shí)可靠,加密測(cè)量也未能改善推斷結(jié)果。雖然不能從1條測(cè)線的數(shù)據(jù)推斷出水平走向角,但不難發(fā)現(xiàn),94% HDI的上下界(分別對(duì)應(yīng)97% HDI和3% HDI)大致確定了目標(biāo)走向的可能取值β0。因此,根據(jù)1條測(cè)線的推斷結(jié)果,基本能夠確定走向角介于±β0之間。

為推斷走向角,使用2條測(cè)線的數(shù)據(jù)重復(fù)上一實(shí)驗(yàn)。從表3中可以看出,2條測(cè)線的觀測(cè)數(shù)據(jù)能夠很好地推斷出包括走向角在內(nèi)的所有參數(shù)信息。若進(jìn)一步將表3中剩余密度的先驗(yàn)分布均值設(shè)置為-1 600 kg/m3或-2 300 kg/m3,推斷結(jié)果并未產(chǎn)生明顯的變化。由此說明,在信噪比較大的情況下,即使地下目標(biāo)的剩余密度未知,也可在反演中作出相對(duì)合理的估計(jì)。

以表3中2條測(cè)線、間距4 m、噪聲標(biāo)準(zhǔn)差為3 E的后驗(yàn)結(jié)果繪制出參數(shù)的邊緣分布及y0、z0的聯(lián)合概率密度分布(圖7)。結(jié)合表3中94% HDI結(jié)果可以看出,圖7中各參數(shù)的75% HDI和概率區(qū)間均值的可靠性更高。因此,在具體應(yīng)用中可根據(jù)需要設(shè)定不同的HDI參考值。

圖7 模型參數(shù)的貝葉斯后驗(yàn)概率分布

以上實(shí)驗(yàn)結(jié)果表明,即使只有1條與目標(biāo)相交的測(cè)線,也能對(duì)地下目標(biāo)的位置和形狀參數(shù)作出相對(duì)保守、合理的估計(jì);當(dāng)測(cè)線與目標(biāo)走向不垂直時(shí),至少需要2條與目標(biāo)相交的測(cè)線才能確定地下目標(biāo)的走向。

3 結(jié) 語(yǔ)

本文通過正演分析計(jì)算規(guī)則形狀的地下目標(biāo)產(chǎn)生的重力/重力梯度異常信號(hào),進(jìn)而驗(yàn)證貝葉斯推斷用于地下目標(biāo)探測(cè)的可行性和有效性。實(shí)驗(yàn)結(jié)果表明:

1)地面垂直重力梯度測(cè)量較微重力測(cè)量具有更高的信噪比,使用貝葉斯推斷能對(duì)地下目標(biāo)的位置、形狀、剩余密度等參數(shù)作出合理估計(jì)。結(jié)合測(cè)區(qū)的坐標(biāo)范圍合理設(shè)定目標(biāo)位置、形狀等參數(shù)的先驗(yàn)信息,可使反演結(jié)果更趨于真值。

2)對(duì)于未知走向的地下目標(biāo),至少需要2條與目標(biāo)相交的測(cè)線才能確定目標(biāo)的走向參數(shù)。

本文研究初步驗(yàn)證了貝葉斯推斷用于地下目標(biāo)探測(cè)的可行性,但對(duì)于實(shí)驗(yàn)中出現(xiàn)的深度參數(shù)反演結(jié)果“上漂”的問題尚未解決。盡管通過模擬實(shí)驗(yàn)分析了地下目標(biāo)在一定高度上的重力場(chǎng)信號(hào)及其貝葉斯推斷結(jié)果,但針對(duì)復(fù)雜背景場(chǎng)相關(guān)噪聲及地形起伏影響下地下目標(biāo)信號(hào)的探測(cè)識(shí)別仍有待進(jìn)一步研究。