張恒磊，耿美霞，胡祥云*

（1.中國地質(zhì)大學(xué)（武漢）地質(zhì)探測與評估教育部重點實驗室，湖北武漢 430074； 2.中國地質(zhì)大學(xué)（武漢）地球物理與空間信息學(xué)院，湖北武漢 430074； 3.哈利法科學(xué)技術(shù)大學(xué)，阿聯(lián)酋阿布扎比 127788）

0 引言

隨著研究不斷深入，根據(jù)重磁異常計算地下物性（密度、磁性）分布的反演方法在地質(zhì)勘查［1-5］與深部構(gòu)造［6-12］等領(lǐng)域廣泛應(yīng)用。前人圍繞反演算法開展了深入研究，取得了一系列研究成果［13-15］，不斷豐富、完善重磁反演理論體系。由于三維反演的對象是面積性重磁異常，數(shù)據(jù)量遠遠大于剖面反演，巨大的核函數(shù)存儲需求是制約重磁三維反演實踐應(yīng)用的難題之一。比如：對于數(shù)據(jù)量為100×100 的觀測數(shù)據(jù)，若模型剖分網(wǎng)格數(shù)為100×100×50，其雙精度核函數(shù)的內(nèi)存需求約為37 GB［16］；對于一個數(shù)據(jù)量為200×200的觀測數(shù)據(jù)，若模型剖分網(wǎng)格數(shù)為200×200×100，其雙精度核函數(shù)的內(nèi)存需求高達1192 GB，遠遠超出目前普通計算機的硬件配置。針對重磁異常三維反演，國內(nèi)外學(xué)者開展了廣泛的研究。姚長利等［16］提出獨特的核函數(shù)等效壓縮存儲技術(shù)，降低存儲需求，提高反演效率，為三維重磁異常遺傳反演算法和隨機子域算法提供了有效技術(shù)途徑；陳召曦等［17］提出基于GPU 的任意三維復(fù)雜形體重磁異?？焖僬莘椒?，極大地縮短了計算時間，為三維重磁數(shù)據(jù)快速反演奠定了基礎(chǔ)；李澤林等［18］提出了數(shù)據(jù)空間磁異常模量三維反演方法，降低剩磁影響；戴世坤等［19］提出一種空間—波數(shù)混合域數(shù)值模擬方法，明顯提高了三維重力異常正演效率；袁洋等［20］基于Toeplitz 結(jié)構(gòu)矩陣提出三維模型磁場快速正演方法。學(xué)者們主要從核矩陣小波壓縮［21］以及核矩陣的結(jié)構(gòu)性特點［22］等角度出發(fā)，研究快速三維反演方法。小波壓縮反演算法在實踐中應(yīng)用廣泛，但鮮有關(guān)于壓縮率如何影響反演精度的討論。大規(guī)模重磁數(shù)據(jù)的三維反演對計算機硬件要求較高，因而在實踐中大多盡量采用較少的三維剖分網(wǎng)格數(shù)量，通常為3000～50000［9，23-37］。

本文針對海量重磁異常數(shù)據(jù)壓縮反演開展研究，提出了基于曲波壓縮的重磁異常三維反演方法，分析了不同壓縮率時小波壓縮與曲波壓縮對重磁異常三維反演的影響。結(jié)合理論分析與較大規(guī)模數(shù)據(jù)（模型剖分網(wǎng)格數(shù)為50 萬）的反演測試，指出在相似反演效果條件下，基于曲波壓縮的重磁異常三維反演可以實現(xiàn)更小的壓縮比，有助于提高大規(guī)模重磁異常的三維反演效率。

1 方法原理

1.1 重磁異常三維反演

重磁異常三維物性反演通常是將地下空間劃分為若干小長方體，每個長方體的物性參數(shù)值不同，由觀測數(shù)據(jù)反演每個小長方體的剩余密度或者磁化強度（或磁化率）。其正演問題可以表示成如下形式

式中：N表示觀測點總數(shù)；M表示剖分長方體個數(shù)；m表示剩余密度或磁化強度；d表示觀測異常；Ai，j是核矩陣A的元素，表示第j個長方體對第i個觀測點處產(chǎn)生的異常響應(yīng)。式（1）可表示成矩陣形式

反演問題即是尋找一個合適的模型，使該模型的正演異常與實際觀測異常達到一定的吻合度，同時還需考慮反演模型符合地質(zhì)構(gòu)造規(guī)律。一般目標(biāo)函數(shù)可表示為

式中：φd表示數(shù)據(jù)擬合差；φm表示模型粗糙度；α為正則化因子，作用是平衡數(shù)據(jù)與模型的權(quán)重。Li等［38-39］提出采用有限差分計算模型粗糙度

式中：Wx、Wy、Wz分別表示x、y、z方向的有限差分算子；Ws表示單位矩陣；ax、ay、az、as是權(quán)系數(shù)?？梢?，過小的正則化因子有助于擬合觀測數(shù)據(jù)，但模型粗糙度大，導(dǎo)致反演結(jié)果偏離真實值；而較大的正則化因子使得反演結(jié)果過于平緩而偏離觀測數(shù)據(jù)［21，40］。三個方向的權(quán)系數(shù)ax、ay、az相對于as的大小直接影響反演結(jié)果沿該方向的平滑性［38-40］。為了保證反演效果的可對比性，本文統(tǒng)一取正則化因子α=0.1，加權(quán)因子ax=ay=as=1，az=0.1。

根據(jù)上述公式直接反演得到的密度或磁化強度趨于地表而不是按照異常體的真實深度分布，即重磁反演的“趨膚效應(yīng)”，這是由于構(gòu)造模型的核函數(shù)是線性的，重磁異常幅值隨著距離增大而快速衰減，核矩陣元素的值隨深度增加而急劇減小。因此，可以通過引入深度加權(quán)函數(shù)來克服核函數(shù)隨深度增大的衰減特征，以降低“趨膚效應(yīng)”。Li等［38-39］在重磁異常反演中引入深度加權(quán)函數(shù)

式中：z為塊體單元中心點埋深；z0表示深度加權(quán)控制參數(shù)，其值取決于塊體單元的尺寸以及觀測面的高度；β為深度加權(quán)因子，磁異常反演取β=3，重力異常反演取β=2。

此外，考慮重磁反演的體積效應(yīng)，本文采用統(tǒng)一的約束條件：將每次反演迭代中小于模型增量最大幅值20%的部分歸零，也即忽略模型增量中的小幅值參數(shù)，使得反演結(jié)果趨于模型真實位置。

1.2 小波壓縮核函數(shù)

在三維反演中，核矩陣A對內(nèi)存的需求隨著模型數(shù)量的增大而迅速增大，這也是制約三維反演應(yīng)用的重要因素之一。比如對一個32×32 的網(wǎng)格數(shù)據(jù)體，以等網(wǎng)格距方式將模型剖分為20層，則雙精度矩陣A對內(nèi)存的需求是0.2 GB；對于一個100×100 的網(wǎng)格數(shù)據(jù)體，以等網(wǎng)格距方式將模型剖分為100層，則雙精度矩陣A對內(nèi)存的需求是74.5 GB，考慮解反演方程（式（2））涉及AT以及ATA等矩陣的計算與存儲，其運算量遠超普通計算機的性能。為了減小矩陣A的存儲需求，Li等［21］提出了小波壓縮

式中X表示小波變換算子，其轉(zhuǎn)置矩陣XT即是小波逆變換算子。根據(jù)Li等［21］的研究成果，本文采用db4 小波進行分析和對比。

本文定義壓縮后核矩陣A中非零元素個數(shù)與原核矩陣非零元素總數(shù)之比為壓縮比K?；谛〔ǖ南∈璞磉_性能，Li等［21］指出可以將核矩陣A壓縮至原有的10%左右，即K=0.1。

1.3 曲波變換壓縮核函數(shù)

前人研究表明，因小波僅能表示水平、垂直和對角三個方向的特征，它適于處理點畸奇信號，而難以處理線畸奇信號。曲波變換是一種多尺度變換，是小波變換的擴展，是具有方向性的小波，能夠?qū)?shù)據(jù)更好地進行稀疏表達［41］，因此有助于保持信號中的細節(jié)特征，在地球物理信號去噪、大地電磁二維反演中的邊緣刻畫等領(lǐng)域得到廣泛應(yīng)用［42-52］。

曲波變換是將信號f(x)與曲波函數(shù)φ(x)在二維空間R2內(nèi)進行卷積得到曲波變換系數(shù)［41-42］

因此，與小波變換的系數(shù)分布特征相似，曲波變換系數(shù)中較大的部分表示有效信號的能量。本文將曲波變換引入重磁異常三維反演，用以壓縮核矩陣規(guī)模，即根據(jù)給定的壓縮比K，將小于K的系數(shù)置零，達到矩陣壓縮的目的。

針對小波壓縮和曲波壓縮在核矩陣壓縮精度、正演精度及反演效果三個方面的差別，設(shè)計如下試驗進行對比分析。將模型空間剖分為14637 個水平分布的棱柱體，其中心深度均為200 m，棱柱體大小為20 m×20 m×20 m，磁化強度為1 A/m，地磁場傾角、偏角分別是-1.7°、-2.0°。圖1 是包含119 個測點（觀測面z=0）的南北向剖面磁異常響應(yīng)核矩陣A（矩陣規(guī)模119×14637）及其正演磁異常（d=A×m）的剖面顯示。

圖1 核矩陣A（上）及其正演計算的磁異常（下）

圖2表示不同壓縮比（K=0.10，0.05，0.02）時的壓縮效果?？梢钥闯?，曲波壓縮誤差顯著小于小波壓縮誤差。當(dāng)K=0.10時，小波壓縮和曲波壓縮都可以較好地恢復(fù)核矩陣，核矩陣壓縮相對平均誤差分別為1.06%、0.61%（圖2a）；當(dāng)K=0.05時，壓縮結(jié)果出現(xiàn)一些畸變，相對平均誤差分別為4.01%、1.19%（圖2b）；當(dāng)K=0.02 時，兩種方法的壓縮核矩陣均出現(xiàn)較大誤差（圖2c），相對平均誤差分別為12.72%、2.95%，尤其是小波壓縮損失了很多核函數(shù)的變化特征（如圖2c中的局部放大圖）。

圖2 K=0.10（a）、0.05（b）、0.02（c）時核矩陣壓縮效果

圖3 為對不同壓縮比的結(jié)果進行正演得到的磁異常及其誤差（正演值與真實值之差），可見當(dāng)K=0.10、0.05、0.02時，小波壓縮所對應(yīng)的相對誤差分別是1.46%、5.34%、22.19%，曲波壓縮對應(yīng)的相對誤差分別是0.29%、1.30%、4.63%。因此，采用曲波壓縮能夠更好地保持核矩陣壓縮精度及異常的正演精度，在K取值較小時（比如0.05，0.02），其壓縮效果顯著優(yōu)于小波壓縮。

圖3 不同K 時的磁異常正演曲線（a）及誤差（b）

2 理論模型壓縮反演

為驗證本文方法的壓縮反演效果，設(shè)計圖4 所示的磁異常模型：模型網(wǎng)格數(shù)為100×100，網(wǎng)格間距為20 m；地磁傾角為45°，地磁偏角為0°；兩個鐵礦體A、B分布于向斜兩翼，其上頂埋深分別為200、300 m。不考慮剩磁，該模型產(chǎn)生的磁異常最大值約為1100 nT。

圖4 理論磁異常模型示意圖

為研究地質(zhì)體在地下空間的三維分布特征，將地下空間剖分為100×100×50 個三維網(wǎng)格，剖分單元為邊長20 m 的立方體。據(jù)此得到核函數(shù)A和AT的雙精度存儲量約為74 GB，ATA雙精度存儲量也為74 GB。分別采用Li等［21］的小波壓縮算法及本文曲波壓縮算法進行反演，取K=0.04，即存儲雙精度的A和AT僅需要占用內(nèi)存2.96 GB。反演在內(nèi)存為16 GB、主頻為2.4 GHz 的筆記本電腦上運行。30 次迭代反演用時3.3 h，迭代均方根誤差（RMS）見圖5，基于小波壓縮和曲波壓縮結(jié)果反演的磁化強度分布見圖6。

圖5 圖4 模型反演迭代均方根誤差曲線

圖6 圖4 模型分別基于小波壓縮（上）和曲波壓縮（下）

從圖5 可以看出，基于兩種方法壓縮結(jié)果的正演擬合都可以較快速地收斂于觀測數(shù)據(jù)，曲波擬合精度稍優(yōu)于小波壓縮。當(dāng)壓縮比較低（K=0.04）時，由于小波壓縮引起的核函數(shù)能量損失，反演結(jié)果存在顯著的拖尾振蕩現(xiàn)象（圖6上）；相對而言，基于曲波壓縮的反演結(jié)果具有較高的反演精度，沒有出現(xiàn)顯著的畸變（圖6下）。

考慮Li等［21］在小波壓縮反演中采用的壓縮比為0.12，其反演結(jié)果較光滑，因而不存在圖6（上）的鋸齒振蕩現(xiàn)象。為了分析該鋸齒振蕩是否由較低壓縮率引起，對圖4模型采用壓縮比為0.10進行小波壓縮反演（存儲雙精度的A和AT需要內(nèi)存7.4 GB），在戴爾PowerEdge-T640 工作站上開展計算，用時約12 h。反演結(jié)果沿PP′剖面的垂直切片如圖7 所示，可見與圖6b 下較相似，因此可認為圖6b 上的鋸齒現(xiàn)象是由較低的壓縮比導(dǎo)致信息失真引起的。

圖7 圖4 模型基于小波壓縮反演的磁化強度PP′剖面

3 應(yīng)用實例

根據(jù)中國西部尕林格鐵礦勘探區(qū)的前期地質(zhì)、鉆探研究結(jié)果，該區(qū)為富鐵礦集區(qū)，磁異常ΔT分布見圖8a。該區(qū)巖礦石磁性特征相對簡單：淺部150～200 m 巨厚的第四系覆蓋層無磁性，強磁性鐵礦體（磁化率為6.296 SI）主要賦存于無磁性或弱磁性的安山巖、透輝石矽卡巖（磁化率為0.005～0.013 SI）中，呈透鏡狀、斜列平行產(chǎn)出（圖8b）［53］，為磁法勘探提供了良好的前提條件。該區(qū)地磁場傾角為56.5°、磁偏角為0.2°，磁異常南側(cè)為正，北側(cè)伴生負異常，呈北西西走向。前期見礦鉆孔與未見礦鉆孔難以揭示礦體位置。該區(qū)磁性體剩磁特征不明顯［54］，因而本文不考慮剩磁。

圖8 尕林格鐵礦勘探區(qū)磁異常ΔT（a）以及鉆孔剖面（b）

前期通過改進的Euler 反褶積［53］、基于深度約束的化極磁異常邊界識別Tilt 方法［54］等厘定了該區(qū)磁異常分布及礦體平面位置、埋深等信息（圖9）。推測除主異常①外，在測區(qū)東部存在次級異常②。由于該區(qū)異常平緩，Euler 反褶積對兩個地質(zhì)體未能進行明確的區(qū)分，有關(guān)礦體的三維空間分布及礦體磁性特征亟待深入研究。

圖9 尕林格鐵礦勘探區(qū)基于深度約束的化極磁異常邊界識別結(jié)果（彩色底圖為Tilt）及改進的Euler 反褶積結(jié)果（圖中散點）

圖8a 顯示的磁異常網(wǎng)格數(shù)為100×100，網(wǎng)格間距為20 m（東西向坐標(biāo)范圍為-360～1620 m，南北向坐標(biāo)范圍為-360～1620 m）。地下三維空間剖分為100（東西向）×100（南北向）×50（縱向）個網(wǎng)格單元，單元邊長均為20 m。分別采用小波壓縮算法和本文曲波壓縮算法進行反演，取K=0.03（存儲雙精度的A和AT僅需內(nèi)存2.22 GB），在內(nèi)存為16 GB、主頻為2.4 GHz 的筆記本電腦上運行。30 次迭代反演用時約2.5 h，迭代均方誤差見圖10。

圖10 尕林格鐵礦勘探區(qū)反演迭代均方根誤差曲線

為了直觀地反映并比較磁異常的分布特征，截取了 AA′剖面（位置見圖8a），基于小波壓縮和本文曲波壓縮數(shù)據(jù)的反演結(jié)果分別見圖11a、圖11b?？梢钥闯觯瑑煞N方法在300 m 以淺區(qū)域呈現(xiàn)相似的磁化強度分布特征，在300 m 以深的區(qū)域小波壓縮反演結(jié)果顯示出較強的磁化強度異常。該深部異常是客觀存在的，還是較小的壓縮比導(dǎo)致壓縮失真而引起的反演畸變？為了解答這個問題，采用K=0.10重新進行反演（存儲雙精度的A和AT需要占用內(nèi)存7.4 GB）?？紤]矩陣計算量較大，在戴爾PowerEdge-T640 工作站上運行，用時約13 h。反演結(jié)果見圖11c?？梢?，較大壓縮比（K=0.1）下的小波壓縮反演結(jié)果與較小壓縮比（K=0.03）下的曲波壓縮反演結(jié)果相似，沒有發(fā)現(xiàn)深部異常體。因此，圖11a 中的深部磁化強度異常是由于小波壓縮采用了較小的壓縮比（K=0.03）導(dǎo)致信號失真而出現(xiàn)的虛假異常。

圖11 尕林格鐵礦勘探區(qū)三維磁化強度反演結(jié)果沿AA′剖面的垂直切片

圖12 是基于兩種壓縮數(shù)據(jù)的磁化強度反演結(jié)果在260 m 深度（主礦體深度）上的水平切片?？梢钥闯?，基于曲波壓縮數(shù)據(jù)的三維反演即使在較小的壓縮比（K=0.03）下依然可以保持較穩(wěn)定的反演結(jié)果，反演結(jié)果的異常形態(tài)與采用較大壓縮比的小波壓縮數(shù)據(jù)反演結(jié)果相似，可清晰地反映旁側(cè)異常圈閉②。

圖12 尕林格鐵礦勘探區(qū)磁化強度三維反演結(jié)果在260 m 深度上的水平切片

4 結(jié)論

自小波壓縮引入重磁異常三維反演，重磁三維反演快速發(fā)展并廣泛應(yīng)用于地質(zhì)勘查、深部地球物理等領(lǐng)域，不斷有學(xué)者從計算效率和反演效果兩個方面推陳出新，發(fā)展了一系列重磁三維反演方法，其中小波壓縮反演方法被國內(nèi)外學(xué)者廣泛應(yīng)用并商業(yè)化。盡管如此，較大的核矩陣計算難題依然是制約海量重磁異常三維反演應(yīng)用的瓶頸，也是當(dāng)前學(xué)者關(guān)注的熱點問題之一。本文從稀疏壓縮的角度出發(fā)，將曲波變換應(yīng)用于重磁異常三維反演的核矩陣大規(guī)模壓縮，著重討論了傳統(tǒng)小波壓縮反演方法與曲波壓縮反演方法在壓縮率與壓縮精度方面的差別，分析了較小壓縮比時反演失真的原因，指出相對于傳統(tǒng)小波壓縮反演10%的壓縮比，曲波壓縮反演在3%～4%壓縮比時仍能保持較好的反演穩(wěn)定性，即便普通計算機亦能支持大規(guī)模的三維反演，推動重磁數(shù)據(jù)三維反演的實用化。

感謝加州理工學(xué)院應(yīng)用和計算數(shù)學(xué)中心Curve-Lab 實驗室提供曲波變換源代碼！