侯振隆 袁 園 王恩德 付建飛 鄭玉君

(①東北大學(xué)深部金屬礦山安全開采教育部重點(diǎn)實(shí)驗(yàn)室，遼寧沈陽 110819; ②東北大學(xué)資源與土木工程學(xué)院，遼寧沈陽 110819; ③自然資源部第二海洋研究所國家海洋局海底科學(xué)重點(diǎn)實(shí)驗(yàn)室，浙江杭州 310012)

0 引言

利用重力勘探數(shù)據(jù)反演地下物性、幾何位置等參數(shù)是一種常用的數(shù)據(jù)處理解釋方法[1-2]。對于計(jì)算地質(zhì)體在水平方向上的分布范圍，邊界識別是一種快速、有效的手段。常用的邊界識別方法包括總水平導(dǎo)數(shù)法[3]、解析信號法[4]、傾斜角法及其總水平導(dǎo)數(shù)法[5-6]、Theta圖法[7]、結(jié)構(gòu)張量特征值法[8]等。重力梯度數(shù)據(jù)與重力異常相比具有更高的信噪比，使用梯度數(shù)據(jù)識別邊界效果較好[9]。

目前，邊界識別方法還存在一些不足之處，如對深部地質(zhì)體邊界識別不清晰，結(jié)果中存在虛假邊界等。因此，學(xué)者們提出了不同的改進(jìn)方式，例如：方向總水平導(dǎo)數(shù)(Directional Total Horizontal Derivative,DTHD)法[10-11]，改進(jìn)結(jié)構(gòu)張量法[12]，均衡總水平導(dǎo)數(shù)法與正則化解析信號法[13]，全張量梯度數(shù)據(jù)的特征值法[14]，水平方向Theta圖法[15]，高階導(dǎo)數(shù)法[16]，全張量梯度角度法[17]，偽總梯度余弦值法[18]等，這些研究在一定程度上改善了以上問題。其中，DTHD法比解析信號法具有更高的分辨率[19]。DTHD法可通過計(jì)算x、y等方向上的總水平導(dǎo)數(shù)，進(jìn)一步提高檢測能力，在實(shí)際應(yīng)用中能夠從這兩個(gè)方向上較好地識別地質(zhì)體邊界。

相關(guān)系數(shù)在位場數(shù)據(jù)的處理解釋中可用于快速成像[20-21]和反演[22-24]。相關(guān)系數(shù)不僅能夠表示觀測異常與地下空間幾何函數(shù)的相關(guān)程度，還可以表示不同物理量的相似程度。根據(jù)這一性質(zhì)，馬國慶等[25]通過總水平導(dǎo)數(shù)與垂直導(dǎo)數(shù)的相關(guān)系數(shù)識別邊界，該方法聯(lián)合兩種物理量增強(qiáng)了識別效果；徐夢龍等[26]還利用各方向均方差的相關(guān)系數(shù)實(shí)現(xiàn)了邊界識別，通過聯(lián)合多個(gè)方向窗口中的均方差數(shù)據(jù)增強(qiáng)了相關(guān)性特征。

本文利用重力梯度數(shù)據(jù)計(jì)算x和y方向總水平導(dǎo)數(shù)的相關(guān)系數(shù)識別邊界，結(jié)合多方向窗口技術(shù)改善識別效果，并通過改進(jìn)的二值化和閾值處理(Binarization and Threshold Processing,BTP)方法[27]突出解釋結(jié)果中的邊界位置，提高識別精度。

1 多方向窗口聯(lián)合相關(guān)系數(shù)方法原理

1.1 x和y方向總水平導(dǎo)數(shù)的相關(guān)系數(shù)

根據(jù)文獻(xiàn)[10]，x和y方向的總水平導(dǎo)數(shù)分別定義為

(1)

(2)

式中：V表示重力場；Vxy、Vxz、Vyx與Vyz表示重力場梯度，且Vxy=Vyx。THDx與THDy的相關(guān)系數(shù)定義為

(3)

式中： THDx,i和THDy,i分別表示窗口內(nèi)第i個(gè)THDx和THDy數(shù)據(jù)；N表示窗口內(nèi)的總點(diǎn)數(shù)。

式(3)中，-1≤C≤1，屬無量綱的物理量。C的絕對值越接近于1，表示THDx與THDy的相關(guān)程度越高；越接近于0，表示相關(guān)程度越低。C值為正數(shù)表示正相關(guān)；若為負(fù)數(shù)，表示負(fù)相關(guān)；當(dāng)C等于0，表示不相關(guān)。根據(jù)式(1)和式(2)，地質(zhì)體在x、y方向的邊界可分別根據(jù)THDx和THDy的極大值分布特征進(jìn)行解釋。C越接近于-1，即THDx與THDy呈負(fù)相關(guān)時(shí)，該點(diǎn)為地質(zhì)體邊界的可能性越大。這種利用極小值分布特征解釋地質(zhì)體邊界的方法即是相關(guān)系數(shù)邊界識別(Correlation Coefficients Edge Detection,CCED)法。

1.2 多方向窗口的聯(lián)合方法

式(3)的計(jì)算使用了滑動窗口，通常窗口的中心點(diǎn)就是測點(diǎn)位置。徐夢龍等[26]為獲得測點(diǎn)與其鄰近各方向上的數(shù)據(jù)在波動性方面的相關(guān)性，使用了8個(gè)方向的窗口，即測點(diǎn)分別在窗口的右下角、下方中心、左下角、左側(cè)中心、左上角、上方中心、右上角、右側(cè)中心，并分別定義測點(diǎn)位于中心和邊緣的8個(gè)窗口為窗口0和A～H(圖1)。

對于任意測點(diǎn)(x,y)，使用窗口0、A～H計(jì)算相關(guān)系數(shù)能夠得到9個(gè)Ck(k=0，A,B,C,D,E,F,G,H)，這些值共同表征以該測點(diǎn)為中心的區(qū)域上THDx與THDy的相關(guān)度。利用不同觀測場或不同計(jì)算條件下的相關(guān)系數(shù)可以計(jì)算聯(lián)合相關(guān)系數(shù)[28]，即多方向窗口聯(lián)合相關(guān)系數(shù)(Multi-directional Window Joint Correlation-coefficient,MWJC)

Cjoint=C0×CA×CB×CC×CD×CE×

CF×CG×CH

(4)

Cjoint可以增強(qiáng)THDx與THDy的相關(guān)性特征，進(jìn)一步提高邊界識別效果。

1.3 改進(jìn)的二值化與閾值處理方法

利用式(4)計(jì)算的Cjoint結(jié)果中可能存在干擾，影響地質(zhì)體邊界的識別和成圖效果。根據(jù)數(shù)字圖像處理中的BTP原理，設(shè)定一閾值α，并據(jù)其對數(shù)據(jù)劃分區(qū)間，對不同區(qū)間內(nèi)的數(shù)據(jù)賦值為0或1。這樣整個(gè)數(shù)據(jù)被轉(zhuǎn)化為只包含0或1的數(shù)據(jù)集，圖像中的細(xì)節(jié)被弱化，邊界將更清晰地顯示出來。MWJC法是利用THDx與THDy的負(fù)相關(guān)性判斷邊界的，即相關(guān)系數(shù)的分布區(qū)間為[-1,1]。為了不破壞結(jié)果的幅值分布規(guī)律并保留邊界細(xì)節(jié)，本文提出改進(jìn)的BTP(Improved BTP,IBTP)方法。具體步驟包括：

圖1 多方向窗口示意圖

(1)將Cjoint的取值區(qū)間[-1,1]劃分為若干個(gè)長度相等的子區(qū)間，并統(tǒng)計(jì)各子區(qū)間上的數(shù)量；

(2)設(shè)定閾值α等于某個(gè)子區(qū)間的最小值，且α≤0(大于α的所有子區(qū)間中包含的結(jié)果數(shù)量占較大比例)；

(3)將大于α的子區(qū)間標(biāo)注為1。這樣既保留了具有負(fù)相關(guān)性特征的幅值的細(xì)節(jié)，也將非邊界結(jié)果的正值全部轉(zhuǎn)化為1，一定程度上增加了圖像中邊界顯示的對比度。

2 理論模型試驗(yàn)

為了驗(yàn)證本文提出的MWJC法在邊界識別中的效果，設(shè)計(jì)了兩個(gè)理論模型進(jìn)行計(jì)算分析。模型的測點(diǎn)均勻分布在z=0平面上，x、y方向觀測范圍均為-1000～1000m。

2.1 模型一

模型一包括兩個(gè)大小相同的正方體，平面位置見圖2。異常體1和異常體2的邊長均為400m，埋深分別為200m和300m，剩余密度均為1.0g/cm3。觀測點(diǎn)距、線距均為10m，總計(jì)201×201個(gè)測點(diǎn)。

該模型正演模擬的Vxy、Vxz和Vyz如圖2所示，進(jìn)而計(jì)算得到THDx與THDy(圖3)及不同窗口的CCED法識別結(jié)果(圖4)和MWJC分布(圖5)。這里CCED的計(jì)算窗口為5×5，即包含5×5個(gè)觀測數(shù)據(jù)，尺寸為40m×40m。由圖3可見，THDx和THDy能夠分別檢測到兩個(gè)異常體的x和y方向的邊界，但檢測不到另一個(gè)方向上的邊界，且強(qiáng)異常所覆蓋范圍較大。由圖4可見，使用窗口0的CCED識別結(jié)果(圖4e)兼有THDx與THDy的優(yōu)點(diǎn)，即均可很好地識別兩個(gè)異常體的x和y方向的邊界，其中，淺部異常體1的邊界最準(zhǔn)確，檢測的深部異常體2的范圍略大于實(shí)際范圍。MWJC計(jì)算結(jié)果(圖5)進(jìn)一步增強(qiáng)了異常體邊界位置的負(fù)相關(guān)性，雖檢測到的深部異常體2的邊界雖仍然大于實(shí)際范圍，但較圖4e的檢測結(jié)果更接近于實(shí)際情況。需要注意的是，圖5中識別的深部地質(zhì)體的范圍比實(shí)際模型范圍略大，且形態(tài)也略顯發(fā)散；相對而言，對淺部異常體1的邊界識別更準(zhǔn)確?；贑CED及MWJC的識別結(jié)果均在一定程度上改善了由于異常體埋深增加而導(dǎo)致的邊界識別分辨率下降問題； 甚至從圖4和圖5還可見，深部地質(zhì)體邊界的負(fù)相關(guān)程度高于淺部地質(zhì)體。

同一測點(diǎn)在不同窗口中的相對位置會不同，因此基于窗口A～H的CCED識別結(jié)果也有所不同：識別的地質(zhì)體邊界會隨窗口方向的變化而移動。例如，在圖1中位于窗口A中的右下角的測點(diǎn)，在圖3a中則位于窗口A的右上角，即識別出的地質(zhì)體邊界向圖中實(shí)際邊界的右上方移動。基于MWJC的識別結(jié)果進(jìn)一步突出了邊界中四個(gè)頂點(diǎn)的位置，提高了識別的邊界與背景的對比度。根據(jù)統(tǒng)計(jì)結(jié)果確定α=-0.2，經(jīng)過IBTP處理的MWJC結(jié)果(圖6)完全消除了圖5a中的干擾，邊界顯示更清晰。從圖6可見，每個(gè)正方體的四個(gè)頂點(diǎn)位置均被準(zhǔn)確地識別，結(jié)合圖4e，能夠判斷不同深度正方體的邊界。

圖2 模型一重力梯度數(shù)據(jù)

圖3 模型一的THDx(a)和THDy(b)分布

圖4 模型一的Ck分布

圖5 模型一的Cjoint分布(a)及Cjoint子區(qū)間數(shù)量統(tǒng)計(jì)柱狀圖(b)

圖6 圖5a經(jīng)IBTP處理的結(jié)果

由以上分析可見，基于MWJC的識別結(jié)果能夠綜合THDx與THDy在y和x方向的邊界識別能力。與其他方法[8,12,16]相比，無需考慮矩陣的秩、歸一化等問題，原理和計(jì)算方法都很簡單，易于編程實(shí)現(xiàn)； 通過IBTP處理后，MWJC的分布特征能夠很好地突出地質(zhì)體邊界。

2.2 模型二

該模型包含三個(gè)剩余密度均為1.0g/cm3的立方異常體，埋深均為200m，大小分別為1000m×1000m×5000m、300m×300m×300m和300m×300m×300m，在水平方向上三個(gè)異常體組呈“凹”形分布，相對位置見圖7。觀測點(diǎn)距、線距均為10m，共201×201個(gè)測點(diǎn)。在理論數(shù)據(jù)中加入5%的高斯噪聲。計(jì)算Cjoint的窗口大小為3×3， 即尺寸為20m×20m， 并采取高斯濾波的方法壓制噪聲。

由于模型中的異常體2和異常體3的體積較小，所以不論是從觀測數(shù)據(jù)(圖7)還是DTHD分布(圖8)都不易分辨；但是從圖9所示的C分布特征能夠識別到異常體邊界，圖中的極小值揭示了異常體的邊界位置。

根據(jù)Cjoint分布(圖10)可以清晰地識別出三個(gè)異常體的邊界拐點(diǎn)，與模型相符。通過統(tǒng)計(jì)Cjoint在各個(gè)子區(qū)間中的個(gè)數(shù)，設(shè)閾值α=-0.6對圖10a進(jìn)行IBTP處理，結(jié)果如圖11所示。由圖可見正值干擾被去除，結(jié)合模型一的計(jì)算結(jié)果，發(fā)現(xiàn)MWJC法對模型邊界中拐角的識別結(jié)果較好，但對于模型邊界中的銳角，只能識別分布趨勢。即便如此，仍可認(rèn)為該結(jié)果對于邊界識別具有實(shí)際意義。

圖7 模型二重力梯度數(shù)據(jù)(含5%噪聲)

圖8 模型二THDx(a)和THDy(b)分布

圖9 模型二的Ck分布

圖10 模型二的Cjoint分布(a)及Cjoint子區(qū)間數(shù)量統(tǒng)計(jì)柱狀圖(b)

圖11 圖10a經(jīng)IBTP處理的結(jié)果

2.3 影響識別結(jié)果的因素分析

為了進(jìn)一步檢驗(yàn)方法的實(shí)用性，以模型二為例，討論算法中影響邊界識別結(jié)果的因素。

由于算法采用的是正方形窗，窗口邊長為(n-1)×Δx，其中Δx表示點(diǎn)距或線距，n≥3且為奇數(shù)。當(dāng)Δx為常數(shù)時(shí)，調(diào)節(jié)窗口中的測點(diǎn)數(shù)即可改變窗口尺寸。將窗口增大至11×11，即尺寸為100m×100m。參數(shù)C及Cjoint和經(jīng)IBTP處理后的結(jié)果分別如圖12和圖13所示。與圖9～圖11對比可見，窗口尺寸增大時(shí)，CCED的識別效果無明顯改善；經(jīng)聯(lián)合多窗口后，MWJC識別的邊界信息減少； 選取α=-0.2進(jìn)行IBTP處理，同樣不能識別出較準(zhǔn)確的邊界。這是由于增大窗口尺寸會使包含的測點(diǎn)數(shù)增大，改變方向時(shí)，更多對應(yīng)方向的觀測點(diǎn)匯入窗口，計(jì)算結(jié)果會有明顯的方向性。聯(lián)合多窗口的CCED結(jié)果降低了結(jié)果中的負(fù)相關(guān)信息量，故出現(xiàn)了部分邊界缺失等問題。

繼續(xù)增加點(diǎn)距和線距至50m，其他條件不變，測點(diǎn)數(shù)減至41×41，處理結(jié)果如圖14和圖15所示。與圖9～圖11對比發(fā)現(xiàn)，盡管增加窗口大小至11×11降低了識別精度，但CCED結(jié)果仍能識別出模型在水平方向上的分布范圍和形狀；不同方向窗口CCED結(jié)果的方向性進(jìn)一步增強(qiáng)，即結(jié)果中負(fù)值區(qū)域與真實(shí)位置在對應(yīng)方向上的偏差更大，無法識別異常體的部分邊界，即MWJC結(jié)果也無法增強(qiáng)負(fù)相關(guān)性，說明此時(shí)算法失效。

上述分析說明：算法選用的窗口尺寸與觀測數(shù)據(jù)的點(diǎn)距、線距均會對邊界識別結(jié)果產(chǎn)生影響。較大的窗口尺寸會包含較多的觀測信息，但不一定會提高分辨率，甚至?xí)p少邊界信息。對于本文算法，若觀測數(shù)據(jù)的點(diǎn)距、線距過大，將導(dǎo)致MWJC及IBTP處理結(jié)果無法檢測到異常體邊界。因此，建議在實(shí)際工作中應(yīng)用多方向窗口時(shí)，應(yīng)保持窗口內(nèi)盡量少的測點(diǎn)數(shù)，例如3×3，即選取待計(jì)算測點(diǎn)附近的窗，這樣既減少了單個(gè)窗口內(nèi)的計(jì)算量，也避免了大窗口帶來的方向性影響。另外，觀測數(shù)據(jù)的點(diǎn)距和線距一般是根據(jù)具體的勘探任務(wù)與工作標(biāo)準(zhǔn)制定的，無法通過數(shù)據(jù)處理來增加點(diǎn)距、線距。當(dāng)數(shù)據(jù)量較大時(shí)(106量級以上)，如考慮數(shù)據(jù)抽稀，建議數(shù)據(jù)量保持量級在104以上，避免重新網(wǎng)格化導(dǎo)致的點(diǎn)距、線距的增加。

圖12 窗口大小為11×11時(shí)模型二的Ck分布

圖13 窗口大小為11×11時(shí)模型二的Cjoint分布(a)及經(jīng)IBTP處理的結(jié)果(b)

圖14 測點(diǎn)數(shù)為41×41時(shí)模型二的Ck分布

圖15 測點(diǎn)數(shù)為41×41時(shí)模型二的Cjoint分布圖

3 實(shí)測數(shù)據(jù)試驗(yàn)

為檢驗(yàn)MWJC法在實(shí)際應(yīng)用中的效果，將其應(yīng)用于加拿大圣喬治灣地區(qū)的實(shí)測航空重力梯度數(shù)據(jù)(圖16)，該地區(qū)位于紐芬蘭西南海濱的圣勞倫斯海灣。

飛機(jī)的平均飛行高度為129m，經(jīng)網(wǎng)格化后的點(diǎn)距、線距均為100m。根據(jù)3.3節(jié)的分析，MWJC計(jì)算時(shí)窗口設(shè)置為3×3，尺寸為200m×200m，即選取窗口的邊界與待計(jì)算測點(diǎn)的距離不超過200m。利用高斯濾波進(jìn)行去噪。圖17為DTHD分布，可見對地下構(gòu)造邊界識別效果較差，只能初步看出三個(gè)極大值區(qū)域，不能確定異常體或者構(gòu)造的具體范圍。通過計(jì)算相關(guān)系數(shù)C(圖18)，可識別出四個(gè)南西—北東向條帶狀負(fù)值區(qū)；聯(lián)合9窗口的C值計(jì)算Cjoint(圖19a)，發(fā)現(xiàn)圖18中的部分極小值消失，結(jié)合模型試驗(yàn)結(jié)果，認(rèn)為被保存下來的極小值分布區(qū)域即是被突出的構(gòu)造邊界的拐角位置或邊界彎曲部分。根據(jù)統(tǒng)計(jì)結(jié)果(圖19b)設(shè)閾值α=0，經(jīng)IBTP處理后的結(jié)果(圖20)中保留了圖18中極小值的條帶狀分布特征，且進(jìn)一步突出了構(gòu)造邊界的拐角區(qū)域。與文獻(xiàn)[11-12]中其他方法的處理結(jié)果對比，可見識別的邊界分布特征基本一致，進(jìn)一步說明MWJC法具有一定的實(shí)用性，適用于實(shí)測數(shù)據(jù)的邊界識別，且效果良好。

圖16 實(shí)測重力梯度數(shù)據(jù)

圖17 實(shí)測數(shù)據(jù)的THDx(a)和THDy(b)分布

圖18 實(shí)測數(shù)據(jù)的Ck分布

圖19 實(shí)測數(shù)據(jù)的Cjoint分布(a)及Cjoint子區(qū)間數(shù)量統(tǒng)計(jì)柱狀圖(b)

圖20 圖19a經(jīng)IBTP處理的結(jié)果

4 結(jié)論

本文通過x、y方向總水平導(dǎo)數(shù)的相關(guān)系數(shù)實(shí)現(xiàn)了基于多方向窗口聯(lián)合相關(guān)系數(shù)的地質(zhì)體邊界識別。多方向窗口技術(shù)與改進(jìn)二值化與閾值處理方法相結(jié)合，使地下結(jié)構(gòu)的邊界成像更加清晰，有助于識別邊界中的拐角、增強(qiáng)邊界顯示效果。理論模型試驗(yàn)證明：過大的窗口尺寸、觀測點(diǎn)距和線距會降低識別效果，甚至使算法失效，合理選擇相關(guān)參數(shù)有助于更加準(zhǔn)確地識別地質(zhì)體邊界；本文方法能夠較準(zhǔn)確地同時(shí)識別不同深度異常體的邊界。將本文方法應(yīng)用于圣喬治灣地區(qū)實(shí)測重力梯度數(shù)據(jù)，對該地區(qū)的地下構(gòu)造邊界得到較為準(zhǔn)確的識別結(jié)果。

感謝加拿大自然資源部(Natural Resources Canada)為本項(xiàng)研究提供了實(shí)測重力梯度數(shù)據(jù)，允許本文使用此數(shù)據(jù)進(jìn)行研究并發(fā)表研究成果。