地下目標體重力張量數據的邊界識別與增強算法的研究

張 超，王慶賓，黃佳喜，馮進凱

(信息工程大學 地理空間信息學院，河南 鄭州 450052)

地下目標體重力張量數據的邊界識別與增強算法的研究

張 超，王慶賓，黃佳喜，馮進凱

(信息工程大學 地理空間信息學院，河南 鄭州 450052)

邊界增強與識別在重力數據處理中占據重要地位，與傳統重力異常數據相比,重力張量及其高階分量對于直接反映異常體的邊界具有更高的精度。當異常數據中的所有網格點的值均較低時,通過Sigmoid變換,可以實現高異常值網格數據的拉升,同時壓縮低灰度級像素,從而凸顯地質體邊界,提高邊界增強后圖像的識別效果。文中利用張量及其分量構建常用的邊界識別算法,通過組合體模型進行多種邊界識別算法的試算，以比較分析各自的效果,并對結果進行Sigmoid變換。結果表明:對于張量高階分量組合形式,水平梯度模、解析信號能基本反映淺異常體的邊界,gzz水平梯度模能較好反映淺異常體邊界,但三者均不能識別深異常體邊界;Tilt梯度、Theta和ITA3效果不理想;ITA2能在有效均衡不同強度異常信號的同時,清晰地識別不同深度異常體的邊界;采用Sigmoid變換,可以明顯提高邊界識別的識別效果。

重力張量;邊界識別;水平梯度;解析信號;Sigmoid變換

重力勘探方法屬于應用地球物理地質勘探與礦藏尋找技術的一個重要分支。由于礦產資源或地下地質構造的密度以及磁性一般同附近的巖體有著顯著的不同，因此，重力勘探在尋找礦產資源與地質構造研究中具有重要的應用價值。重力勘探具有快速、經濟、范圍大等優點，是進行構造劃分、異常圈定和發掘礦產資源等快速而有效的手段[1]。地質體邊界是地質體的斷裂帶、相異地質體之間的分界線。在這些位置處，由于重力能在一定程度上表現出數值的不連續性，呈現出重力異常的較大變化。利用這一特點，可以通過重力異常的較大變化率尋找對應的地質體邊界線，從而為地下地質結構的分析提供支持。

邊界增強與識別在重力數據處理中占據重要地位。原始實測重力異常數據對地下異常體信息反映并不突出，直接根據此數據并不能準確劃分出地質體的邊界。因此需要對原始異常數據進行相應濾波處理來增強對地質體邊界的識別。重力異常水平梯度的極大值和垂直梯度的零值基本反映地質體邊界，邊界提取算法也是基于此提出的[2]。常見的位場數據邊界提取算法主要包括Pedersen等提出的兩個張量不變量組合、水平梯度模法(The total horizontal derivatives,THDR)[3]、解析信號法(The analytic signal amplitude,ASA)[4]、Tilt梯度法[5]、Tilt梯度的水平導數法[6]、Theta圖法[7]、gzz的水平梯度模法等等，然后利用圖中的極值點、零點或其他特征點來識別邊界。本文利用張量及其高階分量構建了常用的邊界識別算法，通過組合體模型進行多種邊界識別算法的試算，比較分析各自的效果，并對結果進行Sigmoid變換。

1 重力異常及梯度張量

建立一個局部坐標系，以地面某一點O為坐標原點，坐標軸ζ軸沿重力方向垂直向下，ξ，η軸在水準面上，以下所有的計算都在局部坐標系下進行。設異常體剩余密度為σ，則引力位

(1)

重力異常就是引力位沿ζ軸方向的導數，即

(2)

(3)

由于在地球外部，位場滿足Laplace方程Uxx+Uyy+Uzz=0 ，以及重力場的無旋性，即Uxy=Uyx，Uxz=Uzx，Uyz=Uzy，因此重力梯度張量中只有5個獨立分量。

2 邊界增強與識別算法

2.1 邊界識別算法

常用的異常數據處理濾波器主要都是基于水平梯度和垂直導數的各種組合來構建，以增強異常地質體的邊界特征，從而實現邊界識別。與傳統重力異常數據相比，重力張量數據對于直接反映異常體的邊界具有更高的精度。通過梯度張量及其高階分量的不同組合，把多個分量的信息集中增強，突出邊界信息，從而輔助數據處理和解釋。

對于各張量組合形式的表達式見表1和表2。

表1 各張量組合形式的表達式

表2 常見張量高階組合形式的表達式

2.2 重力張量數據對比度增強

采用對比度增強方法對重力張量網格數據進行對比度拉伸，增強數據的對比度，擴展其中的強異常與弱異常區域的對比效果，可以增強邊界識別和檢測的效果。這里主要運用Sigmoid變換對網格值進行處理。

Sigmoid變換方程式原型為

(4)

Sigmoid曲線呈現為S狀，通過將其作用于異常數據，可以實現高異常值網格數據的拉伸，同時壓縮低灰度級像素[8]。這在一定意義上符合異常體邊界初步增強后，進一步提高邊界對比度的要求，從而凸顯地質體的邊界，提高邊界增強后圖像的視覺效果。通過引入合適的參數，可以在一定程度上改變響應曲線的形狀，從而使變換算法更加實用。引入參數后的模型為

(5)

式中：fmax為異常數據最大值；I(i,j) 為原位于網格(i,j)的異常值；f(i,j)為網格(i,j)處理結果異常值；α，β為調整變換模型的參數。其中α為用于控制形狀的參數，β則用于控制相位。在實際應用時，可以根據當前處理的異常數據的不同特性選擇適當的α，β值，從而達到有針對性的改善異常數據分布質量，實現對比度拉升。參數中的α值決定變換函數的陡峭程度，而β值則確定變換函數在坐標水平軸上的平移分量。分析可得：當β取值很小時，可以實現壓縮低異常值部分，拉伸高異常值部分。當β取值較大時，則全圖的所有網格點的值均會得以拉升，其中高異常值的像素的拉升度要強于低灰階的拉升度。由該函數性質特點可知，當異常數據中的所有網格點的值均較低時，采用這種方式得到的增強效果比較明顯。

3 模型試算

為充分驗證表2邊界識別算法對地下目標體的識別效果，這里通過建立規則的地下空洞模型進行計算分析。本文所選地面觀測數據區域為401×401點距(點距單位為m)的平面網格面積，并通過規則模型正演計算得到重力異常數據。本文布設一個由球體和長方體構成的組合體模型：其中球體半徑為30 m，中心點坐標為(-110,-110),上頂面埋深為30 m,與圍巖密度差為2.67 g/cm2；長方體的長寬高為120 m×30 m×30 m，中心點坐標為(50,50)，上頂面埋深為15 m，與圍巖密度差為2.67 g/cm2。

3.1 常見張量高階組合形式的邊界識別效果分析

分別利用水平梯度模法、解析信號法、Tilt梯度法、Theta圖法、gzz的水平梯度模法、改進二階傾斜角等方法ITA2對重力異常進行處理，結果見圖1。

圖1 模型重力張量組合形式邊界識別效果

從以上試算結果中可以看出：水平梯度模和解析信號能基本反映出淺異常體平面位置和邊界，但基本不能反映深異常體的邊界；雖然Tilt梯度和Theta圖能有效均衡反映不同深度的異常體信號，但Tilt梯度不能反映異常體的邊界，而Theta圖反映的邊界模糊且發散，并且兩者均出現邊界效應；gzz的水平梯度模能夠較好地反映出淺異常體的邊界，但不能反映深異常體的邊界；改進二階傾斜角ITA2能有效均衡反映不同深度的異常體信號，識別出的邊界略大于真實邊界，對異常體邊界反映效果最好。

3.2 改進二階傾斜角ITA2和改進三階傾斜角ITA3對比

由3.1中試算結果可知，ITA2對不同深度的異常體邊界識別效果最好。理論上，重力張量高階導數對異常體的邊界識別具有更高的分辨率，能更加突出異常體邊界。這里構建ITA3，并比較ITA2和ITA3對不同深度的異常體模型的邊界識別效果。

(6)

試驗模型與3.1中模型一樣，結果如圖2所示。

圖2 ITA2和ITA3比較結果

從圖2可以看出：與ITA2相比，ITA3所識別的淺異常體的邊界形狀清晰準確，識別的深異常體邊界清晰，但卻明顯小于真實異常體邊界，而且產生較大的噪聲和邊界效應。總的來說，識別效果較差，明顯不如ITA2識別效果，可能是由于高階導數的計算增大了噪聲的干擾。

3.3 Sigmoid變換分析

分別對前文中水平梯度模、解析信號、Tilt梯度、Theta圖法、gzz的水平梯度模、改進二階傾斜角ITA2和ITA3進行Sigmoid變換，結果如圖3所示。

圖3 Sigmoid變換結果

從圖3對比可以看出：水平梯度模、解析信號、Tilt梯度、Theta圖法、gzz的水平梯度模、改進二階傾斜角ITA2和ITA3等異常數據的對比度得到明顯的提高，凸顯異常體的邊界，提高邊界識別圖像的識別效果；但是從(d)和(f)可以看出，Sigmoid變換在增強有用異常信號的同時也使得數據中的噪聲被增大。

4 結 論

本文通過對組合體模型試算分析，全面比較多種張量高階組合形式的邊界識別方法的應用效果，并對各種邊界識別算法得到的張量異常數據進行Sigmoid變換以增強數據對比度。分析結果表明：水平梯度模和解析信號能基本反映出淺異常體平面位置和邊界，但不能反映出深異常體的邊界；Tilt梯度和Theta圖能有效均衡反映不同深度的異常體信號，但Tilt梯度不能反映異常體的邊界，而Theta圖反映的邊界模糊且發散，并且兩者均出現邊界效應；gzz的水平梯度模能夠較好地反映出淺異常體的邊界，但不能反映深異常體的邊界；ITA3所識別的淺異常體的邊界形狀清晰準確，但由于高階導數的計算增大噪聲的干擾，從而產生較大的噪聲和邊界效應；改進二階傾斜角ITA2能有效均衡反映不同深度的異常體信號，識別出的邊界略大于真實邊界，對異常體邊界反映效果最好；當異常數據中的所有網格點的值均較低時，采用Sigmoid變換可以凸顯異常體的邊界，提高邊界圖像的識別效果。ITA2和Sigmoid變換能夠在異常體的邊界識別工作中發揮更大的作用，從而促進重力反演解釋工作的順利進行。

[1] 馬國慶．位場(重&磁)及其梯度異常自動解釋方法研究[D]．長春：吉林大學，2013．

[2] 吳亮，劉長弘，王慶賓，等．重力異常對地下異常體邊界的識別[J]．測繪科學，2015，40(12)：257-262．

[3] 馬國慶，李麗麗，杜曉娟．磁張量數據的邊界識別和解釋方法[J]．石油地球物理勘探，2012，47(5)：815-821．

[4] NABIGHIAN M N.The analytic signal of two dimensional magnetic bodies with polygonal cross section：its properties and use for automated anomaly interpretation [J]．Geophysics，2012，37(3)：507．

[5] 王明，郭志宏，駱遙，等．Tilt-Euler方法在位場數據處理及解釋中的應用[J]．物探與化探，2012(2)：126-132．

[6] VERDUZCO B，FAIRHEAD J D，GREEN C M，et al．The meter reader-new insights into magnetic derivatives for structural mapping [J]，Leading Edge，2012 (2)：116．

[7] WIJNS C，PEREZ C，KOWALCZYK P．Theta map：edge detection in magnetic data[J]．Geophysics，2005，70(4)：39-43．

[8] 高云軍，劉召芹，萬文輝，等.DEM網格間距對重力近區地形改正精度影響驗證[J].測繪工程，2016，25(8)：14-18.

[9] 賈魯，王世忠，王慶賓，等.基于Android平臺的兩種重力垂直梯度解算方法研究[J].測繪工程，2016，25(7)：57-61.

[責任編輯：張德福]

Research on underground anomaly boundary recognition and enhancement algorithms based on gravitational tensor

ZHANG Chao, WANG Qingbin, HUANG Jiaxi, FENG Jinkai

(School of Surveying and Mapping, Information Engineering University, Zhengzhou 450052,China)

Boundary enhancement and recognition occupies the important position in magnetic data processing. Compared with the traditional gravity anomaly data, gravitational tensor data has a higher accuracy in directly reflecting the boundary of the anomaly. High abnormal values can be pulled up with low gray-scale pixels compressed at the same time by Sigmoid transformation, when all the values of the grid points are lower, thus highlighting geological boundary and enhancing image recognition effect. This paper calculates many kinds of gravity anomaly processing algorithms based on tensor components, including Sigmoid transformation, and then analyzes their effects through combination model. The result shows that: the horizontal gradient，analytic signal and the horizontal gradient of,gzzcan reflect the boundary of the shallow anomaly with the horizontal gradient of,gzza better effect, but none of them can recognize the boundary of deep anomaly. The result of Tilt gradient and Theta diagram method is not ideal, so it is with the ITA3. By contrast, the ITA2can effectively balance the different intensity anomaly and at the same time clearly identify the different depth of anomaly boundary. Sigmoid transformation is proposed, which can distinctly improve the recognition effect of boundary identification.

gravitational tensor;boundary recognition;horizontal gradient;analytical signal;sigmoid transformation

引用著錄：張 超，王慶賓，黃佳喜，等.地下目標體重力張量數據的邊界識別與增強算法的研究[J].測繪工程，2017，26(4)：16-21.

10.19349/j.cnki.issn1006-7949.2017.04.004

2016-06-09

973國家高科技計劃資助項目(6132220202)

張 超(1990-)，男，碩士研究生.

P223

A

1006-7949(2017)04-0016-06