重力梯度數(shù)據(jù)的聯(lián)合相關系數(shù)邊界識別方法

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

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

0 引言

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

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

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

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

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

1.1 x和y方向總水平導數(shù)的相關系數(shù)

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

(1)

(2)

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

(3)

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

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

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

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

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

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

CF×CG×CH

(4)

Cjoint可以增強THDx與THDy的相關性特征，進一步提高邊界識別效果。

1.3 改進的二值化與閾值處理方法

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

圖1 多方向窗口示意圖

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

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

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

2 理論模型試驗

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

2.1 模型一

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

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

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

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

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

圖4 模型一的Ck分布

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

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

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

2.2 模型二

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

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

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

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

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

圖9 模型二的Ck分布

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

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

2.3 影響識別結果的因素分析

為了進一步檢驗方法的實用性，以模型二為例，討論算法中影響邊界識別結果的因素。

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

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

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

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

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

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

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

3 實測數(shù)據(jù)試驗

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

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

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

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

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

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

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

4 結論

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

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