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

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

1 海軍工程大學(xué)電氣工程學(xué)院,武漢市解放大道717號,430033 2 中國地質(zhì)大學(xué)(武漢)地質(zhì)探測與評估教育部重點實驗室,武漢市魯磨路388號,430074

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

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

1 研究方法

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

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

假設(shè)在空間直角坐標(biāo)系下觀測點的坐標(biāo)為(xp,yp,zp),長方體的范圍為x1≤x≤x2,y1≤y≤y2,z1≤z≤z2,z軸向下為正,則相應(yīng)的計算公式[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)體的走向長度比截面的尺度和埋藏深度大得多,且觀測剖面位于目標(biāo)長軸中部時,可將其視為二度體,無需考慮長軸方向的重力場變化,此時可對式(1)和式(2)進(jìn)行進(jìn)一步簡化[10]。

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

(3)

圖1 測量坐標(biāo)系與長方體坐標(biāo)系示意圖

(4)

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

1.2.1 貝葉斯公式

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

(5)

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

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

(6)

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

1.2.2 似然函數(shù)

測線上第i個點的觀測值qi可表示為:

qi=si+nb+ninst

(7)

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

(8)

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

1.2.3 先驗信息

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

圖2 貝葉斯推斷基本流程

2 實驗結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 語

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

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

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

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