曲波變換在位場信號提取中的應(yīng)用研究

張楊，王君恒，曹煉鵬，馮裕華，朱江皇，付強(qiáng)

(1.中國地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院，北京 100083； 2.深圳市地質(zhì)環(huán)境監(jiān)測中心，廣東 深圳 518034)

0 引言

重磁勘探技術(shù)是地球物理勘探中重要的方法之一。近年來，隨著深部探測和航空重磁勘探技術(shù)的發(fā)展，對有效信號分離提取的需求越來越大。大區(qū)域產(chǎn)生的背景場中，經(jīng)常包含微小的局部地質(zhì)體產(chǎn)生的異常，因此，如何在密度差異或磁性差異不明顯和存在強(qiáng)干擾的情況下，劃分出不同地質(zhì)體的邊界、消除干擾的背景場、去除噪聲的影響，己成為急需解決的問題[1-4]。

Curvelet變換是Candès和Donoho在1999年最先提出的方法[5]，為了讓Curvelet變換更容易被人們理解和使用，Candès等先后提出第二代Curvelet變換框架體系，提出兩種以第二代Curvelet變換理論為基礎(chǔ)的快速離散實(shí)現(xiàn)方法。新算法理論的提出，使得Curvelet變換的實(shí)現(xiàn)更簡單、快速，減少了初代算法帶來的復(fù)雜冗余[6-9]。

在國內(nèi)外，許多學(xué)者將Curvelet變換運(yùn)用于各種數(shù)據(jù)處理，如圖像數(shù)據(jù)、地震資料數(shù)據(jù)等[10-21]。國內(nèi)學(xué)者將Curvelet變換運(yùn)用于位場數(shù)據(jù)處理方面的研究較少，陳召曦通過建立簡單、復(fù)雜組合模型對Curvelet變換的多尺度、多方向性進(jìn)行了分析[22]。郇恒飛研究了Curvelet變換和小波閾值去噪方法在位場數(shù)據(jù)處理中的應(yīng)用，通過模型試驗(yàn)對兩種去噪方法的效果進(jìn)行對比[23]。高鐵等進(jìn)行了 Curvelet 變換去噪和小波閾值去噪的實(shí)驗(yàn)對比，驗(yàn)證了 Curvelet 變換去噪的有效性[24]。楊斯涵研究Curvelet變換的多尺度分析性質(zhì)及其去噪原理，得出Curvelet多尺度分解后各層系數(shù)重構(gòu)數(shù)據(jù)可以與地下不同深度規(guī)模異常體相對應(yīng)的結(jié)論[25]。

筆者首先介紹了Curvelet變換的基本原理，然后在前人研究的基礎(chǔ)上開發(fā)了相關(guān)Curvelet變換軟件，并且同時利用重力位場理論模型、加噪模型和南嶺東部地區(qū)實(shí)際資料分析了Curvelet變換的分解和去噪能力。綜合前人的研究成果，結(jié)果表明該方法可同時適用于位場數(shù)據(jù)的多尺度分解和去噪處理研究，為實(shí)際重磁數(shù)據(jù)的多尺度分析處理提供參考。

1 第二代Curvelet變換

第二代曲波變換是在小波變換、脊波變換和第一代曲波變換的基礎(chǔ)上構(gòu)建出來的。

二維空間R2中，x為空間的位置變量，ω為頻率域變量，r、θ為對應(yīng)頻率域下的極坐標(biāo)。假設(shè)存在平滑、實(shí)值非負(fù)的“角度窗”V(t)(支撐域?yàn)閠∈[-1,1])和“半徑窗”W(r)(支撐域區(qū)間為r∈(1/2,2))，且滿足容許性條件:

(1)

對于每一個尺度j≥j0可以引入頻率窗口Uj，定義頻率域窗函數(shù)為：

(2)

式中：[j/2]表示j/2整數(shù)部分；頻率域窗函數(shù)Uj為極坐標(biāo)下受“半徑窗”W和“角度窗”V支撐區(qū)間限定的楔形區(qū)域。

如果φj(x)的傅立葉變換滿足等式

(3)

則定義φj(x)為母曲波函數(shù)。所以，在尺度2-j上所有的曲波函數(shù)都可以由φj(x)經(jīng)過旋轉(zhuǎn)與平移得到。

(4)

所以任何一個函數(shù)f∈L2(R2)的曲波變換系數(shù)c(j,l,k)可以通過此函數(shù)f與某一尺度上的曲波函數(shù)φj,l,k的內(nèi)積來表示：

(5)

將地球物理數(shù)據(jù)應(yīng)用于曲波變換時，地球物理數(shù)據(jù)可看作函數(shù)f，如本文使用的重力位場數(shù)據(jù)；將重力位場數(shù)據(jù)函數(shù)f與曲波函數(shù)φj,l,k進(jìn)行內(nèi)積c(j,l,k)=〈f,φj,l,k〉，就可以得到不同尺度下的曲波系數(shù)c(j,l,k)；對得到的曲波系數(shù)c(j,l,k)進(jìn)行處理并反曲波變換重構(gòu)，就可得到處理后的重力位場數(shù)據(jù)函數(shù)f′。

2 Curvelet變換應(yīng)用

2.1 位場數(shù)據(jù)模型分離重構(gòu)試驗(yàn)

二維數(shù)據(jù)函數(shù)f經(jīng)過Curvelet多尺度變換分解后，在不同方向和不同分解尺度上都有對應(yīng)的一組系數(shù)：

(6)

式中：Cj,0(k1,k2)為分解后原結(jié)果的低頻系數(shù)，相當(dāng)于數(shù)據(jù)的粗尺度部分，在該頻帶上無方向性；Cj,l(k1,k2)為尺度j，方向l下的高頻系數(shù)，相當(dāng)于數(shù)據(jù)的精細(xì)尺度部分。按照需要對不同尺度的系數(shù)進(jìn)行處理，以達(dá)到實(shí)現(xiàn)數(shù)據(jù)的融合、去噪、多尺度、多方向等處理的要求，然后對處理后的系數(shù)進(jìn)行選擇性的反Curvelet變換重構(gòu)，就可以得到處理后的結(jié)果。對于512×512的二維數(shù)據(jù)，經(jīng)過Curvelet變換后，可以得到Curvelet系數(shù)如表1所示。

表1 512×512數(shù)據(jù)的Curvelet系數(shù)結(jié)構(gòu)

隨著層數(shù)的增加，反Curvelet變換重構(gòu)的圖像越精細(xì)，反映的越是高頻成分和數(shù)據(jù)的細(xì)節(jié)部分，層數(shù)越低，反映的越是粗糙成分和低頻成分。使用Curvelet變換就是把原始數(shù)據(jù)按照數(shù)據(jù)的大小多尺度分解為不同的層數(shù)，默認(rèn)的層數(shù)大小設(shè)置為 ceil(log2min(M,N)-3)。對于512×512的二維數(shù)據(jù)，Curvelet變換分解為6層系數(shù)。

將Curvelet變換應(yīng)用到重力數(shù)據(jù)，設(shè)置了5個場源體(直立六面體)和觀測面組成的理論模型，觀測面為z=0 km的平面，觀測面z坐標(biāo)方向向下為正。5個直立六面體場源分別編號為A1、A2、A3、B1和B2，其幾何參數(shù)及密度參數(shù)見表2所示。將B1、B2看作深部區(qū)域異常，A1、A2、A3看作淺部局部異常，進(jìn)行信號的分離提取。重力異常見圖1。

表2 直立六面體場源參數(shù)統(tǒng)計

圖1 直立六面體組合重力異常Fig.1 Vertical hexahedron combined gravity anomaly

對圖1a重力數(shù)據(jù)進(jìn)行Curvelet變換，分別提取6層的系數(shù)進(jìn)行反Curvelet變換重構(gòu)， Curvelet變換的多尺度分解性質(zhì)能夠完成位場數(shù)據(jù)的多尺度分解，可以對不同規(guī)模地質(zhì)體所產(chǎn)生的不同尺度的重力異常進(jìn)行分離，來達(dá)到重力位場多尺度分解的目的。規(guī)模大且埋藏深的地質(zhì)體，重力異常通常表現(xiàn)為平緩的區(qū)域異常，經(jīng) Curvelet變換后得到低層系數(shù)；規(guī)模小、埋藏淺的地質(zhì)體重力異常通常表現(xiàn)為變化劇烈的局部異常，經(jīng) Curvelet變換后得到高層系數(shù)(圖2)。

圖2 Curvelet變換分解后的系數(shù)重構(gòu)重力異常Fig.2 Reconstructed graph of gravity anomaly with coefficients decomposed by Curvelet transform

第5層和第6層系數(shù)反映Curvelet變換重構(gòu)圖像數(shù)據(jù)近似為0，不能反映數(shù)據(jù)信息。將第1層圖象定義為分解后得到的深層重力異常，重力平面圖反映了埋深最深、規(guī)模最大的兩個模型B1、B2的異常，消除了A1、A2模型的影響，但包含A3模型的影響，因?yàn)锳3模型產(chǎn)生的重力異常幅值大，其重力異常在深層與淺層都有所體現(xiàn)。將第2層到第4層Curvelet系數(shù)重構(gòu)圖定義為淺層重力異常，包含埋深最淺、規(guī)模最小A1、A2、A3模型的異常。與圖1對比，可以看到Curvelet系數(shù)第1層重構(gòu)圖體現(xiàn)了區(qū)域重力異常，反映深部數(shù)據(jù)，Curvelet系數(shù)第2層到第4層重構(gòu)圖體現(xiàn)了局部重力異常，反映淺部數(shù)據(jù)。

將A3模型的z坐標(biāo)z1=0.5 km、z2=1.0 km改為z1=1.0 km、z2=1.4 km，原始疊加模型重力異常Curvelet變換重構(gòu)圖像見圖3。

圖3 Curvelet變換分解后的重構(gòu)重力異常Fig.3 The reconstructed gravity anomaly map after the decomposition of the Curvelet transform

第5層和第6層系數(shù)重構(gòu)圖像近似為0，不能反映數(shù)據(jù)信息。將第1層圖象定義為Curvelet分解后得到的深層重力異常，重力平面圖反映了埋深最深、規(guī)模最大的兩個模型B1、B2的重力異常，A3模型的重力異常影響被削弱，消除了A2、A1模型的影響，反映的是區(qū)域異常。將第2層系數(shù)重構(gòu)圖定義為中層重力異常，去掉了B1、B2模型體的重力異常，包含A1、A2、A3模型體產(chǎn)生的重力異常。將第3、4層系數(shù)重構(gòu)結(jié)果定義為淺層重力異常，包含埋深最淺、規(guī)模最小A1、A2模型體的重力異常，去掉了A3模型體的影響。從分離結(jié)果來看，Curvelet變換有效地將不同深度層上的重力數(shù)據(jù)進(jìn)行了分離，從而完成了位場數(shù)據(jù)的分解。

圖3與圖2對比，A3模型可以看成中層重力模型，較B1、B2模型埋深淺，較A2、A1模型埋深深，故圖3第1層重構(gòu)圖與圖2第1層重構(gòu)圖相比，A3模型影響減弱；圖3第3層和第4層重構(gòu)圖與圖2第3層和第4層重構(gòu)圖相比，主要反映了A1、A2模型重力異常，A3模型重力異常被消除。所以，通過建立不同埋深的異常體，驗(yàn)證了該方法的分層效果。該方法分層結(jié)果與異常體大小、埋深和密度有關(guān)，分解的系數(shù)層數(shù)越高，反映的越是淺層異常或局部異常；系數(shù)層數(shù)越低，反映的越是深層異?；騾^(qū)域異常。

2.2 閾值去噪模型試驗(yàn)

Curvelet變換去噪的原理是，數(shù)據(jù)經(jīng)過Curvelet變換之后，噪聲通常在變換后會表現(xiàn)為絕對值很小的系數(shù)，所以，我們可以引入閾值分割的理念，對分解后每一層系數(shù)進(jìn)行閾值分割，將數(shù)值較小的系數(shù)設(shè)置為0，對數(shù)值較大的系數(shù)則進(jìn)行保留。實(shí)際數(shù)據(jù)中的有效信息通常在Curvelet變換后表現(xiàn)為數(shù)值較大的系數(shù)。因此，通過這種閾值判斷，可以保留原始數(shù)據(jù)中的重要信息，進(jìn)而消除噪聲的干擾，再將處理后的系數(shù)進(jìn)行反Curvelet變換重構(gòu)，就達(dá)到了去噪的目的。此外，由于曲波變換實(shí)現(xiàn)了不同尺度、層次數(shù)據(jù)的分解，對分解后各層系數(shù)進(jìn)行分析，可以發(fā)現(xiàn)，噪聲通常集中在最高幾層系數(shù)中，有時可以直接將具有嚴(yán)重噪聲的系數(shù)層去掉，便可實(shí)現(xiàn)噪聲濾除。

閾值去噪方法包含硬閾值、軟閾值、自適應(yīng)閾值等。此外，還可以根據(jù)不同情況對Curvelet變換系數(shù)進(jìn)行處理來達(dá)到濾波的目的。

本文采用硬閾值去噪方法來處理位場數(shù)據(jù)資料，硬閾值去噪處理的方法原理是將大于設(shè)定閾值的系數(shù)保留，而小于閾值的系數(shù)設(shè)置為零。

(7)

式中：Cp為原始數(shù)據(jù)的Curvelet變換系數(shù)；Tm是設(shè)

對原始數(shù)據(jù)進(jìn)行去噪可以分為3個步驟:

1) 對原始加噪數(shù)據(jù)進(jìn)行Curvelet變換多尺度分解;

2) 對Curvelet變換分層系數(shù)采用硬閾值法，對每層系數(shù)進(jìn)行處理;

3) 將處理過后的Curvelet變換系數(shù)進(jìn)行反Curvelet變換，就可以得到經(jīng)過閾值處理后的數(shù)據(jù)重構(gòu)圖。

將表2中重力位場數(shù)據(jù)加5%重力幅值的隨機(jī)噪聲進(jìn)行Curvelet變換處理(圖4)，可以看出噪聲集中在4、5層，即高層系數(shù)重構(gòu)數(shù)據(jù)中，第3層中也有少部分噪聲干擾。

圖4 加噪數(shù)據(jù)曲波變換重構(gòu)Fig.4 The reconstructed image with noise-added data using Curvelet transform

去噪是通過對Curvelet系數(shù)進(jìn)行硬閾值處理來實(shí)現(xiàn)的，最精細(xì)的尺度層被設(shè)置為4×sigma_jl閾值，其他層數(shù)選擇3×sigma_jl的閾值，這里 sigma_jl是規(guī)模為j和角度l的系數(shù)的噪聲水平(等于噪聲水平乘以相應(yīng)曲波系數(shù)的L2范數(shù))。

對原始重力異常數(shù)據(jù)加噪聲，然后應(yīng)用Curvelet變換閾值去噪(見圖5～7)，由以上模型的去噪實(shí)驗(yàn)可以看出，Curvelet變換閾值去噪方法效果良好，可以去掉不同幅值噪聲的影響。

圖5 加噪數(shù)據(jù)Curvelet變換閾值去噪(加10%幅值隨機(jī)噪聲,sigma=0.8)Fig.5 Curvelet transform threshold denoising of noise-added data(add 10% amplitude random noise，sigma=0.8)

圖6 加噪數(shù)據(jù)Curvelet變換閾值去噪(加5%幅值隨機(jī)噪聲,sigma=0.45)Fig.6 Curvelet transform threshold denoising of noise-added data(add 5% amplitude random noise，sigma=0.45)

圖7 加噪數(shù)據(jù)Curvelet變換閾值去噪(加1%幅值隨機(jī)噪聲,sigma=0.15)Fig.7 Curvelet transform threshold denoising of noise-added data(add 1% amplitude random noise，sigma=0.15)

2.3 實(shí)際資料處理

南嶺位于華南地區(qū)的中南部，廣泛發(fā)育不同成因、時代的花崗巖，盛產(chǎn)與花崗巖有密切成因關(guān)系的鎢、錫稀有金屬礦產(chǎn)，其成礦條件獨(dú)特，而且成礦潛力巨大，專家學(xué)者對該區(qū)域的研究歷史十分悠久[26-30]。

本次處理的重力位場實(shí)際資料是南嶺東部區(qū)域EGM2008布格重力異常，數(shù)據(jù)范圍：經(jīng)度114°～116°，緯度24°30′～ 26°，網(wǎng)格化成大概1′× 1′的間距(見圖8)。研究區(qū)地形見圖9。

圖8 南嶺東部EGM2008布格重力異常Fig.8 Bouguer gravity anomaly map of EGM2008 in the eastern part of Nanling

圖9 研究區(qū)地形Fig.9 Topographic map of the study area

南嶺東部EGM2008布格重力異常是該區(qū)域地質(zhì)體信息的綜合反映，包括了殼內(nèi)各種偏離正常密度分布的礦體與構(gòu)造的影響，也包括了地殼下界面起伏的影響等。從原始資料(圖8)可以看出，中部區(qū)域以及東南角區(qū)域布格重力異常值高，對應(yīng)地形圖中地形高程相對較低，西北角區(qū)域布格重力異常值最低，對應(yīng)地形圖中地形高程相對較高，布格重力異常范圍為-115～-15 mGal，布格重力異常能與地形圖高程信息聯(lián)系起來。

對南嶺東部EGM2008布格重力異常數(shù)據(jù)進(jìn)行Curvelet變換處理，分解為4層Curvelet系數(shù)，將第1層曲波系數(shù)設(shè)為深部層，第2層曲波系數(shù)設(shè)為中部層，第3層曲波系數(shù)設(shè)為淺部層，第4層曲波系數(shù)設(shè)為噪聲干擾層，對深部層、中部層、淺部層和干擾層曲波系數(shù)進(jìn)行反Curvelet變換重構(gòu)，就可以得到深層、中層、淺層和噪聲干擾層重力異常平面(圖10)。

圖10 Curvelet變換重構(gòu)重力異常Fig.10 The reconstructed gravity anomaly map using Curvelet transform

綜合前人研究成果，如參考羅凡華南地區(qū)衛(wèi)星布格重力異常和華南地區(qū)僅反映莫霍面變化的衛(wèi)星布格重力異常[31]，結(jié)合圖10a可以認(rèn)為第1層Curvelet系數(shù)重構(gòu)重力異常圖近似為深部區(qū)域異常，消除了大部分局部異常的干擾。大體上中部區(qū)域和東南角區(qū)域重力異常高，對應(yīng)該區(qū)域地形高程值相對較小，西北角區(qū)域重力異常值低，對應(yīng)該區(qū)域地形高程值相對較大，能與地形圖高程信息聯(lián)系起來。重力異常范圍為-95～-32 mGal，紅色為重力異常高值，藍(lán)色為重力異常低值。

對于圖10b中，第2層系數(shù)重構(gòu)重力異常圖近似為中層深度局部異常體或構(gòu)造產(chǎn)生的重力異常。圖10c中，第3層系數(shù)重構(gòu)重力異常圖近似為淺層局部異常體或各種地質(zhì)構(gòu)造產(chǎn)生的重力異常。圖10d為高頻干擾異常。

將圖10c中第3層系數(shù)重構(gòu)重力異常圖與該區(qū)域地質(zhì)圖(圖11)進(jìn)行對比，對比圖見圖12，可以看到地質(zhì)圖中標(biāo)定的鎢礦與錫礦的巖體分布(深藍(lán)色和棕色圓圈)與圖10c中重力異常體的分布和走向大致相似，可以認(rèn)為圖10c主要反映了淺層密度不均勻礦體、花崗巖巖體邊界及各種地質(zhì)構(gòu)造等產(chǎn)生的局部重力異常，對淺地表礦體分布或地質(zhì)構(gòu)造等信息具有一定指示作用。

圖11 南嶺東部地質(zhì)簡圖Fig.11 Geological sketch map of eastern Nanling

圖12 南嶺東部地質(zhì)簡圖(a)與重力異常(b)對比Fig.12 Comparison of geological sketch map (a) and gravity anomaly map (b) of eastern Nanling

3 結(jié)論

1) Curvelet變換根據(jù)原始數(shù)據(jù)的大小將其分解為不同尺度的系數(shù)，各個尺度的系數(shù)又可以經(jīng)過處理后反Curvelet變換為重構(gòu)數(shù)據(jù)，并且各個尺度的系數(shù)相互之間可以疊加、增刪。

2) Curvelet變換閾值去噪中，噪聲信號為高頻細(xì)節(jié)信號，一般在分解重構(gòu)的高層系數(shù)中，對高層系數(shù)進(jìn)行閾值處理可以有效去掉噪聲的干擾。

3) Curvelet變換運(yùn)用于位場數(shù)據(jù)中，可以近似的認(rèn)為，低層系數(shù)重構(gòu)數(shù)據(jù)表示深層區(qū)域異常信號，中層系數(shù)重構(gòu)數(shù)據(jù)表示中層異常體產(chǎn)生的信號， 高層系數(shù)重構(gòu)圖表示淺層局部異常體信號。由南嶺實(shí)測數(shù)據(jù)處理看出，低層Curvelet系數(shù)重構(gòu)重力異?？梢越普J(rèn)為深部區(qū)域重力異常，高層Curvelet系數(shù)重構(gòu)重力異?？梢越普J(rèn)為淺層密度不均勻礦體、花崗巖巖體邊界及各種地質(zhì)構(gòu)造等產(chǎn)生的局部重力異常，對淺地表地質(zhì)體或各種地質(zhì)構(gòu)造信息具有一定指示作用。

致謝：感謝審稿人和編輯對本文提出的寶貴意見。感謝中國地質(zhì)大學(xué)(北京)同學(xué)和深圳市地質(zhì)環(huán)境監(jiān)測中心的幫助。