楊 劍， 王 橋， 郭 鏡， 曾琴琴

(中國地質調查局 成都地質調查中心，成都 610081)

0 引言

地球物理勘探得到的位場數據是地下不同深度層次、規(guī)模、形態(tài)的場源產生異常的疊加，不同深度地質體產生的異常疊加到區(qū)域場之上，使得觀測的數據具有非線性、非平穩(wěn)性的特點[1]。在實際工作過程中，位場數據異常存在著線性特征，這往往能引起地質工作者的興趣，由于這類異常對應著地下斷裂構造、不同巖性地質體的邊界接觸帶或者其他具有一定密度或磁性差異的構造特征。對這些線性特征進行增強、提取并進行半定量地解釋是重磁資料處理的主要內容。然而，位場數據中若混入噪聲干擾，再利用總水平導數、總水平梯度傾斜角、解析信號振幅等方法進行高次求導運算會把噪聲放大，導致提取的線性構造位置發(fā)生偏離甚至出現錯誤[2-3]。這里以小波理論為基礎，提出一種新的位場邊界識別方法—小波模極大值法。從而脫離對位場數據進行導數的運算，避免了求位場數據高階導數帶來的不穩(wěn)定性，其抗噪性使得該方法在重磁異常線性特征提取結果更準確、穩(wěn)定[4]。

1 小波模極大值方法原理

模極大值法，是根據信號和噪聲在多尺度空間上小波變換系數的模極值傳播規(guī)律的不同而發(fā)展起來的一種去噪算法[5]。

設二維平滑函數θ(u,v)滿足：

(1)

由θ(u,v)定義兩個二維小波：

(2)

則f(u,v)在尺度s上的二維小波變換包括兩個分量：

(3)

(4)

在尺度s=2j的離散情況下，則對于任意函數f(u,v)∈L2(R2)在尺度2j下沿想，x、y方向的小波變換定義如下：

(5)

離散小波變換的模定義如下：

(6)

其中:j為分解層數(0≤j≤log2n)；n為信號長度。多尺度小波模極大值法邊界檢測就是通過追蹤小波變換模極大值點來實現的。

為了快速地查詢小波變換后的模極大值點，定義相角為模的梯度方向與水平方向夾角：

(7)

該值對數據中模極大值的快速判定有很大作用。

2 模型仿真試驗

為了從理論上驗證小波模極大值處理方法的有效性，建立了理論模型仿真試驗，模型由7個組合長方體構成(模型參數見表1)，它是由不同形態(tài)、不同埋深、不同剩余密度的長方體組成，正演得到的重力異常見圖1(a)，添加2%高斯白噪聲得到加有噪聲模型重力異常見圖1(b)。

表1 長方體組合模型參數

圖1 組合模型正演重力異常Fig.1 Forward gravity anomaly of model(a)未加噪聲模型重力異常；(b)加有2.5%高斯白噪聲模型重力異常

圖2 TAHG邊界識別結果Fig.2 Boundary recognition results with TAHG method(a)未加噪聲異常；(b)加有噪聲異常

圖3 小波模極大值邊界識別結果Fig.3 Boundary recognition results with wavelet modular maximum(a)未加噪聲異常；(b)加有噪聲異常

表2 兩種不同方法求得的信噪比和均方根誤差

在利用小波模極大值方法處理平面重磁異常時，為了檢驗其抗噪聲干擾能力和線性構造位置的準確度，利用TAHG和小波模極大值方法對未加噪聲和加有噪聲的正演重力異常分別進行處理，檢測結果分別見圖2、圖3，為了從數值上定量的反映小波模極大值去噪效果，比較了TAHG方法和小波模極大值方法的信噪比(SNR)，均方根誤差(RMSE)，信噪比越高，原始信號與估計信號的均方根誤差越小，則消噪信號越接近于原始信號，去噪效果越好。表2列出了用這2種方法進行去噪后，信號的SNR和RMSE的比較，從信噪比和均方根誤差來看，小波模極大值法的去噪性能最好。

圖4 北衙礦區(qū)地質圖Fig.4 Geological map of the Beiya deposit

從圖2、圖3中可以看出，對未添加噪聲的重力異常，TAHG(圖2(a))和小波模極大值(圖3(a))均可以識別出7個長方體模型的異常邊界；在添加噪聲的情況下，TAHG(圖2(b))識別的結果已經面目全非，線性特征不明顯，而小波模極大值方法(圖3(b))仍能夠準確、清晰地獲取組合長方體模型的邊界線性特征，與未添加噪聲計算的結果大同小異，只是在細節(jié)上存在差異。

由上述模型試驗可知，小波模極大值方法是一種有效的位場數據濾波算法，它能夠抵抗噪聲干擾來提高識別線性特征的準確度，進而可以準確地獲取地下斷裂構造或地質體邊界的位置、走向及分布范圍。

3 小波模極大值方法的應用

3.1 礦區(qū)地質、地球物理背景

北衙礦區(qū)地處揚子準地臺西緣的麗江臺緣拗陷褶皺帶，屬云南揚子成礦區(qū)中南段的大理-麗江成礦集中區(qū)(圖4)，是典型的矽卡巖-熱液型礦床，礦區(qū)出露的地層有二疊系上統(tǒng)峨眉山組、三疊系下統(tǒng)青天堡組砂巖、三疊系中統(tǒng)北衙組碳酸鹽巖、第四系。北衙礦區(qū)出露巖漿巖以喜馬拉雅山期形成的淺成侵入富堿斑巖為主，與圍巖北衙組碳酸鹽巖接觸交代，形成鐵金礦體，地層近南北向斷裂組為礦區(qū)內控礦斷層，近東西向斷裂組為破礦斷裂[6-8]，成礦受巖體-地層-構造控制，成“三位一體”。

從礦區(qū)地表地質調查及萬硐山礦段露天開采發(fā)現，北衙礦區(qū)鐵金礦床呈現隱伏特征，是典型的隱伏矽卡巖性礦床，通過有效的地球物理方法尋找礦化異常及巖體低重異常，是該區(qū)重要的找礦手段[9-10]。

通過北衙礦區(qū)巖礦石物性統(tǒng)計得知(表3，圖5)，斑巖的密度小于沉積巖的密度，在研究區(qū)屬低密度體，可產生低重力異常，利用巖體的低密度特性開展重力調查，以圈定巖體在平面的位置。鐵礦石和玄武

表3 巖(礦)石物性標本成果統(tǒng)計表

圖5 北衙礦區(qū)物性柱狀圖Fig.5 Physical histogram of Beiya deposit

巖是高磁體，沉積巖及斑巖磁性較弱，利用高精度磁法可劃分玄武巖地層和沉積巖地層邊界，運用磁異常信息進一步確定礦床的有利富集部位[11-12]。

3.2 應用效果

為了圈定巖體產生的低重異常以及磁鐵礦體產生的高磁異常，在北衙礦區(qū)開展了1:10 000地面高精度重力調查和地面高精度磁測，測區(qū)總面積42 km2。測點重力值計算進行了固體潮改正和零點位移改正，計算方法按規(guī)范執(zhí)行。各單項計算均計算到0.001×10-5m/s2。計算結果經100%復算無誤差后提供使用。野外各種原始資料進行了嚴格驗收和及時整理，各項外部改正及精度滿足后期處理要求。

圖6 滑動平均求取剩余重力異常Fig.6 Residual gravity anomaly with moving average

圖7 小波模極大值求取剩余重力異常Fig.7 Regional gravity anomaly with wavelet modular maximum

在重力資料處理過程中對比了兩種不同的處理方法，利用滑動平均法(圖6)以及小波模極大值法(圖7)求取了剩余重力異常，圖6中萬硐山礦段已經露天開采，低重異常A已經證實是斑巖體引起，根據地球物理從“已知”到“未知”的原則[13-15]，推斷異常B-D，但是后期工程驗證異常B，28ZK01并未打到隱伏巖體，從小波模極大值法求取的剩余重力異常(圖7)上看，28ZK01處低重異常并不是很明顯，此處低重異常B是由于地表第四系較厚，在滑動平均求取剩余重力異常過程中混入了噪聲干擾，在進行高次求導運算中把噪聲放大，形成了低重的 “假異?！?，根據小波模極大值求取的剩余重力異常重新圈定了巖體在地表的位置走向(圖7中粉色虛線)，并切去了兩條剖面進行二維反演推斷了深部巖體走向，后期經過工程驗證與我們推斷的巖體基本一致，其中，B-B'剖面與55號勘探線基本重合，從勘探線剖面圖(圖10)看出，大沙地巖體和紅泥塘巖體與B-B'剖面重力反演推斷的結果基本吻合，證明了我們采用小波模極大值法處理的剩余重力異常是行之有效的。

為了獲取磁鐵礦體產生的局部磁異常以及劃分玄武巖邊界，在磁法資料處理過程中對比了趨勢分析法與小波模極大值方法[16-17]。由于礦區(qū)內人文干擾嚴重，受玄武巖以及地表紅土型第四系剩磁影響，增加了磁鐵礦產生的高磁異常的篩選難度，通過對照地表地質調查結果以及萬硐山露天開采勘探結果，發(fā)現常規(guī)的趨勢分析方法處理的局部異常(圖11)分辨率較低，已知的萬硐山礦段和紅泥塘礦段分布大面積的高磁異常且連為一體，雖然異常較完整且峰值高達3 000 nT，但是根據前期地質路線調查及槽探工程控制情況，紅泥塘礦段并沒有如此大的規(guī)模。另一方面，地層解譯過程中發(fā)現，礦區(qū)東部出露的玄武巖產生的高磁異常邊界在常規(guī)的趨勢分析方法識別的較模糊，界線不明顯。而通過小波模極大值方法處理的局部磁異常(圖12)，在萬硐山和紅泥塘已知礦段呈現更加明顯的“干凈、完整”高磁異常，異常中心明顯，峰值高。利用小波模極大值方法處理的局部磁異常結果圈定的礦化體邊界位置，為工程鉆探提供了合理依據。從后期搜集紅泥塘礦段勘探線55線剖面圖(圖13)看出，55ZK10控制了礦體的中部頂底面，0 m～3.56 m為黃褐色第四系坡積層，主要由粘土，灰?guī)r碎塊等組成；3.56 m～145.48 m為淺灰色鐵化灰?guī)r，細晶結構，塊狀構造，巖石礦物成分由方解石，少量砂屑，少量鐵質和白云石組成；145.48 m～160.79 m為黃褐色構造角礫巖，碎裂結構，角礫狀構造，呈棱角狀；160.79 m～200.54 m為灰白色石英正長斑巖，斑狀結構，塊狀構造，巖石斑晶由石英，長石組成，粒徑為1 mm～5 mm不等，巖石裂隙發(fā)育，可見沿裂隙發(fā)育的褐鐵礦，局部見黃鐵礦的集合體。呈星點狀，脈狀產出，巖芯較破碎，褐鐵礦化較強，于195.96 m處見長0.37 m的灰褐色褐鐵礦，品位較高；200.54 m～237.37 m為灰白色，黃色構造角礫巖，角礫主要為石英正長斑巖和少量灰?guī)r，并有少量褐鐵礦碎塊，角礫粒徑為5 mm～40 mm不等，大小不均勻，與泥質一同產出，呈棱角狀，于217.6 m處見長1.63 m的白云巖，質地較純；237.37 m～330.07 m為灰黑色，灰綠色磁鐵礦化，黃鐵礦化，黃銅礦化矽卡巖，變晶結構，塊狀構造，巖石主要礦物成分為輝石，方解石，石榴子石以及黃銅礦，黃鐵礦和磁鐵礦等金屬礦物，整體蝕變較強，黃鐵礦呈星點狀團包狀產出，黃銅礦成團包狀分布在黃鐵礦邊緣，于328.85 m開始磁鐵礦含量較高；以下至580.86 m為灰色蝕變雜砂巖，粒狀結構，塊狀構造，巖石礦物主要有石英、長石。另外從55ZK16和55ZK14鉆孔數據上看，控制的礦體東西邊界與平面磁異常位置對應性非常的好，證明了磁異常的來源正是深部磁鐵礦體產生；同時，礦區(qū)東部玄武巖產生的高磁異常界線也更加清晰(圖中黑色虛線框為地表地質填圖界線)，這也為地表地質工作提供了有效依據[18]。

圖8 A—A’反演擬合計算剖面圖Fig.8 Density inversion profile map of prospecting A—A’

圖9 B—B’反演擬合計算剖面圖Fig.9 Density inversion profile map of prospecting B—B’

通過以上模型試驗及實際資料處理表明，小波模極大值方法是一種有效的位場數據線性處理方法，在位場分離和邊界識別中的應用是成功的，該方法的優(yōu)越性體現在它能夠抵抗噪聲干擾來提高提取線性特征的準確度，進而可以準確地獲取地下斷裂構造或地質體邊界的位置、走向及分布范圍，提高了

圖10 55勘探線剖面圖Fig.10 Prospecting 55 line profile map

圖11 趨勢分析處理局部磁異常Fig.11 Magnetic anomalies with trend analysis

圖12 小波模極大值處理局部磁異常Fig.12 Magnetic anomalies with wavelet modular maximum

異常篩選的分辨率，增強了異常邊界識別度，在北衙礦區(qū)找礦中起到了非常好的效果。

4 結論

針對位場數據中混入噪聲導致地質體邊界識別不清楚的問題，我們以小波理論為基礎，提出一種新的位場邊界識別方法—小波模極大值法。通過模型試驗及實測數據應用，該方法脫離對位場數據進行導數的運算，避免了求位場數據高階導數帶來的不穩(wěn)定性，抗噪性能強，提高提取線性特征及地質體邊界的準確度，進而可以準確地獲取地下斷裂構造或地質體邊界的位置、走向及分布范圍，是一種具有實際應用效果的位場分離處理方法。

感謝

云南黃金集團有限公司提供原始數據，感謝中國地質大學(武漢)的合作，感謝中國地震局第一監(jiān)測中心張雙喜工程師的協助，感謝審稿老師提出的寶貴意見。

圖13 55線勘探線剖面圖(礦體段)Fig.13 Prospecting 55 line profile map