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

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

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

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

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

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

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

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

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

1 重力異常及梯度張量

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

(1)

重力異常就是引力位沿ζ軸方向的導(dǎo)數(shù)，即

(2)

(3)

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

2 邊界增強與識別算法

2.1 邊界識別算法

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

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

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

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

2.2 重力張量數(shù)據(jù)對比度增強

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

Sigmoid變換方程式原型為

(4)

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

(5)

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

3 模型試算

為充分驗證表2邊界識別算法對地下目標體的識別效果，這里通過建立規(guī)則的地下空洞模型進行計算分析。本文所選地面觀測數(shù)據(jù)區(qū)域為401×401點距(點距單位為m)的平面網(wǎng)格面積，并通過規(guī)則模型正演計算得到重力異常數(shù)據(jù)。本文布設(shè)一個由球體和長方體構(gòu)成的組合體模型：其中球體半徑為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對重力異常進行處理，結(jié)果見圖1。

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

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

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

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

(6)

試驗?zāi)Ｐ团c3.1中模型一樣，結(jié)果如圖2所示。

圖2 ITA2和ITA3比較結(jié)果

從圖2可以看出：與ITA2相比，ITA3所識別的淺異常體的邊界形狀清晰準確，識別的深異常體邊界清晰，但卻明顯小于真實異常體邊界，而且產(chǎn)生較大的噪聲和邊界效應(yīng)?？偟膩碚f，識別效果較差，明顯不如ITA2識別效果，可能是由于高階導(dǎo)數(shù)的計算增大了噪聲的干擾。

3.3 Sigmoid變換分析

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

圖3 Sigmoid變換結(jié)果

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

4 結(jié) 論

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

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

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

[3] 馬國慶，李麗麗，杜曉娟．磁張量數(shù)據(jù)的邊界識別和解釋方法[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方法在位場數(shù)據(jù)處理及解釋中的應(yīng)用[J]．物探與化探，2012(2)：126-132．

[6] VERDUZCO B，F(xiàn)AIRHEAD 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網(wǎng)格間距對重力近區(qū)地形改正精度影響驗證[J].測繪工程，2016，25(8)：14-18.

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

[責(zé)任編輯：張德福]

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

引用著錄：張 超，王慶賓，黃佳喜，等.地下目標體重力張量數(shù)據(jù)的邊界識別與增強算法的研究[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