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

張 超，王慶賓，黃佳喜，馮進(jìn)凱

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

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

張 超，王慶賓，黃佳喜，馮進(jìn)凱

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

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

重力張量;邊界識(shí)別;水平梯度;解析信號(hào);Sigmoid變換

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

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

1 重力異常及梯度張量

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

(1)

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

(2)

(3)

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

2 邊界增強(qiáng)與識(shí)別算法

2.1 邊界識(shí)別算法

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

對(duì)于各張量組合形式的表達(dá)式見表1和表2。

表1 各張量組合形式的表達(dá)式

表2 常見張量高階組合形式的表達(dá)式

2.2 重力張量數(shù)據(jù)對(duì)比度增強(qiáng)

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

Sigmoid變換方程式原型為

(4)

Sigmoid曲線呈現(xiàn)為S狀，通過將其作用于異常數(shù)據(jù)，可以實(shí)現(xiàn)高異常值網(wǎng)格數(shù)據(jù)的拉伸，同時(shí)壓縮低灰度級(jí)像素[8]。這在一定意義上符合異常體邊界初步增強(qiáng)后，進(jìn)一步提高邊界對(duì)比度的要求，從而凸顯地質(zhì)體的邊界，提高邊界增強(qiáng)后圖像的視覺效果。通過引入合適的參數(shù)，可以在一定程度上改變響應(yīng)曲線的形狀，從而使變換算法更加實(shí)用。引入?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ù)，β則用于控制相位。在實(shí)際應(yīng)用時(shí)，可以根據(jù)當(dāng)前處理的異常數(shù)據(jù)的不同特性選擇適當(dāng)?shù)摩粒轮?，從而達(dá)到有針對(duì)性的改善異常數(shù)據(jù)分布質(zhì)量，實(shí)現(xiàn)對(duì)比度拉升。參數(shù)中的α值決定變換函數(shù)的陡峭程度，而β值則確定變換函數(shù)在坐標(biāo)水平軸上的平移分量。分析可得：當(dāng)β取值很小時(shí)，可以實(shí)現(xiàn)壓縮低異常值部分，拉伸高異常值部分。當(dāng)β取值較大時(shí)，則全圖的所有網(wǎng)格點(diǎn)的值均會(huì)得以拉升，其中高異常值的像素的拉升度要強(qiáng)于低灰階的拉升度。由該函數(shù)性質(zhì)特點(diǎn)可知，當(dāng)異常數(shù)據(jù)中的所有網(wǎng)格點(diǎn)的值均較低時(shí)，采用這種方式得到的增強(qiáng)效果比較明顯。

3 模型試算

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

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

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

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

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

3.2 改進(jìn)二階傾斜角ITA2和改進(jìn)三階傾斜角ITA3對(duì)比

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

(6)

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

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

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

3.3 Sigmoid變換分析

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

圖3 Sigmoid變換結(jié)果

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

4 結(jié) 論

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

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

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

[3] 馬國(guó)慶，李麗麗，杜曉娟．磁張量數(shù)據(jù)的邊界識(shí)別和解釋方法[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方法在位場(chǎng)數(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] 高云軍，劉召芹，萬(wàn)文輝，等.DEM網(wǎng)格間距對(duì)重力近區(qū)地形改正精度影響驗(yàn)證[J].測(cè)繪工程，2016，25(8)：14-18.

[9] 賈魯，王世忠，王慶賓，等.基于Android平臺(tái)的兩種重力垂直梯度解算方法研究[J].測(cè)繪工程，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

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

10.19349/j.cnki.issn1006-7949.2017.04.004

2016-06-09

973國(guó)家高科技計(jì)劃資助項(xiàng)目(6132220202)

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

P223

A

1006-7949(2017)04-0016-06