優(yōu)化濾波方法及其在中國(guó)大陸布格重力異常數(shù)據(jù)處理中的應(yīng)用

郭良輝，孟小紅，石 磊，陳召曦

1 地下信息探測(cè)技術(shù)與儀器教育部重點(diǎn)實(shí)驗(yàn)室（中國(guó)地質(zhì)大學(xué)，北京），北京 100083

2 中國(guó)地質(zhì)大學(xué)（北京）地球物理與信息技術(shù)學(xué)院，北京 100083

3 中國(guó)地震局地球物理研究所，北京 100081

1 引 言

根據(jù)重力觀測(cè)值整理得到的重力異常，包含了從地表到深部所有密度不均勻體引起的重力效應(yīng)，信息非常豐富.實(shí)測(cè)重力異常是所有這些重力效應(yīng)的疊加.要根據(jù)實(shí)測(cè)重力異常反演解釋某個(gè)地質(zhì)體（比如巖礦體、地質(zhì)構(gòu)造或物性分界面），必須首先從疊加異常中分離出單純由這個(gè)地質(zhì)體引起的異常，然后用這個(gè)異常進(jìn)行反演解釋.然而，從疊加異常中分離出某個(gè)地質(zhì)體引起的異常，或者把疊加異常分解為幾個(gè)地質(zhì)體引起的單一異常，是比較困難的.重力異常分離是目前重力資料處理解釋中沒(méi)有很好解決的難題之一.

當(dāng)前，重力異常分離方法眾多，如最小二乘平滑法、帶通濾波法、向上延拓法、匹配濾波法［1］、補(bǔ)償圓滑濾波法［2］、正則化濾波法［3］、維納濾波法［4-5］、小波變換法［6-8］、熵濾波法［9-10］、非線性濾波法［11］等.這些方法均有各自的特點(diǎn)和優(yōu)勢(shì)，但同時(shí)又有各自的局限性.

傳統(tǒng)向上延拓法存在的一個(gè)主要問(wèn)題是在進(jìn)行向上延拓時(shí)，淺源短波長(zhǎng)信息和深源長(zhǎng)波長(zhǎng)信息同時(shí)延拓，在淺源信息被壓制的同時(shí)，深源信息也受到壓制，將實(shí)測(cè)異常減去上延后的區(qū)域異常所得到的剩余異常中，仍包含有區(qū)域成分，異常分離不徹底［12-13］.針對(duì)該問(wèn)題，Pawlowski［12］根據(jù)維納濾波和格林等效層原理提出了優(yōu)選延拓法，應(yīng)用該方法進(jìn)行向上延拓時(shí)，可以在壓制淺源短波長(zhǎng)信息的同時(shí)，保持深源長(zhǎng)波長(zhǎng)信息不衰減.此后，許德樹(shù)等［14］、孟小紅等［15］在優(yōu)選延拓法基礎(chǔ)上給出了基于優(yōu)選向上延拓的延拓差值法，實(shí)現(xiàn)對(duì)某一指定頻段或尺度的局部異常有選擇地提取或分離，即類似于帶通濾波.傳統(tǒng)向上延拓法存在的另一個(gè)重要問(wèn)題是延拓高度必須已知，然而實(shí)際應(yīng)用中延拓高度往往是未知的［13］，優(yōu)選延拓法同樣存在此問(wèn)題.針對(duì)這個(gè)問(wèn)題，Zeng等［13］提出了向上延拓法的最佳延拓高度估計(jì)方法，該方法分別以一系列遞增的延拓高度對(duì)實(shí)測(cè)異常做向上延拓，得到相應(yīng)的延拓場(chǎng)，計(jì)算各前后不同延拓高度的延拓場(chǎng)的互相關(guān)，這些互相關(guān)值的最大偏差點(diǎn)所對(duì)應(yīng)的高度即為最佳延拓高度.郭良輝等［16］和 Meng等［15］給出了優(yōu)選向上延拓的最佳延拓高度估計(jì)辦法.

由于延拓高度往往是未知的，這使得傳統(tǒng)向上延拓法和優(yōu)選延拓法的實(shí)際應(yīng)用受到一定限制.本文在優(yōu)選延拓法的理論基礎(chǔ)上，進(jìn)一步研究提出基于維納濾波和格林等效層原理的優(yōu)化濾波法，該方法分離重力異常不需要已知延拓高度，可對(duì)重力異常進(jìn)行指定頻段的低通、帶通和高通濾波.

本文首先給出優(yōu)化濾波法的方法原理和計(jì)算步驟，然后利用該方法對(duì)理論模型的含噪重力異常進(jìn)行去噪和分離異常試驗(yàn)，并與傳統(tǒng)向上延拓法和帶通濾波法的處理結(jié)果進(jìn)行對(duì)比，最后利用優(yōu)化濾波法對(duì)中國(guó)大陸重力異常數(shù)據(jù)進(jìn)行去噪和異常分離，并對(duì)分離出的區(qū)域重力異常做密度界面約束反演，得到中國(guó)大陸莫霍面深度分布特征.

2 優(yōu)化濾波方法

2.1 優(yōu)化濾波算子

假設(shè)實(shí)測(cè)重力異常gobs（x，y）是深源長(zhǎng)波長(zhǎng)異常gd（x，y）與淺源短波長(zhǎng)異常gs（x，y）的疊加：

那么實(shí)測(cè)重力異常的功率譜密度函數(shù)P可看成是深源異常功率譜密度函數(shù)Pd與淺源異常功率譜密度函數(shù)Ps之和［5］，

Pawlowski［5］指出實(shí)測(cè)重力異常的傅里葉功率譜可通過(guò)引入格林等效層概念來(lái)模擬，即利用處于不同深度的m＋n層格林等效源薄層來(lái)建立實(shí)測(cè)重力異常的徑向平均對(duì)數(shù)功率譜模型（假定各等效層互不相關(guān)），

其中，深源場(chǎng)成分用m層格林等效層來(lái)模擬，而淺源場(chǎng)分量用n層格林等效層來(lái)模擬，誤差E來(lái)源于對(duì)實(shí)測(cè)數(shù)據(jù)的功率譜的不完全擬合.顯然，對(duì)于不同地區(qū)，充分?jǐn)M合時(shí)等效層的數(shù)目將有所不同.

按照 Naidu［17］和 Dampney［18］的定義，任意深度h處格林等效層的功率譜密度函數(shù)可寫成

其中，s是與期望的等效層強(qiáng)度成正比的常數(shù)，k是角頻率.

設(shè)目標(biāo)層為深部i→j（n＋m≥j≥i≥1）等效層，則優(yōu)化濾波法的期望是保持目標(biāo)層場(chǎng)源信息不變，而對(duì)其它層場(chǎng)源信息做壓制，即，

設(shè)維納濾波器的輸入是實(shí)測(cè)異常gobs（x，y），期望輸出是gpf（x，y），按照維納濾波器設(shè)計(jì)的要求，濾波器的頻率響應(yīng)［19-20］為

其中，G代表頻率域的頻譜，〈〉代表數(shù)學(xué)期望，＊代表復(fù)共軛.式中的分子是濾波器期望輸出和輸入信號(hào)（實(shí)測(cè)異常）的互功率頻譜密度函數(shù)，分母是濾波器輸入的功率頻譜密度函數(shù).

結(jié)合公式（3），式（6）的維納濾波轉(zhuǎn)換函數(shù)可改寫為

其中，P（i→j）為目標(biāo)層場(chǎng)源信息的功率譜密度函數(shù)，P′（i→j）是目標(biāo)層場(chǎng)源信息與其它層場(chǎng)源信息的互功率譜密度函數(shù).

假設(shè)目標(biāo)層與其它深度層場(chǎng)源信息互不相關(guān)，則P′（i→j）＝0，那么公式（7）的維納濾波轉(zhuǎn)換函數(shù)可進(jìn)一步改寫為

稱公式（8）的維納濾波轉(zhuǎn)換函數(shù)為優(yōu)化濾波算子.當(dāng)i＝1時(shí)，i→j即為1→j，此時(shí)的優(yōu)化濾波是分離1→j層場(chǎng)源信息的低通濾波；當(dāng)j＝m＋n時(shí)，i→j即為i→m＋n，此時(shí)的優(yōu)化濾波是分離i→m＋n層場(chǎng)源信息的高通濾波；當(dāng)n＋m＞j≥i＞1時(shí)，此時(shí)的優(yōu)化濾波是分離i→j層場(chǎng)源信息的帶通濾波.因此，利用公式（8）可對(duì)實(shí)測(cè)重力異常進(jìn)行不同形式的優(yōu)化濾波，分離出目標(biāo)層場(chǎng)源在原始觀測(cè)面上的異常信息.

2.2 優(yōu)化濾波計(jì)算步驟

優(yōu)化濾波法的計(jì)算步驟如下：

（1）對(duì)實(shí)測(cè)異常做傅里葉變換，計(jì)算徑向平均對(duì)數(shù)功率譜；

（2）用分段線性擬合法擬合實(shí)測(cè)徑向平均對(duì)數(shù)功率譜；

首先分析實(shí)測(cè)徑向平均對(duì)數(shù)功率譜形狀特征，確定格林等效源層數(shù)和各等效層分布的頻率范圍.一般地，徑向頻率越低，頻段分割越細(xì)窄，段數(shù)越多，而徑向頻率越高，頻段分割越粗糙，段數(shù)越少.每頻段對(duì)應(yīng)著一個(gè)格林等效層，低頻段主要對(duì)應(yīng)區(qū)域深源等效層，中頻段主要對(duì)應(yīng)淺源等效層，而高頻段主要對(duì)應(yīng)噪聲［1，4-5］.然后利用直線擬合實(shí)測(cè)徑向平均對(duì)數(shù)功率譜的每個(gè)分段功率譜.最后利用（9）式估計(jì)各等效層的近似深度［1，21-22］：

其中，r1和r2分別是擬合直線的徑向頻率起點(diǎn)和終點(diǎn)，P（r1）和P（r2）分別是r1和r2對(duì)應(yīng)的實(shí)測(cè)徑向平均對(duì)數(shù)功率譜.

（3）利用公式（4）表示的功率譜密度函數(shù)模型和所估計(jì)出的各等效層深度值來(lái)擬合實(shí)測(cè)徑向平均對(duì)數(shù)功率譜；

設(shè)a＝e－2kh，公式（4）對(duì)于第i個(gè)徑向頻率點(diǎn)可寫成

其中，Pi是第i個(gè)徑向頻率點(diǎn)的徑向平均對(duì)數(shù)功率譜，ei是功率譜擬合誤差.

對(duì)于l個(gè)功率譜樣點(diǎn)，有l(wèi)線性方程構(gòu)成的方程組，寫成矩陣如下：

通過(guò)共軛梯度法等迭代優(yōu)化算法對(duì)公式（12）求解可得到向量S，把它代入公式（4）可計(jì)算任意第i個(gè)徑向頻率點(diǎn)或第i個(gè)格林等效層的理論徑向平均對(duì)數(shù)功率譜.

（4）利用公式（8）構(gòu)建優(yōu)化濾波算子；

（5）對(duì)實(shí)測(cè)異常進(jìn)行優(yōu)化濾波處理；

（6）進(jìn)行反傅里葉變換，得到空間域的目標(biāo)異常.

3 理論模型數(shù)據(jù)試驗(yàn)

圖1 （a）理論模型各長(zhǎng)方體在水平面位置分布，（b）理論重力異常等值線圖，（c）含噪重力異常等值線圖.重力異常單位：mGal.Fig.1 The outline of the model on the horizontal plane（a），the theoretical gravity anomaly of the model（b），and the noised gravity anomaly（c）.The unit of gravity anomaly is mGal.

理論模型由三個(gè)不同深度層、不同大小和剩余密度值的11個(gè)直立長(zhǎng)方體組合而成，圖1a顯示了各直立長(zhǎng)方體在水平面上的投影分布.假設(shè)深層長(zhǎng)方體A1、A2和A3產(chǎn)生區(qū)域重力異常，中層長(zhǎng)方體B1、B2、B3、B4、B5和淺層長(zhǎng)方體C1、C2和C3共同產(chǎn)生局部重力異常（或剩余重力異常），理論模型的理論重力異常等值線圖見(jiàn)圖1b，對(duì)理論重力異常加入基準(zhǔn)值4%的高斯隨機(jī)噪聲見(jiàn)圖1c.

對(duì)含噪重力異常計(jì)算徑向?qū)?shù)功率譜（圖2），通過(guò)功率譜形狀分析，用0～0.2745cycles·km-1（頻段1）、0.2745～1.2157cycles·km-1（頻段2）和1.2157～5cycles·km-1（頻段3）三個(gè)等效層分段擬合（圖2）.這里，頻段1的功率譜主要對(duì)應(yīng)A層長(zhǎng)方體的低頻異常信息，頻段2的功率譜主要對(duì)應(yīng)B和C層長(zhǎng)方體的中高頻信息，而頻段3的功率譜主要對(duì)應(yīng)高斯噪聲的高頻干擾信息.

下面對(duì)含噪重力異常進(jìn)行優(yōu)化濾波去噪和分離異常試驗(yàn).首先進(jìn)行頻段1和2的優(yōu)化濾波（低通濾波），目的是壓制高斯噪聲，優(yōu)化濾波算子的頻率響應(yīng)見(jiàn)圖3a，濾波結(jié)果見(jiàn)圖3d.從圖可見(jiàn)，優(yōu)化濾波法有效壓制了高斯噪聲干擾，除C層有效信號(hào)保留的不是很好之外（因?yàn)槠漕l率與噪聲的較接近），A和B層的有效信號(hào)都較好地保留了.然后進(jìn)行頻段1的優(yōu)化濾波（低通濾波），目的是分離出A層各長(zhǎng)方體引起的區(qū)域異常，優(yōu)化濾波算子的頻率響應(yīng)見(jiàn)圖3b，濾波結(jié)果見(jiàn)圖3e.從圖可見(jiàn)，優(yōu)化濾波法有效壓制了高斯噪聲干擾及B和C層的局部異常，A層的有效信號(hào)被較好地保留.最后進(jìn)行頻段2的優(yōu)化濾波（帶通濾波），目的是分離出B和C層各長(zhǎng)方體共同引起的局部異常，優(yōu)化濾波算子的頻率響應(yīng)見(jiàn)圖3c，濾波結(jié)果見(jiàn)圖3f.從圖可見(jiàn)，優(yōu)化濾波法有效壓制了高斯噪聲干擾及A層的區(qū)域異常，B和C層的有效信號(hào)得到較好的保留.

圖2 含噪重力異常徑向?qū)?shù)功率譜及其分段擬合黑色：含噪重力異常功率譜；藍(lán)色：分頻段擬合；紅色：模型擬合的功率譜.Fig.2 The radially averaged logarithm power spectrum of the noised gravity anomaly（black line）and its fitting by using piece－wise linearizationThe blue lines show the fitted straight lines in each radial frequency section，and the red line indicates the fitted radial logarithm power spectrum by using power spectrum density function models.

圖3 壓制高斯噪聲的優(yōu)化濾波算子頻率響應(yīng)（a）和濾波結(jié)果（d），分離A層長(zhǎng)方體異常的優(yōu)化濾波算子頻率響應(yīng)（b）和濾波結(jié)果（e），分離B和C層長(zhǎng)方體異常的優(yōu)化濾波算子頻率響應(yīng)（c）和濾波結(jié)果（f）.重力異常單位：mGal.Fig.3 The amplitude response of the preferential filtering operator for suppressing Gaussian noise （a），separating regional anomalies of layer A （b），and separating local anomalies of layers B and C（c）.（d），（e）and（f）show the corresponding anomalies of（a），（b）and（c）.The unit of gravity anomaly is mGal.

圖4 帶通濾波法的壓制高斯噪聲（a）、A層的區(qū)域異常（b）、B和C層的局部異常（c），傳統(tǒng)向上延拓法的壓制高斯噪聲（d）、區(qū)域異常（e）和局部異常（f）.重力異常單位：mGal.Fig.4 The denoised anomalies（a），the separated regional anomalies（b）and the separated local anomalies（c）by using the low－pass filtering，and those（d，e，f）of the upward continuation.The unit of gravity anomaly is mGal.

下面，分別利用帶通濾波法和傳統(tǒng)向上延拓法對(duì)理論含噪重力異常進(jìn)行同樣的去噪和異常分離試驗(yàn).壓制高斯噪聲時(shí)，帶通濾波法的濾波波長(zhǎng)寬度為800m，向上延拓高度為300m，結(jié)果分別見(jiàn)圖4（a，d）.從圖可見(jiàn)，帶通濾波法有效壓制了高斯噪聲干擾，保留有效信號(hào)，而傳統(tǒng)向上延拓法在壓制高斯噪聲干擾的同時(shí)，也部分壓制了有效信號(hào).分離A層的區(qū)域異常時(shí)，帶通濾波法的濾波波長(zhǎng)寬度為4000m，向上延拓高度為1000m，結(jié)果分別見(jiàn)圖4b和（e）.從圖可見(jiàn)，帶通濾波法和傳統(tǒng)向上延拓法在有效壓制高頻噪聲干擾及B和C層局部異常的同時(shí)，也部分壓制了A層的區(qū)域異常，尤其是傳統(tǒng)向上延拓法的區(qū)域異常被壓制較多.分別用兩種方法壓制高斯噪聲后的異常減去區(qū)域異常，即可得到分離B和C層的局部異常，結(jié)果分別見(jiàn)圖4（c，f）.從圖可見(jiàn)，帶通濾波法、傳統(tǒng)向上延拓法都能分離出B和C層的局部異常，但也遺留了A層的部分區(qū)域異常，尤其是傳統(tǒng)向上延拓法的區(qū)域異常遺留較多，也就是說(shuō)異常分離不徹底.

圖5顯示了這三種不同方法分離出的異常與理論重力異常（圖1b，沒(méi)有噪聲干擾）沿X＝5680m剖面的對(duì)比圖.從圖4和圖5可見(jiàn)，在壓制高斯噪聲方面，三種方法分離效果基本一致，但傳統(tǒng)向上延拓法部分壓制了有效信號(hào)；在分離A層區(qū)域異常方面，優(yōu)化濾波法較好保留了區(qū)域異常，帶通濾波法部分壓制了區(qū)域異常，而傳統(tǒng)向上延拓法較多壓制了區(qū)域異常；在分離B和C層局部異常方面，優(yōu)化濾波法結(jié)果最接近理論值，而帶通濾波法部分遺留了區(qū)域異常，傳統(tǒng)向上延拓法較多遺留了區(qū)域異常.

圖5 不同方法在X＝5680m剖面的異常分離結(jié)果與理論重力異常對(duì)比（a）壓制高斯噪聲；（b）A層的區(qū)域異常；（c）B和C層的局部異常.藍(lán)色：理論值；紅色：優(yōu)化濾波法；綠色：帶通濾波法；黃色：傳統(tǒng)向上延拓法.Fig.5 Comparisons between the theoretical gravity anomaly（blue lines）and the separated anomalies by using the three methods along profile X＝5680m（a）Denoised anomalies；（b）Separated regional anomalies；（c）Separated local anomalies.Red lines correspond to the preferential filtering，green lines to the low－pass filtering，and yellow lines to the upward continuation.

4 中國(guó)大陸重力異常數(shù)據(jù)試驗(yàn)

依托“深部探測(cè)技術(shù)與試驗(yàn)研究”項(xiàng)目—第五子課題（Sinoprobe－01－05），開(kāi)展了中國(guó)大陸重力異常數(shù)據(jù)的應(yīng)用實(shí)驗(yàn)研究.利用最新的地球重力場(chǎng)模型EGM2008（Earth Gravitational Model 2008）［23］計(jì)算了中國(guó)大陸自由空氣重力異常.EGM2008是由美國(guó)國(guó)家地理空間情報(bào)局（NGA）推出的官方版高階地球重力場(chǎng)模型，該模型的階次完全至2159（另外球諧系數(shù)的階擴(kuò)展至2190，次為2159），空間分辨率約為5′.由于采用了GRACE衛(wèi)星跟蹤數(shù)據(jù)、衛(wèi)星測(cè)高數(shù)據(jù)及地面重力數(shù)據(jù)等，使得該模型無(wú)論在精度還是在分辨率方面均取得了巨大的進(jìn)步.該模型在海域的精度可達(dá)到1∶100萬(wàn)比例尺重力勘探的要求［23］，在中國(guó)大陸大部分地區(qū)精度可達(dá)10mGal［24］.因此，EGM2008可用于中等比例尺的海洋地質(zhì)與資源調(diào)查和小比例尺的大陸重力編圖與構(gòu)造研究.

利用EGM2008，計(jì)算了中國(guó)大陸自由空氣重力異常（網(wǎng)度為10km），然后按照區(qū)域重力調(diào)查規(guī)范做地形校正和布格校正（中間層校正），最終得到中國(guó)大陸布格重力異常數(shù)據(jù)（網(wǎng)度為10km），見(jiàn)圖6所示.校正所用的地形高程數(shù)據(jù)（圖10）來(lái)自于美國(guó)斯克里普斯海洋研究所網(wǎng)站（http：／／topex.ucsd.edu／cgi－bin／get＿data.cgi），空間分辨率為1′.從圖6可見(jiàn)，校正后得到的布格重力異常在中國(guó)大陸范圍內(nèi)普遍存在高頻噪聲干擾，經(jīng)分析認(rèn)為，這些噪聲干擾一部分是由EGM2008計(jì)算的自由空氣重力異常數(shù)據(jù)自身攜帶而遺留下來(lái)的，另一部分來(lái)源于地形校正和布格校正，而后者主要是由于自由空氣重力異常數(shù)據(jù)和地形數(shù)據(jù)的分辨率與精度不一致造成的.顯然，高頻噪聲干擾對(duì)后續(xù)數(shù)據(jù)處理解釋及應(yīng)用等將產(chǎn)生嚴(yán)重影響，應(yīng)首先對(duì)重力異常數(shù)據(jù)進(jìn)行去噪處理.

對(duì)圖6所示的中國(guó)大陸布格重力異常數(shù)據(jù)嘗試采用了傳統(tǒng)向上延拓法、帶通濾波法、巴特沃斯濾波法等多種常規(guī)方法進(jìn)行去噪處理，但結(jié)果都不理想，要么噪聲沒(méi)有壓制干凈，要么有效信號(hào)部分損失，利用本文的優(yōu)化濾波法對(duì)中國(guó)大陸布格重力異常數(shù)據(jù)進(jìn)行去噪試驗(yàn)，取得較滿意結(jié)果.

圖6 EGM2008計(jì)算并經(jīng)過(guò)校正得到的中國(guó)大陸布格重力異常圖（網(wǎng)度為10km）Fig.6 The Bouguer gravity anomaly of Chinese continent calculated from EGM 2008 with terrain correction and Bouguer correction.The grid size is 10km.

首先計(jì)算中國(guó)大陸布格重力異常的徑向平均對(duì)數(shù)功率譜（圖7），根據(jù)功率譜形狀特征，用0～0.0014cycles·km-1（頻段 1）、0.0014～0.0039 cycles·km-1（頻段2）、0.0039～0.0099cycles·km-1（頻段3）、0.0099～0.0174cycles·km-1（頻段4）、0.0174～0.0249cycles·km-1（頻段5）、0.0249～0.05cycles·km-1（頻段6）6個(gè)等效層分段擬合.分析認(rèn)為，頻段5和頻段6主要對(duì)應(yīng)高頻噪聲干擾，頻段1～4主要對(duì)應(yīng)有效信號(hào).然后進(jìn)行頻段1～4的優(yōu)化濾波（低通濾波），目的是壓制高頻噪聲干擾，得到的去噪布格重力異常見(jiàn)圖8所示.從圖可見(jiàn)，優(yōu)化濾波法有效壓制了噪聲干擾，又較好地保留了深部有效信息.通過(guò)對(duì)比，去噪后的布格重力異常與前人成果［25］基本一致.之后進(jìn)行頻段1～2的優(yōu)化濾波（低通濾波），目的是分離出中國(guó)大陸區(qū)域重力異常，結(jié)果見(jiàn)圖9所示.

從圖8和9可見(jiàn)，中國(guó)境內(nèi)布格重力異常值自東向西逐漸減小.布格重力異常值在海岸線為0mGal左右；由海岸線向西，異常值緩慢遞減，并進(jìn)入負(fù)值區(qū)；沿大興安嶺—太行山—武陵山一帶顯示為過(guò)渡帶，這是我國(guó)一條重要的重力梯級(jí)帶，走向?yàn)镹E向；再往西，至青藏高原周邊地區(qū)，異常值迅速減小，這是另一條重要的重力梯級(jí)帶，圍繞青藏高原呈弧形展布；青藏高原大部分地區(qū)的異常值小于-400mGal.可見(jiàn)，由優(yōu)化濾波法去噪和分離出的中國(guó)大陸布格重力異常，能較好反映出全國(guó)重力場(chǎng)變化特征，提供關(guān)于地殼和上地幔中物質(zhì)質(zhì)量分布的信息.其中，區(qū)域布格重力異常（圖9）主要由莫霍面引起的，反映了莫霍面起伏的基本特征.

圖7 中國(guó)大陸布格重力異常徑向?qū)?shù)功率譜及其分段擬合黑色：中國(guó)大陸布格重力異常功率譜；藍(lán)色：分頻段擬合；紅色：擬合的功率譜.Fig.7 The radially averaged logarithm power spectrum of the Bouguer gravity anomaly of Chinese continent（black line）and its fitting by using piece－wise linearization The blue lines show the fitted straight lines in each radial frequency section，and the red line indicates the fitted radial logarithm power spectrum by using the power spectrum density function models.

圖8 優(yōu)化濾波去噪后的中國(guó)大陸布格重力異常（網(wǎng)度為10km）Fig.8 The denoised Bouguer gravity anomaly of Chinese continent by the preferential filtering method.The grid size is 10km

圖9 優(yōu)化濾波分離出的中國(guó)大陸區(qū)域重力異常（網(wǎng)度為10km）Fig.9 The separated regional gravity anomaly of Chinese continent by the preferential filtering method.The grid size is 10km

為認(rèn)識(shí)中國(guó)大陸莫霍面分布特征，利用優(yōu)化濾波法分離出的區(qū)域重力異常（圖9）反演了中國(guó)大陸莫霍面深度，其中反演方法采用帶先驗(yàn)信息約束的非線性回歸法［26］，以減少反演的多解性.自20世紀(jì)60年代以來(lái)我國(guó)開(kāi)展了大量的深地震測(cè)深、深地震反射和寬頻地震觀測(cè)等地震探測(cè)，其具有較大的探測(cè)深度和較高的垂向分辨率，為研究中國(guó)大陸深部結(jié)構(gòu)和構(gòu)造格架提供了及其重要的信息.對(duì)這些人工地震探測(cè)推斷的莫霍面深度圖進(jìn)行了收集、整理和矢量化，最終得到中國(guó)大陸莫霍面深度控制點(diǎn)數(shù)據(jù)庫(kù)［27］，數(shù)據(jù)點(diǎn)分布見(jiàn)圖10所示，它們將作為本文莫霍面深度約束反演的已知控制點(diǎn)信息.

以圖10所示的已知控制點(diǎn)信息作為約束，采用非線性回歸法對(duì)優(yōu)化濾波法分離出的區(qū)域重力異常（圖9）進(jìn)行莫霍面深度約束反演，得到中國(guó)大陸莫霍面深度分布，見(jiàn)圖11所示.反演結(jié)果的均方差為4.58km.通過(guò)對(duì)比分析，本文反演得到的莫霍面深度分布與前人成果［28］基本一致，但細(xì)節(jié)稍多些.從圖11可見(jiàn)，以東經(jīng)105°～110°的南北地震帶為界，東西兩側(cè)的中國(guó)大陸莫霍面深度明顯不同，以東地區(qū)相對(duì)較淺，介于30～45km，莫霍面等值線主要呈北東向展布，以西地區(qū)相對(duì)較深，在45～74km之間，莫霍面等值線主要呈東西向展布，體現(xiàn)出東西不同的深部結(jié)構(gòu)和構(gòu)造動(dòng)力學(xué)背景.由于篇幅限制，重力約束反演得到的中國(guó)大陸莫霍面深度特征（圖11）將另文詳細(xì)描述.

圖10 中國(guó)大陸地形圖及已知莫霍面深度控制點(diǎn)分布（紅點(diǎn)）Fig.10 The topography of Chinese continent and positions of the Moho depth control points deduced by the deep seismic detection（red points）

圖11 重力反演得到的中國(guó)大陸莫霍面深度分布（網(wǎng)度為10km）Fig.11 The Moho depth distribution of Chinese continent from gravity inversion.The grid size is 10km

5 結(jié) 論

本文在優(yōu)選延拓法的理論基礎(chǔ)上，研究提出了基于格林等效層概念和維納濾波器的優(yōu)化濾波法，用于對(duì)重力異常去噪和分離，給出了優(yōu)化濾波算子的構(gòu)建和優(yōu)化濾波的計(jì)算步驟.理論模型數(shù)據(jù)試驗(yàn)表明本文方法有效，異常分離效果優(yōu)于傳統(tǒng)的帶通濾波法和向上延拓法.通過(guò)中國(guó)大陸重力異常數(shù)據(jù)的優(yōu)化濾波去噪和異常分離，得到有效的布格重力異常數(shù)據(jù)和區(qū)域重力異常數(shù)據(jù)，對(duì)后者進(jìn)行密度界面約束反演得到中國(guó)大陸莫霍面深度分布.

與優(yōu)選延拓法相比，優(yōu)化濾波法與延拓高度無(wú)關(guān)，其異常分離不受延拓高度困擾，具有一定的優(yōu)勢(shì).當(dāng)然，優(yōu)化濾波法也有其缺點(diǎn)，如利用格林等效層近似地下場(chǎng)源，且假設(shè)目標(biāo)層與其他深度層場(chǎng)源信息互不相關(guān)，而實(shí)際地質(zhì)中，不同深度層場(chǎng)源可能存在部分相關(guān)，比如淺層場(chǎng)源有可能向下延伸到深部；另外，不同深度層場(chǎng)源的頻譜都是寬頻帶的，且互相重疊，這樣，很難從疊加的頻譜中徹底分離出目標(biāo)層場(chǎng)源的單一頻譜.

（References）

［1］Spector A，Grant F S.Statistical models for interpreting aeromagnetic data.Geophysics，1970，35（2）：293－302.

［2］侯重初.補(bǔ)償圓滑濾波方法.石油物探，1981，（2）：22－29.Hou Z C.Filtering of smooth compensation.Geophysical Prospecting for Petroleum （in Chinese），1981，（2）：22－29.

［3］安玉林，管志寧.濾除高頻干擾的正則化穩(wěn)定因子.物探化探計(jì)算技術(shù)，1985，7（1）：13－23.An Y L，Guan Z N.The regularized stable factors of removing high frequency disturbances.Computing Techniques for Geophysical and Geochemical Exploration （in Chinese），1985，7（1）：13－23.

［4］Pwalowski R S，Hansen R O.Gravity anomaly separation by Wiener filtering.Geophysics，1990，55（5）：539－548.

［5］Pwalowski R S.Green′s equivalent－layer concept in gravity band－pass filter design.Geophysics，1994，55（5）：69－76.

［6］侯遵澤，楊文采.中國(guó)重力異常的小波變換與多尺度分析.地球物理學(xué)報(bào)，1997，40（1）：85－95.Hou Z Z，Yang W C.Wavelet transform and multi－scale analysis on gravity anomalies of China.Chinese J.Geophys.（in Chinese），1997，40（1）：85－95.

［7］Fedi M，Quarta T.Wavelet analysis for the regional－residual and local separation of potential field anomalies.Geophysical Prospecting，1998，46（5）：507－525.

［8］楊文采，施志群，侯遵澤等.離散小波變換與重力異常多重分解.地球物理學(xué)報(bào)，2001，44（4）：534－541.Yang W C，Shi Z Q，Hou Z Z，et al.Discrete wavelet transform for multiple decomposition of gravity anomalies.Chinese J.Geophys.（in Chinese），2001，44（4）：534－541.

［9］Nikitin A A，Vasov O K，Belov A P，et al.Vozmozhnosti kompleksnoy geofizicheskoy interpretatsii na baze entropiynogo fil′tra.Izvestiya Akademii Nauk Turkmenskoy SSR.Seriya Fiziko－Tekhnicheskikh，Khimicheskikhi Geologicheskikh Nauk（in Russia），1984，2：79－82.

［10］郭良輝，孟小紅，石磊.磁異常ΔT三維相關(guān)成像.地球物理學(xué)報(bào)，2010，53（2）：435－441.Guo L H，Meng X H，Shi L.3Dcorrelation imaging for magnetic anomaly ΔT data.Chinese J.Geophys. （in Chinese），2010，53（2）：435－441.

［11］Keating P，Pinet N.Use of non－linear filtering for the regional－residual separation of potential field data.Journal of Applied Geophysics，2011，73（4）：315－322.

［12］Pwalowski R S.Preferential continuation for potential－field anomaly enhancement.Geophysics，1995，60（2）：390－398.

［13］Zeng H，Xu D，Tan H.A model study for estimating optimum upward－continuation height for gravity separation with application to a Bouguer gravity anomaly over a mineral deposit，Jilin province，northeast China.Geophysics，2008，72（4）：145－150.

［14］許德樹(shù)，曾華霖.優(yōu)選延拓技術(shù)及其在中國(guó)布格重力異常圖處理上的應(yīng)用.現(xiàn)代地質(zhì)，2000，14（2）：215－222.Xu D S，Zeng H L.Preferential continuation and its application to Bouguer gravity anomaly in China.Geoscience（in Chinese），2000，14（2）：215－222.

［15］Meng X H，Guo L H，Chen Z X，et al.A method for gravity anomaly separation based on preferential continuation and its application.Applied Geophysics，2009，6（3）：217－225.

［16］郭良輝，孟小紅，陳召曦.優(yōu)選向上延拓及其延拓高度估計(jì).CPS／SEG Beijing 2009International Geophysical Conference＆Exposition，2009，ID1164.Guo L H，Meng X H，Chen Z X.Preferential upward continuation and the estimation of its continuation height.CPS／SEG Beijing 2009International Geophysical Conference＆Exposition，2009，ID1164.

［17］Naidu P.Spectrum of the potential field due to randomly distributed sources.Geophysics，1968，33（2）：337－345.

［18］Dampney C N G. The equivalent source technique.Geophysics，1969，34（1）：39－53.

［19］Wiener N.Extrapolation，interpolation，and smoothing of stationary time series.New York：John Wiley ＆Sons Inc，1949.

［20］Clarke G K C.Optimum second－derivative and downward continuation filters.Geophysics，1969，34（3）：424－437.

［21］Boler F M.Aeromagnetic measurements，magnetic source depths，and the Curie point isotherm in the Vale－Owyhee，Oregon［Master thesis］.Corvallis：Oregon State University，1978.

［22］Connard G，Couch R，Gemperle M.Analysis of aeromagnetic measurements from the Cascade Range in central Oregon.Geophysics，1983，48（3）：376－390.

［23］張明華，張家強(qiáng).現(xiàn)代衛(wèi)星測(cè)高重力異常分辨能力分析及在海洋資源調(diào)查中應(yīng)用.物探與化探，2005，29（4）：295－298，303.Zhang M H，Zhang J Q.Resolution of modern satellite altimetric gravity anomaly and its application to marine geological survey.Geophysical ＆ Geochemical Exploration（in Chinese），2005，29（4）：295－298，303.

［24］楊金玉，張訓(xùn)華，張菲菲等.EGM2008地球重力模型數(shù)據(jù)在中國(guó)大陸地區(qū)的精度分析.地球物理學(xué)進(jìn)展，2012，27（4）：1298－1306.Yang J Y，Zhang X H，Zhang F F，et al.On the accuracy of EGM2008earth gravitational model in Chinese Mainland.Progress in Geophysics （in Chinese），2012，27（4）：1298－1306.

［25］袁學(xué)誠(chéng).中國(guó)地球物理圖集.北京：地質(zhì)出版社，1996.Yuan X C.Chinese Geophysical Atlas（in Chinese）.Beijing：Geological Publishing House，1996.

［26］楊永.利用重、磁異常研究東海南部中生界分布［碩士論文］.西安：長(zhǎng)安大學(xué)，2010.Yang Y.Research on Mesozoic strata in the southern part of east China sea by gravity and magnetic anomalies［Master thesis］（in Chinese）.Xi′an：Chang′an University，2010.

［27］熊小松.中國(guó)大陸莫霍面深度與變化特征及其地球動(dòng)力學(xué)意義［博士論文］.北京：中國(guó)地質(zhì)科學(xué)院，2010.Xiong X S.Moho depth and variation of the continent in China and its geodynamic implications［Ph.D.thesis］（in Chinese）.Beijing：China Academy of Geological Science，2010.

［28］Li S，Mooney W D，F(xiàn)an J.Crustal structure of mainland China from deep seismic sounding data.Tectonophysics，420（1－2）：239－252.