重力異常分離的小波域優(yōu)化位變?yōu)V波方法

劉彩云， 姚長(zhǎng)利， 鄭元滿

1 長(zhǎng)江大學(xué) 信息與數(shù)學(xué)學(xué)院， 湖北荊州 4340232 中國(guó)地質(zhì)大學(xué) 地球物理與信息技術(shù)學(xué)院， 北京 100083

?

重力異常分離的小波域優(yōu)化位變?yōu)V波方法

劉彩云1, 2， 姚長(zhǎng)利2*， 鄭元滿2

1 長(zhǎng)江大學(xué) 信息與數(shù)學(xué)學(xué)院， 湖北荊州 4340232 中國(guó)地質(zhì)大學(xué) 地球物理與信息技術(shù)學(xué)院， 北京 100083

在重力異常分離中，頻率域?yàn)V波分離方法是以全局?jǐn)?shù)據(jù)頻譜特征設(shè)計(jì)針對(duì)性的濾波器實(shí)現(xiàn)的.濾波器參數(shù)與空間位置無(wú)關(guān)，因此無(wú)法針對(duì)局部數(shù)據(jù)頻譜特征動(dòng)態(tài)調(diào)整濾波器參數(shù).針對(duì)該局限性，在小波域?yàn)V波理論和優(yōu)化濾波方法的基礎(chǔ)上，本文研究提出了小波域優(yōu)化位變?yōu)V波法，該方法具有空間變化濾波能力，在不同空間位置實(shí)現(xiàn)不同的濾波器特性，從而能實(shí)現(xiàn)局部數(shù)據(jù)頻譜與全局?jǐn)?shù)據(jù)頻譜存在較大差異的重力異常分離問(wèn)題.理論模型數(shù)據(jù)分離實(shí)驗(yàn)結(jié)果表明，在局部數(shù)據(jù)頻譜與全局?jǐn)?shù)據(jù)頻譜差異較大的情況下，該方法相對(duì)于Butterworth濾波和優(yōu)化濾波等方法具有優(yōu)勢(shì).最后，用一個(gè)實(shí)例進(jìn)行檢驗(yàn)計(jì)算，體現(xiàn)了所提方法技術(shù)的效果和應(yīng)用前景.

異常分離； 小波域； 優(yōu)化濾波； 位變?yōu)V波； 重力異常

1 引言

實(shí)測(cè)重力異常是不同埋深、不同規(guī)模、不同形態(tài)的所有密度不均勻地質(zhì)體重力效應(yīng)的疊加，若要由實(shí)測(cè)重力異常反演某個(gè)特定地質(zhì)體，必須先從疊加的重力異常中分離出單純由該目標(biāo)地質(zhì)體引起的異常，再用分離出的異常進(jìn)行反演和解釋.由于重力場(chǎng)具有固有的多解性，從疊加異常中分離出目標(biāo)地質(zhì)體的單一異常是困難的.重力異常的有效分離是目前重力資料處理解釋中的難題之一.

重力異常分離的方法很多，現(xiàn)在常用的主要包括最小二乘平滑法、解析延拓法、趨勢(shì)分析法、匹配濾波法、補(bǔ)償圓滑濾波法、正則化濾波法、維納濾波法、小波變換法、非線性濾波法、優(yōu)化濾波法等（Spector and Grant, 1970; 侯重初, 1981; 安玉林和管志寧, 1985; Gupta and Ramani, 1980; Pawlowski and Hansen, 1990; Pawlowski, 1994; Fedi and Quarta, 1998; Ridsdill-Smith, 1998; 楊文采等, 2001; Keating and Pinet, 2011; 郭良輝等2012）.這些方法具有各自的特點(diǎn)與優(yōu)勢(shì)，但同時(shí)也存在各自的局限性.

例如，維納濾波法中，為了準(zhǔn)確地分離出目標(biāo)地質(zhì)體引起的異常，需要先估計(jì)出目標(biāo)地質(zhì)體的頻譜，再以此為依據(jù)設(shè)計(jì)維納濾波器.針對(duì)這個(gè)問(wèn)題，Gupta和Ramani（1980）采用徑向譜分析方法，估計(jì)目標(biāo)地質(zhì)體的頻譜；Pawlowski和Hansen（1990）利用鉆孔數(shù)據(jù)等先驗(yàn)地質(zhì)信息構(gòu)建模型，估計(jì)目標(biāo)地質(zhì)體的頻譜，此后，Pawlowski（1994）又提出利用格林等效層理論估計(jì)目標(biāo)地質(zhì)體頻譜的方法；郭良輝等（2012）在優(yōu)選延拓的理論基礎(chǔ)上，進(jìn)一步研究提出基于維納濾波和格林等效層理論的優(yōu)化濾波法，該方法分離重力異常不需要已知延拓高度，可對(duì)重力異常進(jìn)行指定頻段的低通帶通和高通濾波.這些經(jīng)典的頻率域分離方法都是根據(jù)全局?jǐn)?shù)據(jù)頻譜特征設(shè)計(jì)濾波器進(jìn)行濾波分離，對(duì)于局部數(shù)據(jù)頻譜和全局?jǐn)?shù)據(jù)頻譜差異較大的重力異常分離，頻率域分離方法的效果有待改進(jìn).熊光楚（1979）、熊光楚和張福榮（1981）曾針對(duì)線性濾波器局限性問(wèn)題，研究并提出用于重磁異常解釋的空間域位變?yōu)V波器.只是當(dāng)時(shí)考慮的還是異常組合不太復(fù)雜的簡(jiǎn)單情況，還不適合大數(shù)據(jù)量及復(fù)雜異常的分離情況.

小波變換是近年來(lái)獲得廣泛應(yīng)用的一種數(shù)字信號(hào)處理工具，小波變換也可用于估計(jì)信號(hào)的局部頻譜特性.Kaiser（1996）經(jīng)過(guò)理論推導(dǎo)指出，對(duì)信號(hào)進(jìn)行頻率域?yàn)V波操作，可通過(guò)對(duì)信號(hào)的小波變換系數(shù)進(jìn)行線性組合運(yùn)算來(lái)近似實(shí)現(xiàn)，并稱(chēng)之為小波域?yàn)V波.本文針對(duì)常規(guī)頻率域?yàn)V波方法存在的不具有位變?yōu)V波的局限性問(wèn)題，在小波域?yàn)V波理論和優(yōu)化濾波方法的基礎(chǔ)上，研究提出小波域優(yōu)化位變?yōu)V波方法，該方法充分發(fā)揮小波變換的空間定位優(yōu)點(diǎn)，能在不同空間位置表現(xiàn)出不同的濾波器特性，適用于局部數(shù)據(jù)頻譜與全局?jǐn)?shù)據(jù)頻譜存在較大差異的重力異常分離問(wèn)題.

本文首先簡(jiǎn)要闡述小波域?yàn)V波方法，然后給出小波域位變?yōu)V波器設(shè)計(jì)方法和小波域優(yōu)化位變?yōu)V波方法及其計(jì)算步驟，再利用該方法對(duì)理論模型數(shù)據(jù)進(jìn)行異常分離實(shí)驗(yàn)，并與Butterworth濾波（李鐘慎, 2006）和頻率域優(yōu)化濾波法（郭良輝等, 2012）的處理結(jié)果進(jìn)行對(duì)比，最后，用一個(gè)實(shí)例檢驗(yàn)了所提方法技術(shù)的效果.

2 小波域優(yōu)化位變?yōu)V波

2.1 小波域?yàn)V波原理

小波域?yàn)V波是基于小波變換的濾波方法，Kaiser（1996）給出了詳細(xì)的理論推導(dǎo)，這里先做簡(jiǎn)要闡述.

信號(hào)s（x）的連續(xù)小波變換（ContinuousWaveletTransform,CWT）定義為：

（1）

令

（2）

公式（1）可用卷積方式描述如下：

（3）

對(duì)公式（3）兩邊做傅里葉變換，可得：

（4）

為便于后面的描述，簡(jiǎn)記為：

（5）

為實(shí)現(xiàn)濾波功能，將公式（5）兩邊同乘以向量v，v=（v1,v2,…,vk），其第i個(gè)分量vi是第i個(gè)尺度ai的函數(shù)vi=v（ai）（i=1,2,…,k），vi稱(chēng)為濾波器的等效系數(shù).可得：

（6）

（7）

則（6）式簡(jiǎn)寫(xiě)為

（8）

對(duì)公式（8）兩邊做傅里葉反變換，可得到：

（9）

公式（9）表示，信號(hào)s′（x）可按如下步驟得到：先對(duì)原信號(hào)s（x）做小波變換，然后對(duì)小波變換系數(shù)W（ai,b）加權(quán)求和，權(quán)重系數(shù)即為小波域?yàn)V波器等效系數(shù)vi.

2.2 小波域位變?yōu)V波

本文將Kaiser（1996）給出的小波域?yàn)V波器的構(gòu)建方法進(jìn)行推廣，構(gòu)建小波域位變?yōu)V波器.具體構(gòu)建過(guò)程如下：

假設(shè)公式（6）中濾波器等效系數(shù)vi不僅與尺度ai有關(guān)，而且與空間位置b有關(guān)，即為函數(shù)v（ai,b）.則公式（6）可寫(xiě)作：

（10）

（11）

公式（10）簡(jiǎn)寫(xiě)為：

（12）

對(duì)公式（12）兩邊做傅里葉反變換，可得到：

（13）

將公式（11）展開(kāi)，寫(xiě)成如下的方程組形式：

（14）

值得指出的是，求解方程組（14）得到的v（ai,b）是該方程組的最小二乘解，因此，不可避免的會(huì)引入誤差，為后面敘述方便，本文稱(chēng)之為小波域位變?yōu)V波誤差.

小波域位變?yōu)V波仿真實(shí)驗(yàn)結(jié)果見(jiàn)3.1節(jié).

2.3 小波域優(yōu)化位變?yōu)V波

基于小波變換進(jìn)行優(yōu)化位變?yōu)V波計(jì)算的主要步驟如下：

（1） 根據(jù)實(shí)測(cè)異常數(shù)據(jù)的區(qū)域特征和數(shù)據(jù)處理的需要，將實(shí)測(cè)異常數(shù)據(jù)進(jìn)行分塊（對(duì)數(shù)據(jù)做小波變換時(shí)頻分析，確定各頻率成分在空間的分布規(guī)律），相鄰分塊在邊界處可相互重疊或不重疊，本文選擇相鄰分塊間不重疊.

（2） 對(duì)每塊實(shí)測(cè)數(shù)據(jù)進(jìn)行頻率域優(yōu)化濾波法的特征評(píng)價(jià)，設(shè)計(jì)從第j塊異常數(shù)據(jù)中分離出第l層數(shù)據(jù)的優(yōu)化濾波轉(zhuǎn)換函數(shù)Hjl（r）.

i） 進(jìn)行傅里葉變換，計(jì)算第j塊異常數(shù)據(jù)的對(duì)數(shù)功率譜Pj（r）.

ii） 分段線性擬合第j塊異常數(shù)據(jù)的對(duì)數(shù)功率譜.

分析第j塊異常數(shù)據(jù)的對(duì)數(shù)功率譜，確定頻段數(shù)和各頻段的頻率范圍，每一頻段對(duì)應(yīng)一個(gè)格林等效層.利用直線Pjl（r）=Kjlr+Bjl擬合每一頻段異常數(shù)據(jù)對(duì)數(shù)功率譜，其中，Kjl,Bjl分別是擬合第j塊異常數(shù)據(jù)的第l頻段的直線的斜率和截距.

iii） 用格林等效層對(duì)數(shù)功率譜擬合第j塊異常數(shù)據(jù)實(shí)測(cè)對(duì)數(shù)功率譜Pj（r）.

iv） 假設(shè)目標(biāo)層與其他深度層場(chǎng)源信息互不相關(guān)，則可設(shè)計(jì)從第j塊異常數(shù)據(jù)中分離出第l層數(shù)據(jù)的優(yōu)化濾波轉(zhuǎn)換函數(shù).具體為

（15）

其中，Pj（r）是第j塊異常數(shù)據(jù)的功率譜密度函數(shù)，Pjl（r）是第j塊異常數(shù)據(jù)第l層（目標(biāo)層）異常數(shù)據(jù)的功率譜密度函數(shù).

（4） 將每個(gè)分塊異常數(shù)據(jù)計(jì)算得到的局部濾波等效系數(shù)向量v（bj）根據(jù)其空間位置進(jìn)行組合，得到全體重力異常的小波域位變?yōu)V波等效系數(shù)向量v（b）.

（5） 對(duì)全體重力異常數(shù)據(jù)做小波變換，得到小波變換系數(shù)W（a,b）.

（6） 將小波域位變?yōu)V波器等效系數(shù)向量v（b）和W（a,b）代入公式（13），計(jì)算得到小波域優(yōu)化位變?yōu)V波分離出的目標(biāo)層（第l層）異常數(shù)據(jù).

（7） 算法結(jié)束.

基于此，我們就建立了小波域優(yōu)化位變?yōu)V波方法.下面通過(guò)理論數(shù)據(jù)分離實(shí)驗(yàn)和實(shí)測(cè)資料處理，測(cè)試其效果.

3 理論數(shù)據(jù)分離實(shí)驗(yàn)

3.1 小波域位變?yōu)V波實(shí)驗(yàn)

實(shí)驗(yàn)1：小波域位變?yōu)V波器的特點(diǎn)是濾波器參數(shù)可隨空間位置的變化而變化，下面以一維信號(hào)為例，說(shuō)明小波域位變?yōu)V波的功能.給定如下信號(hào)：

（16）其中，f1=0.005 cycles·m-1，f2=0.01 cycles·m-1，f3=0.04 cycles·m-1，f4=0.08 cycles·m-1.實(shí)驗(yàn)1信號(hào)波形及其小波變換時(shí)頻分析如圖1所示.

從圖1a可看出，該信號(hào)含有四個(gè)不同的頻率成分，且具有明顯的位置特征（x≤500 m處是頻率為0.005 cycles·m-1和0.01 cycles·m-1的信號(hào)，x>500 m處是頻率為0.04 cycles·m-1和0.08 cycles·m-1的信號(hào)）.從圖1b所示的小波變換時(shí)頻分析圖（為無(wú)量綱的小波系數(shù)）可以看出，以x=500 m處為界，時(shí)頻分析圖的左右兩半部分具有明顯不同的頻率成分.其中，x≤500 m處具有兩段大尺度的能量異常，表明x≤500 m區(qū)域有兩個(gè)低頻信號(hào)；x>500 m處具有兩段小尺度的能量異常，表明x>500 m處有兩個(gè)高頻信號(hào).

對(duì)實(shí)驗(yàn)1信號(hào)分別采用傳統(tǒng)的頻率域?yàn)V波和小波域位變?yōu)V波方法進(jìn)行對(duì)比實(shí)驗(yàn)，結(jié)果如圖2—4所示.其中，采用截止頻率為0.008 cycles·m-1，階數(shù)為16的Butterworth（李鐘慎，2006）低通/高通濾波器，濾波結(jié)果如圖2所示.

采用截止頻率為0.06 cycles·m-1，階數(shù)為16的Butterworth低通/高通濾波器，濾波結(jié)果如圖3所示.

采用小波域位變?yōu)V波器，該濾波器在x≤500 m的地方截止頻率為0.008 cycles·m-1，在x>500 m的地方截止頻率為0.06 cycles·m-1，小波域位變?yōu)V波器濾波結(jié)果如圖4所示.

從圖2—4可以看出：

（1） 采用頻率域?yàn)V波分離時(shí)，若設(shè)截止頻率為0.008 cycles·m-1，則能分離出x≤500 m處的兩個(gè)信號(hào)：頻率分別為0.005 cycles·m-1和0.01 cycles·m-1，但無(wú)法分離x>500 m處的兩個(gè)信號(hào)：頻率分別為0.04 cycles·m-1和0.08 cycles·m-1（見(jiàn)圖2）；

（2） 若設(shè)截止頻率為0.06 cycles·m-1，則能分離出x>500 m處的兩個(gè)信號(hào)：頻率分別為0.04 cycles·m-1和0.08 cycles·m-1，但無(wú)法分離x≤500 m處的兩個(gè)信號(hào)：頻率分別為0.005 cycles·m-1和0.01 cycles·m-1（見(jiàn)圖3）；

圖1 實(shí)驗(yàn)1信號(hào)及其小波變換時(shí)頻分析圖

圖2 截止頻率0.008 cycles·m-1的Butterworth低通/高通濾波結(jié)果

圖3 截止頻率0.06 cycles·m-1的Butterworth低通/高通濾波結(jié)果

（3） 采用小波域位變?yōu)V波分離時(shí)，設(shè)x≤500 m處截止頻率為0.008 cycles·m-1，x>500 m處截止頻率為0.06 cycles·m-1，能同時(shí)分離出四個(gè)頻率的信號(hào)（見(jiàn)圖4）.

圖4 小波域位變?yōu)V波結(jié)果（x≤500 m處截止頻率0.008 cycles·m-1，x>500 m處截止頻率0.06 cycles·m-1）

小波變換這種能根據(jù)信號(hào)不同位置的特點(diǎn)，設(shè)計(jì)針對(duì)性的適應(yīng)位置變化的位變?yōu)V波功能，顯然是過(guò)去傳統(tǒng)的重磁處理濾波方法手段所不具有的.經(jīng)典的頻率域?yàn)V波實(shí)現(xiàn)重磁異常分離方法是根據(jù)全局?jǐn)?shù)據(jù)的頻譜特征設(shè)計(jì)濾波器進(jìn)行濾波分離的，當(dāng)局部數(shù)據(jù)頻譜與全局?jǐn)?shù)據(jù)頻譜差異較大時(shí)，頻率域分離方法會(huì)遇到和圖2、圖3類(lèi)似的困難，難以進(jìn)行精細(xì)定位分離.小波域位變?yōu)V波方法能針對(duì)不同局部數(shù)據(jù)采用不同濾波器，能有效克服分離不精細(xì)的缺陷.

3.2 理論模型數(shù)據(jù)實(shí)驗(yàn)

實(shí)驗(yàn)2：為了檢驗(yàn)小波域優(yōu)化位變?yōu)V波的效果，我們?cè)O(shè)計(jì)了如下組合模型，該組合模型由四個(gè)水平無(wú)限延伸的長(zhǎng)方體構(gòu)成，各模型參數(shù)如表1所示.在地表測(cè)點(diǎn)間距20 m，剖面長(zhǎng)度為10000 m，重力異常及模型分布剖面如圖5所示.

圖5為組合模型的重力異常.從異常分布可以看出，盡管模型1、2的異常幅值比模型3、4的異常要小，但其異常寬度要更窄，如果簡(jiǎn)單地基于頻譜分析進(jìn)行常規(guī)譜分離，則只能將異常體1、2的異常分離出來(lái)，而無(wú)法分離出相對(duì)較淺的異常體1、3的異常，而這方面正是小波域優(yōu)化位變?yōu)V波能體現(xiàn)的特點(diǎn).

表1 實(shí)驗(yàn)2各模型參數(shù)表

圖5 實(shí)驗(yàn)2重力異常及其剖面圖

下面分別采用頻率域優(yōu)化濾波和小波域優(yōu)化位變?yōu)V波方法進(jìn)行重力異常分離對(duì)比實(shí)驗(yàn).功率譜分段擬合分析對(duì)頻率域優(yōu)化濾波分離效果起到重要作用，功率譜分段擬合分得越細(xì)，優(yōu)化濾波分離結(jié)果越好.為了客觀對(duì)比頻率域優(yōu)化濾波方法和小波域優(yōu)化位變?yōu)V波方法，我們?cè)趦煞N方法中均采用低頻、中頻、高頻三個(gè)頻段分段擬合，分別對(duì)應(yīng)區(qū)域異常、局部異常和擬合誤差（其含義后面闡述）.模型實(shí)驗(yàn)2重力異常功率譜分析如圖6所示.

圖6 頻率域優(yōu)化濾波法功率譜分析與擬合結(jié)果

以此為指導(dǎo)，設(shè)計(jì)頻率域優(yōu)化濾波器，模型實(shí)驗(yàn)2的頻率域優(yōu)化濾波法分離結(jié)果如圖7所示.其中，低通濾波分離結(jié)果為深部區(qū)域異常，帶通濾波分離結(jié)果為淺部局部異常，高通濾波分離結(jié)果為優(yōu)化濾波中采用格林等效層近似實(shí)測(cè)數(shù)據(jù)所產(chǎn)生的擬合誤差（郭良輝等, 2012）.

圖7 實(shí)驗(yàn)2重力異常頻率域優(yōu)化濾波分離結(jié)果

從圖7可以看出，分離出的區(qū)域異常（圖7a）主要是長(zhǎng)方體2、3、4的異常，同時(shí)含有少量的長(zhǎng)方體1的異常；分離出的局部異常（圖7b）主要是長(zhǎng)方體1的異常，同時(shí)含有極少量的長(zhǎng)方體3的異常.值得注意的是，在區(qū)域異常中（圖7a），長(zhǎng)方體3、4的異常疊加在一起，未能有效分離開(kāi).如果單從深淺場(chǎng)源分離的角度看（即將淺部場(chǎng)源1、3從相對(duì)較深的場(chǎng)源2、4中分離出來(lái)），這個(gè)效果也不好，同時(shí)也不可能效果好.因?yàn)楸容^而言，局部深的場(chǎng)源2（相對(duì)于場(chǎng)源1）與局部淺的場(chǎng)源3（相對(duì)于場(chǎng)源4）相比，場(chǎng)源2的頻譜反而更高，常規(guī)的整體頻譜分離方式，不可能將場(chǎng)源2歸算為深部場(chǎng)源，以及將場(chǎng)源3歸算為淺部場(chǎng)源.小波位變?yōu)V波則具有這樣的能力.

對(duì)模型實(shí)驗(yàn)2的重力異常數(shù)據(jù)進(jìn)行小波變換時(shí)頻分析，如圖8所示.這里需要說(shuō)明的是，圖8中的大尺度對(duì)應(yīng)于低頻異常成分，反映的是深部場(chǎng)源信息，為了尺度與深度直觀對(duì)應(yīng)，圖8中縱坐標(biāo)選擇從大到小.

從圖8可以看出，大致以x=3500 m為界，圖的左右部分具有不同的頻率成分，其中，x>3500 m處具有明顯強(qiáng)的低頻成分和中頻成分，而x<3500 m處，低頻成分強(qiáng)度較小，且具有頻率更高的高頻成分（小波系數(shù)延伸到更小尺度處）.以此為依據(jù)，將實(shí)驗(yàn)1的模型重力異常以x=3500 m為界分為前后兩段局部數(shù)據(jù)，分別對(duì)局部數(shù)據(jù)進(jìn)行功率譜分析（如圖9所示），兩段局部數(shù)據(jù)的功率譜與全局?jǐn)?shù)據(jù)的功率譜存在明顯的差異.如果能根據(jù)每段局部數(shù)據(jù)的功率譜特征，有針對(duì)性的設(shè)計(jì)適用于該段局部數(shù)據(jù)的濾波器，則有望改善分離效果.

分別對(duì)分段后的局部數(shù)據(jù)進(jìn)行功率譜分析和分段線性擬合，結(jié)果如圖10所示.

以圖10結(jié)果為依據(jù)，按照本文2.3節(jié)所述的方法設(shè)計(jì)小波域優(yōu)化位變?yōu)V波器，模型重力異常分離結(jié)果如圖11所示.其中，位變低通、帶通濾波分離結(jié)果為深部區(qū)域異常和淺部局部異常，位變高通濾波分離結(jié)果為擬合誤差，由優(yōu)化濾波擬合誤差（郭良輝等, 2012）、小波域位變?yōu)V波誤差和分塊頻譜估計(jì)所引入誤差的疊加而成.

從圖11所示的小波優(yōu)化位變?yōu)V波分離結(jié)果可以看出，僅需采用一次小波域位變?yōu)V波，即可有效分離出四個(gè)長(zhǎng)方體的異常.其中，分離出的區(qū)域異常包含長(zhǎng)方體2、4的異常，局部異常包含長(zhǎng)方體1、3的異常，與給定的模型符合性明顯改進(jìn).

值得指出的是，圖7c中優(yōu)化濾波高通濾波結(jié)果優(yōu)于圖11c中小波優(yōu)化位變?yōu)V波高通濾波結(jié)果，這是因?yàn)閳D11c中除了含有與圖7c相同的優(yōu)化濾波擬合誤差（郭良輝等, 2012）之外，還含有小波域位變?yōu)V波誤差和分塊頻譜估計(jì)引入的誤差.小波域優(yōu)化位變?yōu)V波中，高通濾波結(jié)果的幅值明顯小于低通和帶通濾波結(jié)果的幅值，若我們重點(diǎn)關(guān)注區(qū)域異常/局部異常分離效果的改善，則該方法對(duì)擬合誤差的小幅放大或可接受.另外，小波域優(yōu)化位變?yōu)V波的擬合誤差與選用的母小波有關(guān)，可以通過(guò)選用最佳的母小波等方法減小擬合誤差，由于篇幅限制，將另文詳細(xì)描述.

圖8 實(shí)驗(yàn)2重力異常數(shù)據(jù)小波時(shí)頻分析圖

圖9 全局?jǐn)?shù)據(jù)與局部數(shù)據(jù)功率譜分析對(duì)比圖

圖10 實(shí)驗(yàn)2重力異常局部數(shù)據(jù)功率譜分析結(jié)果圖

圖11 實(shí)驗(yàn)2重力異常小波域優(yōu)化位變?yōu)V波分離結(jié)果

對(duì)小波域優(yōu)化位變?yōu)V波和頻率域優(yōu)化濾波分離出的區(qū)域異常與局部異常，分別計(jì)算其平均絕對(duì)誤差和標(biāo)準(zhǔn)差，結(jié)果如表2所示.從表2可以看出，小波域優(yōu)化位變?yōu)V波方法相對(duì)于頻率域優(yōu)化濾波方法能更好分離出區(qū)域異常和局部異常，前者分離結(jié)果優(yōu)于后者.

表2 實(shí)驗(yàn)2分離結(jié)果對(duì)比分析

實(shí)驗(yàn)2通過(guò)相對(duì)簡(jiǎn)單的模型試驗(yàn)，驗(yàn)證了小波域優(yōu)化濾波方法的有效性.下面我們?cè)O(shè)計(jì)相對(duì)復(fù)雜的組合模型，檢驗(yàn)頻率混疊嚴(yán)重情況下小波域優(yōu)化位變?yōu)V波的效果.

實(shí)驗(yàn)3：組合模型由七個(gè)水平無(wú)限延伸的長(zhǎng)方體構(gòu)成，各模型參數(shù)如表3所示.在地表測(cè)點(diǎn)間距20 m，剖面長(zhǎng)度為10000 m，重力異常及模型分布剖面如圖12所示.

表3 實(shí)驗(yàn)3各模型參數(shù)表

圖12 實(shí)驗(yàn)3重力異常及其剖面圖

從圖12和表3可以看出，長(zhǎng)方體C1從淺部向下延伸到深部，其頂深與A1、A2的相同，底深大于B1、B2的底深，可知其頻譜會(huì)與A組、B組長(zhǎng)方體的頻譜重疊.對(duì)實(shí)驗(yàn)3重力異常數(shù)據(jù)進(jìn)行功率譜分析，結(jié)果如圖13所示，圖中A、B、C、D四組長(zhǎng)方體的功率譜相互重疊嚴(yán)重.

圖13 實(shí)驗(yàn)3重力異常功率譜圖

從圖13可看出，模型實(shí)驗(yàn)3重力異常功率譜分為低頻、中頻、高頻三個(gè)頻段.采用頻率域優(yōu)化濾波方法進(jìn)行線性擬合，其中，低頻段波數(shù)為[0,0.0007]cycles·m-1，對(duì)應(yīng)場(chǎng)源平均深度為840.5 m；中頻段波數(shù)為[0.0007,0.0049]cycles·m-1，對(duì)應(yīng)場(chǎng)源平均深度為141.4 m；高頻段波數(shù)為[0.0049,0.025]cycles·m-1，對(duì)應(yīng)場(chǎng)源平均深度為12.0 m.以此為指導(dǎo)，設(shè)計(jì)頻率域優(yōu)化濾波器.

模型實(shí)驗(yàn)3的頻率域優(yōu)化濾波法分離結(jié)果如圖14所示.

圖14 實(shí)驗(yàn)3重力異常頻率域優(yōu)化濾波分離結(jié)果

從圖14可以看出，分離出的區(qū)域異常（圖14a）主要是B組和D組長(zhǎng)方體的異常，同時(shí)含有少量的A組和C組長(zhǎng)方體的異常；分離出的局部異常（圖14b）主要是A組和C組長(zhǎng)方體的異常.值得注意的是，在分離出的區(qū)域異常中（圖14a），左側(cè)呈單峰特征，不能分辨出兩個(gè)B組的長(zhǎng)方體，右側(cè)由于殘留的C組長(zhǎng)方體的異常較多，曲線形態(tài)變形較大；相應(yīng)的，在分離出的局部異常中（圖14b），左側(cè)曲線形態(tài)變化較大，右側(cè)C組長(zhǎng)方體的異常幅值衰減較大.

對(duì)模型實(shí)驗(yàn)3的重力異常數(shù)據(jù)進(jìn)行小波變換時(shí)頻分析，根據(jù)時(shí)頻分析結(jié)果，將重力異常以x=4000 m為界分為前后兩段局部數(shù)據(jù).分別對(duì)兩段局部數(shù)據(jù)進(jìn)行功率譜分析，可識(shí)別出四個(gè)等效場(chǎng)源，平均場(chǎng)源深度分別為64.5 m、302.0 m、263.6 m、1214.7 m，與理論模型的A、B、C、D四組長(zhǎng)方體相對(duì)應(yīng).以此為指導(dǎo)，設(shè)計(jì)小波域優(yōu)化位變?yōu)V波器，模型重力異常分離結(jié)果如圖15所示.

圖15 實(shí)驗(yàn)3重力異常小波域優(yōu)化位變?yōu)V波分離結(jié)果

從圖15可以看出，分離出的區(qū)域異常（圖15a）主要是B組和D組長(zhǎng)方體的異常，局部異常（圖15b）主要是A組和C組長(zhǎng)方體的異常.在分離出的區(qū)域異常中（圖15a），左側(cè)呈雙峰特征，能夠很好地分辨出B1、B2長(zhǎng)方體，右側(cè)的曲線形態(tài)也與理論值一致；在分離出的局部異常（圖15b）中，左側(cè)A組長(zhǎng)方體的異常曲線形態(tài)與理論值一致，僅幅值有所衰減，右側(cè)C組長(zhǎng)方體的異常曲線形態(tài)也與理論值相一致，且幅值衰減得到明顯的改善.對(duì)比圖14和圖15可以看出，小波優(yōu)化位變?yōu)V波的區(qū)域異常、局部異常分離結(jié)果明顯優(yōu)于頻率域優(yōu)化濾波.

分別計(jì)算兩種方法分離結(jié)果的平均絕對(duì)誤差和標(biāo)準(zhǔn)差，結(jié)果如表4所示.從表4可以看出，小波域優(yōu)化位變?yōu)V波方法相對(duì)于頻率域優(yōu)化濾波方法能更

好分離出區(qū)域異常和局部異常，前者分離結(jié)果優(yōu)于后者.

表4 實(shí)驗(yàn)3分離結(jié)果對(duì)比分析

從實(shí)驗(yàn)2和實(shí)驗(yàn)3重力異常分離實(shí)驗(yàn)結(jié)果可以看出，當(dāng)局部數(shù)據(jù)頻譜特征和全局?jǐn)?shù)據(jù)頻譜特征差異較大時(shí)，全局?jǐn)?shù)據(jù)的頻譜是所有場(chǎng)源信息綜合影響的結(jié)果，經(jīng)典的頻率域分離方法會(huì)出現(xiàn)分離不徹底的問(wèn)題.這種情況下，小波域優(yōu)化位變?yōu)V波法具有一定優(yōu)勢(shì)，能根據(jù)局部數(shù)據(jù)的頻譜特征動(dòng)態(tài)調(diào)整濾波器參數(shù)，有效改善分離效果.

另一方面，當(dāng)局部數(shù)據(jù)頻譜與全局?jǐn)?shù)據(jù)頻譜差異很小時(shí)，小波域優(yōu)化位變?yōu)V波方法將近似地退化為小波域優(yōu)化濾波方法，由于小波域?yàn)V波誤差的存在，其分離效果會(huì)略差于頻率域優(yōu)化濾波方法.

4 實(shí)測(cè)資料處理

實(shí)驗(yàn)4為某一重力異常剖面實(shí)測(cè)資料，在地表測(cè)點(diǎn)間距0.25 km，剖面長(zhǎng)度為380 km，剖面重力異常如圖16所示.

對(duì)實(shí)測(cè)資料數(shù)據(jù)分別用頻率域優(yōu)化濾波、Butterworth濾波和小波域優(yōu)化位變?yōu)V波進(jìn)行異常分離.

（1）頻率域優(yōu)化濾波

對(duì)實(shí)測(cè)剖面重力異常進(jìn)行功率譜分析，結(jié)果如圖17所示.

以此為指導(dǎo)，設(shè)計(jì)優(yōu)化濾波器，頻率域優(yōu)化濾波法分離結(jié)果如圖18所示.

圖16 實(shí)測(cè)數(shù)據(jù)剖面重力異常

圖17 實(shí)測(cè)剖面重力異常功率譜分析與優(yōu)化濾波擬合

圖18 實(shí)測(cè)剖面重力異常頻率域優(yōu)化濾波異常分離結(jié)果

（2）Butterworth濾波

作為對(duì)比，我們根據(jù)實(shí)測(cè)剖面重力異常功率譜分析結(jié)果（圖17），設(shè)計(jì)Butterworth濾波器對(duì)實(shí)測(cè)剖面重力異常進(jìn)行傳統(tǒng)的頻率域?yàn)V波分離，其中，低通截止頻率為0.0053 cycles·km-1，高通截止頻率為0.3 cycles·km-1，濾波器階數(shù)為16，分離結(jié)果如圖19所示.

從圖18和圖19可以看出，Butterworth濾波和頻率域優(yōu)化濾波分離結(jié)果相似，區(qū)域異常形態(tài)簡(jiǎn)單，僅反映出該剖面異常具有兩邊高、中間低的大尺度（波長(zhǎng)大于180 km）特征；相應(yīng)的，局部異常形態(tài)過(guò)于復(fù)雜，包含了波長(zhǎng)從6 km到180 km的所有異常，局部異常分離結(jié)果不夠理想.

（3）小波域優(yōu)化位變?yōu)V波

對(duì)實(shí)測(cè)剖面重力異常進(jìn)行小波變換時(shí)頻分析，如圖20所示.

圖19 實(shí)測(cè)剖面重力異常頻率域Butterworth濾波異常分離結(jié)果

圖20 實(shí)測(cè)剖面重力異常小波變換時(shí)頻分析圖

圖21 全局?jǐn)?shù)據(jù)與局部數(shù)據(jù)功率譜分析對(duì)比圖

對(duì)分段后的實(shí)測(cè)剖面重力異常分別進(jìn)行功率譜分析，并分段線性擬合，按照本文2.3節(jié)所述的方法設(shè)計(jì)小波域優(yōu)化位變?yōu)V波器.實(shí)測(cè)剖面重力異常分離結(jié)果如圖22所示.

從圖22可以看出，小波域優(yōu)化位變?yōu)V波分離結(jié)果中，區(qū)域異常信息豐富，不僅反映出該剖面異常所具有的兩邊高、中間低的大尺度特征，還包含較多的中等尺度異常；相應(yīng)的，局部異常則分辨得更加清晰，僅含有小尺度的淺部場(chǎng)源異常.

圖23是實(shí)測(cè)剖面重力異常密度反演結(jié)果圖.對(duì)比圖18、圖19和圖23則會(huì)發(fā)現(xiàn)，頻率域Butterworth濾波和優(yōu)化濾波分離出的區(qū)域異常過(guò)于簡(jiǎn)單，與反演結(jié)果不能很好的對(duì)應(yīng)，而局部異常中殘留有較多區(qū)域異常，既包含深度小于10 km的淺部場(chǎng)源異常，又包含深度10～50 km的深部場(chǎng)源異常.

對(duì)比圖22和圖23可以看出，小波域優(yōu)化位變?yōu)V波分離出的區(qū)域異常與反演結(jié)果中深度大于10 km的場(chǎng)源分布規(guī)律一致，局部異常與反演結(jié)果中深度小于10 km的場(chǎng)源分布規(guī)律一致，小波域優(yōu)化位變?yōu)V波的分離結(jié)果與反演解釋結(jié)果吻合得相當(dāng)好.

以上三種方法均對(duì)實(shí)測(cè)剖面重力異常分高、中、低三個(gè)頻段進(jìn)行分離，以便客觀對(duì)比三種方法的分離效果，實(shí)驗(yàn)結(jié)果體現(xiàn)了小波域優(yōu)化位變?yōu)V波法的優(yōu)越性.值得指出的是，由于采用的頻段數(shù)較少，分離結(jié)果仍可進(jìn)一步進(jìn)行分解，例如，小波域優(yōu)化位變?yōu)V波分離出的區(qū)域異常（見(jiàn)圖22a）仍然含有更深場(chǎng)源對(duì)應(yīng)的低頻成分.再次采用小波域優(yōu)化位變?yōu)V波方法對(duì)其進(jìn)行濾波分離，得到結(jié)果如圖24所示.

圖22 小波域優(yōu)化位變?yōu)V波分離結(jié)果

圖23 實(shí)測(cè)剖面重力異常反演解釋成果圖

由圖24可以看出，小波域優(yōu)化位變?yōu)V波所得區(qū)域異常的二次低通濾波結(jié)果與深度≥45 km的更深場(chǎng)源相對(duì)應(yīng)；二次高通濾波結(jié)果與反演結(jié)果中深度10～45 km場(chǎng)源之間的對(duì)應(yīng)關(guān)系更加明顯.

對(duì)頻率域優(yōu)化濾波分離所得局部異常（見(jiàn)圖18b）進(jìn)行二次分離，結(jié)果如圖25所示.

圖24 小波域優(yōu)化位變?yōu)V波所得區(qū)域異常的二次分離結(jié)果

圖25 頻率域優(yōu)化濾波所得局部異常的二次分離結(jié)果

從圖25可以看出，頻率域優(yōu)化濾波所得局部異常的二次低通濾波結(jié)果與反演結(jié)果中深度10～50 km的場(chǎng)源分布規(guī)律一致，但幅值太小.因而，在二次高通濾波結(jié)果中，仍然殘留較多的低頻成分，與反演結(jié)果中淺部場(chǎng)源對(duì)應(yīng)關(guān)系不太明顯.

從實(shí)測(cè)資料分離實(shí)驗(yàn)結(jié)果可以看出，不論是采用低、中、高頻三個(gè)頻段進(jìn)行異常分離還是對(duì)分離后的結(jié)果進(jìn)行二次分離，小波域優(yōu)化位變?yōu)V波方法的分離結(jié)果均優(yōu)于頻率域優(yōu)化濾波和Butterworth方法.

5 結(jié)論

本文在小波域?yàn)V波理論和優(yōu)化濾波方法的基礎(chǔ)上，研究提出小波域優(yōu)化位變?yōu)V波方法，該方法具有空間變化濾波能力，在不同空間位置表現(xiàn)出不同的濾波器特性，適用于局部數(shù)據(jù)頻譜與全局?jǐn)?shù)據(jù)頻譜存在較大差異的重力異常分離問(wèn)題.理論模型重力異常分離實(shí)驗(yàn)表明，在局部數(shù)據(jù)頻譜與全局?jǐn)?shù)據(jù)頻譜差異較大時(shí)，該方法相對(duì)于頻率域優(yōu)化濾波法具有一定優(yōu)勢(shì).實(shí)測(cè)資料分離實(shí)驗(yàn)結(jié)果表明，該方法分離結(jié)果優(yōu)于頻率域優(yōu)化濾波和Butterworth濾波方法，與反演解釋成果吻合得更好.

小波域位變優(yōu)化濾波方法提供了設(shè)計(jì)小波域位變?yōu)V波器的有效途徑，該方法可應(yīng)用于其他需要位變?yōu)V波的位場(chǎng)數(shù)據(jù)處理領(lǐng)域.需要指出的是該位變?yōu)V波法會(huì)在分塊邊界處出現(xiàn)一定程度的階躍現(xiàn)象，有的時(shí)候需要對(duì)邊界附近的分離結(jié)果進(jìn)行平滑濾波以消除階躍現(xiàn)象，這是下一步實(shí)用化方面需要完善的地方.

An Y L, Guan Z N. 1985. The regularized stable factors of removing high frequency disturbances.ComputingTechniquesforGeophysicalandGeochemicalExploration（in Chinese）, 7（1）: 13-23.

Fedi M, Quarta T. 1998. Wavelet analysis for the regional-residual and local separation of potential field anomalies.GeophysicalProspecting, 46（5）: 507-525.

Guo L H, Meng X H, Shi L, et al. 2012. Preferential filtering method and its application to Bouguer gravity anomaly of Chinese continent.ChineseJ.Geophys. （in Chinese）, 55（12）: 4078-4088, doi: 10.6038/j.issn.0001-5733.2012.12.020.

Gupta V K, Ramani N. 1980. Some aspects of regional-residual separation of gravity anomalies in a Precambrian terrain.Geophysics, 45（9）: 1412-1426.

Hou Z C. 1981. Filtering of smooth compensation.GeophysicalProspectingforPetroleum（in Chinese）, （2）: 22-29. Kaiser G. 1996. Wavelet filtering in the scale domain.∥ Proceedings of the SPIE 2762, Wavelet Applications III. Orlando, FL: SPIE.

Keating P, Pinet N. 2011. Use of non-linear filtering for the regional-residual separation of potential field data.JournalofAppliedGeophysics, 73（4）: 315-322.

Li Z S. 2006. Study of the new high order Butterworth optimum transfer function.JournalofHuaqiaoUniversity（NaturalScience） （in Chinese）, 27（2）: 174-176.

Pawlowski R S, Hansen R O. 1990. Gravity anomaly separation by Wiener filtering.Geophysics, 55（5）: 539-548.

Pawlowski R S. 1994. Green′s equivalent-layer concept in gravity band-pass filter design.Geophysics, 59（1）: 69-76.

Ridsdill-Smith T A. 1998. Separating aeromagnetic anomalies using wavelet matched filters.∥ Proceedings of the 68th Annual International Meeting, SEG, Expanded Abstracts. 550-553.

Spector A, Grant F S. 1970. Statistical models for interpreting aeromagnetic data.Geophysics, 35（2）: 293-302.

Xiong G C. 1979. Spatially varying filtering and its application to interpretation of gravity and magnetic anomalies.GeophysicalandGeochemicalExploration（in Chinese）, 3（5）: 43-49.

Xiong G C, Zhang F R. 1981. Spatially varying filter.GeophysicalandGeochemicalExploration（in Chinese）, 5（4）: 205-213.

Yang W C, Shi Z Q, Hou Z Z, et al. 2001. Discrete wavelet transform for multiple decomposition of gravity anomalies.ChineseJ.Geophys. （in Chinese）, 44（4）: 534-542.

附中文參考文獻(xiàn)

安玉林, 管志寧. 1985. 濾除高頻干擾的正則化穩(wěn)定因子. 物化探計(jì)算技術(shù), 7（1）: 13-23.

郭良輝, 孟小紅, 石磊等. 2012. 優(yōu)化濾波方法及其在中國(guó)大陸布格重力異常數(shù)據(jù)處理中的應(yīng)用. 地球物理學(xué)報(bào), 55（12）: 4078-4088, doi: 10.6038/j.issn.0001-5733.2012.12.020.

侯重初. 1981. 補(bǔ)償圓滑濾波方法. 石油物探, （2）: 22-29.

李鐘慎. 2006. 新型高階Butterworth最佳傳遞函數(shù). 華僑大學(xué)學(xué)報(bào)（自然科學(xué)版）, 27（2）: 174-176.

熊光楚. 1979. 位變?yōu)V波及其在重磁異常解釋中的應(yīng)用. 物探與化探, 3（5）: 43-49.

熊光楚, 張福榮. 1981. 位變?yōu)V波器. 物探與化探, 5（4）: 205-213.

楊文采, 施志群, 侯遵澤等. 2001. 離散小波變換與重力異常多重分解. 地球物理學(xué)報(bào), 44（4）: 534-542.

（本文編輯 何燕）

Preferential spatially varying filtering method in the wavelet domain for gravity anomaly separation

LIU Cai-Yun1, 2, YAO Chang-Li2*, ZHENG Yuan-Man2

1SchoolofMathematicsandInformation,YangtzeUniversity,HubeiJingzhou434023,China2SchoolofGeophysicsandInformationTechnology,ChinaUniversityofGeosciences,Beijing100083,China

The classical frequency domain filtering method for gravity anomaly separation cannot change its frequency response at different spatial positions to adapt the frequency characteristic of local data, for the reason of lacking spatial information with Fourier transform. A preferential spatially varying filtering method in the wavelet domain （PSVF-WD） is proposed based on the scaling filtering theory and preferential filtering method, in order to overcome the limitation of the classical frequency domain filtering method mentioned above.This method uses a preferential spatially varying filter to separate gravity anomalies. Firstly, it segments gravity anomaly data into several blocks after analyzing the spatial distribution characteristics of frequency components with the wavelet analysis method. Secondly, it obtains the local translation function with the preferential filtering method and calculates the local equivalent coefficients with the method derived in this paper. Thirdly, it combines the local equivalent coefficients to global ones according to the position information of them and achieves the design of PSVF-WD. Lastly, it obtains separated gravity anomalies using the global equivalent coefficients of PSVF-WD and wavelet coefficients of gravity anomalies.We test the PSVF-WD with gravity-anomaly separation experiments of three synthetic data and one real data. In experiment 1, the PSVF-WD separates the composite signal of four different frequency sinusoidal signals using low-high filtering at one time. The results of this experiment show the spatially varying filtering ability of PSVF-WD. In experiment 2, the synthetic model consists of four prisms. Two of them belong to a relatively shallower layer and the other two of them belong to relatively deeper layer. However, one prism of the relatively shallower layer is deeper than another prism of the relatively deeper layer. The PSVF-WD method separates the anomalies of the relative shallower layer and deeper layer pretty well, while the preferential filtering method only separates the anomalies of the shallowest prism and remains anomalies of the rest of three prisms together. The average absolute error and standard deviation of separation results are 0.0500 nT and 0.0540, respectively for SVF-WD, while they are 0.0503 nT and 0.1079, respectively for preferential filtering. In experiment 3, the synthetic model consists of one large horizontal prism, one large vertical prism and five small prisms with different depths. The problem is more complicated because the spectrum aliasing is serious. The regional-residual separation results of PSVF-WD method are pretty well, but regional anomalies separated by the preferential filtering method still have too many residual anomalies and deform severely. The average absolute error and standard deviation of separation result are 0.0705 nT and 0.0531, respectively, while they are 0.0709 nT and 0.0867, respectively for preferential filtering.In the field data separation experiment with the PSVF-WD method, the separated regional anomalies contain large- and middle- scale anomalies, which correspond to the field sources deeper than 10 km described in inversion results, and the separated residual anomalies contain small-scale anomalies, which correspond to the field sources shallower than 10 km described in inversion results. In the experiments with preferential filtering and traditional Butterworth filtering, the separated regional anomalies are both too simple and do not correspond to inversion results very well, and the separated residual anomalies both contain too many regional anomalies.The proposed PSVF-WD method has the ability of spatially varying filtering, and is suitable for separating the gravity anomalies whose spectrum of local data is different from the average spectrum of global data obviously. The results of synthetic and field data separation experiments show that the proposed PSVF-WD method is superior to the classical frequency domain methods such as Butterworth filtering and preferential filtering method when the spectrum of local data is different from the average spectrum of global data obviously. In sum, this paper provides an approach to design a spatially varying filter in the wavelet domain, which can be applied to the other spatially varying filtering fields in potential data processing.

Anomaly separation; Wavelet domain; Preferential filtering; Spatially varying filtering; Gravity anomaly

國(guó)家高技術(shù)研究發(fā)展計(jì)劃（863計(jì)劃）課題（2014AA06A613），國(guó)家自然科學(xué)基金項(xiàng)目（61273179），湖北省教育廳科學(xué)技術(shù)研究項(xiàng)目（D20131206）聯(lián)合資助.

劉彩云，女，1975年生，博士、副教授，主要從事小波分析、重磁勘探和應(yīng)用數(shù)學(xué)等方面的教學(xué)和科研工作. E-mail: liucaiyun01@hotmail.com

*通訊作者 姚長(zhǎng)利，男，1965年生，教授，博士生導(dǎo)師，主要從事重磁勘探方法技術(shù)研究.E-mail: clyao@cugb.edu.cn

10.6038/cjg20151234.

10.6038/cjg20151234

P631

2014-12-01，2015-09-30收修定稿

劉彩云， 姚長(zhǎng)利， 鄭元滿. 2015. 重力異常分離的小波域優(yōu)化位變?yōu)V波方法.地球物理學(xué)報(bào)，58（12）:4740-4755,

Liu C Y, Yao C L, Zheng Y M. 2015. Preferential spatially varying filtering method in the wavelet domain for gravity anomaly separation.ChineseJ.Geophys. （in Chinese）,58（12）:4740-4755,doi:10.6038/cjg20151234.