楊文采，孫艷云，侯遵澤，于常青

1大地構(gòu)造與動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，中國地質(zhì)科學(xué)院地質(zhì)研究所，北京 100037

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

1 引言

區(qū)域地球物理學(xué)旨在應(yīng)用先進(jìn)的物理場觀測技術(shù)，按規(guī)定網(wǎng)站密度觀測區(qū)域地球物理場，取得關(guān)于區(qū)域地殼上地幔組成結(jié)構(gòu)的信息，最終揭示區(qū)域地殼三維結(jié)構(gòu)與構(gòu)造形態(tài)，為資源勘探、地質(zhì)災(zāi)害預(yù)防及生態(tài)均衡提供科學(xué)依據(jù).由于區(qū)域地球物理場觀測網(wǎng)度均勻，不受巖芯分布及巖石出露的限制，而且信息來自各種深度的場源，可以根據(jù)物理學(xué)定律推算地殼三維結(jié)構(gòu)與構(gòu)造形態(tài)，因此區(qū)域地球物理學(xué)研究日益受到重視.

地球物理場包括重力場、磁場、反射地震波場、寬頻地震波場、地溫場等；觀測的方式有地面、航空和人造衛(wèi)星三種.航空磁場是觀測成本最低的地球物理場，但是它隨離開場源的距離的平方較快衰減，因此僅反映上地殼物質(zhì)磁化強(qiáng)度信息，信息量少.反射地震波場是觀測成本最高的地球物理場，雖然可以提供的信息量最多，但因要求太高的投入對大區(qū)域研究不實(shí)用.重力場是觀測成本較低的地球物理場，它隨離開場源的距離衰減，因此可反映全地殼物質(zhì)的密度信息，信息量較多.尤其是我國大部分地區(qū)都完成了1∶20萬區(qū)域地面重力調(diào)查，周邊地區(qū)也有衛(wèi)星重力場作參照，利用區(qū)域重力場的信息，揭示區(qū)域地殼三維密度結(jié)構(gòu)與構(gòu)造分區(qū)，就成為現(xiàn)今區(qū)域地球物理學(xué)研究的主攻方向.不過，要實(shí)現(xiàn)揭示區(qū)域地殼三維密度結(jié)構(gòu)與構(gòu)造的目標(biāo)，還必須研究系統(tǒng)的數(shù)據(jù)處理、信息提取和反演解釋方法.

經(jīng)過多年研究，我們把小波多尺度分析、表面刻痕分析以及位場頻率域解釋理論和密度擾動(dòng)反演方法有機(jī)地結(jié)合起來，形成了應(yīng)用區(qū)域重力場刻畫地殼三維密度結(jié)構(gòu)構(gòu)造的數(shù)據(jù)處理、反演解釋和信息提取的方法系統(tǒng)，分為區(qū)域重力場按場源分層、小波變換多尺度分析、場源分層深度及密度擾動(dòng)反演、分層刻痕分析和構(gòu)造邊界定位四個(gè)子系統(tǒng)，詳見圖1.這一方法系統(tǒng)簡稱為區(qū)域重力場多尺度刻痕分析方法，本文將扼要介紹其方法原理及應(yīng)用效果.

圖1 區(qū)域重力場多尺度刻痕分析方法框圖Fig.1 The flow－chart of the multi－scale scratch analysis for quantitative interpretation of regional gravity field

2 布格重力場小波變換的尺度場源深度轉(zhuǎn)換律

1990年以后，我們首先進(jìn)行了重力異常小波變換多尺度分析研究，取得了較好的效果.我們應(yīng)用Mallat塔式算法（Daubechies，1990；Mallat，1989），構(gòu)造出專門用于重力異常分解的小波基（侯遵澤和楊文采，1997；侯遵澤等1998）.二維離散小波變換的原理和區(qū)域重力場分解應(yīng)用已在前文作過介紹（侯遵澤和楊文采，1997，2011，2012；侯遵澤等，1998；楊文采等，2001），本文不再復(fù)述.區(qū)域重力場可表為二維單值實(shí)函數(shù)Δg（x，y），地殼三維密度結(jié)構(gòu)模型可表為三維單值實(shí)函數(shù)Δρ（x，y）.如果不附加約束性規(guī)律，僅靠地面二維數(shù)據(jù)是不可能取得三維密度結(jié)構(gòu)模型的.首先要附加的約束性規(guī)律稱為布格重力場小波變換的尺度－場源深度轉(zhuǎn)換律，推導(dǎo)如下.

區(qū)域重力場小波多尺度分析方法的理論依據(jù)之一，乃是布格重力異常場的尺度與場源埋藏深度呈現(xiàn)同步增大關(guān)聯(lián).換言之，場源埋藏深度越大，布格重力異常場的水平尺度就越大.例如，對于球形密度擾動(dòng)物體，重力異常表為（曾華霖，2005）

式中G為萬有引力常數(shù)，M為場源質(zhì)量；h為場源中心深度，x為觀測平面上測點(diǎn)離場源中心的水平距離.定義布格重力異常場的水平寬度w＝2x′，x′為異常從極大值衰減為1／2極大值的水平距離，W即為孤立異常的特征尺度.把x′代入（1）式可以導(dǎo)出由此可見，場源埋藏深度h越大，觀測平面上布格重力異常場的水平尺度W就越大，二者之間近似呈正比例.當(dāng)然，由于形狀不同，比例系數(shù)會(huì)有所變化，即α＝0.2～0.9.例如上地殼脈形侵入體可取α＝0.8.

當(dāng)多個(gè)源體產(chǎn)生的重力異常疊加時(shí)，疊加異常便不再具有類似于（2）式定義的特征尺度.例如，上地殼中心埋深4km的球體源的特征尺度為W＝6.13，下地殼中心埋深34km的球體源的特征尺度W＝52.1；它們疊加之后尺度可在6～52km的范圍變化，沒有特征尺度.小波多尺度分析方法利用小波基的特征尺度，把地面重力疊加異常場按尺度分解，使產(chǎn)生的異常集重新具有特征尺度.對上述例子，用尺度為6km，12km，24km和48km的小波基分解上述疊加異常，尺度為6km的一階小波細(xì)節(jié)D1提取出上地殼中心埋深h＝4km的球體源的異常場，尺度為48km的四階小波細(xì)節(jié)提取出下地殼中心埋深h＝34km的球體源的異常場，小波細(xì)節(jié)便重新恢復(fù)了特征尺度.因此，小波多尺度分析方法可用于刻畫地殼三維密度結(jié)構(gòu).

設(shè)均勻方格網(wǎng)重力異常場數(shù)據(jù)的網(wǎng)眼間距為Δ，進(jìn)行常規(guī)小波多尺度分析時(shí)，小波基取4～5個(gè)樣點(diǎn)，它有峰有谷，特征尺度約為峰寬之半，取為Δ／2，小波細(xì)節(jié)的尺度呈2的整數(shù)冪遞增.記小波細(xì)節(jié)階為n，n階小波細(xì)節(jié)特征尺度為

重力異常場小波多尺度分析要求小波細(xì)節(jié)特征尺度Ln和異常的特征尺度W匹配，即Ln＝W.代入（2）式有

代入（3）式有

（5）式表示重力異常源體埋深h與小波細(xì)節(jié)階次n的關(guān)系式，稱為重力場小波變換的尺度－源深度轉(zhuǎn)換律.給出取樣點(diǎn)間距Δ，不同比例尺的地面重力異常場數(shù)據(jù)在小波變換后由（3）式估算的小波細(xì)節(jié)特征尺度見表1.由此可見，取樣點(diǎn)間距大于20km的重力異常場數(shù)據(jù)因小波細(xì)節(jié)特征尺度過大，可能不適合用于刻畫精細(xì)地殼三維密度結(jié)構(gòu)（表1中下劃線數(shù)字）.取樣點(diǎn)間距小于1km的重力異常場數(shù)據(jù)因小波細(xì)節(jié)特征尺度過小，只適合用于刻畫上地殼三維密度結(jié)構(gòu).

表1 小波細(xì)節(jié)特征尺度Table 1 Characterized scale of wavelet details

綜上所述，對于適當(dāng)比例尺的區(qū)域重力數(shù)據(jù)而言，密度擾動(dòng)源埋藏深度和地面重力異常場的水平尺度呈近似正比關(guān)系，這就是重力場小波變換的尺度－源深度轉(zhuǎn)換律.因此，進(jìn)行地面重力異常場多尺度分解以后，小尺度的小波細(xì)節(jié)反映埋藏深度小的場源分布，而大尺度的小波細(xì)節(jié)反映埋藏深度大的場源分布.自20世紀(jì)90年代侯遵澤和楊文采提出并采用小波多尺度分析方法進(jìn)行重力場多尺度分解以來，重力異常的多尺度分析方法技術(shù)取得了良好的效果（侯遵澤和楊文采，1997，2011；侯遵澤等，1998；楊文采等，2001，2005，2014）.究其原因，正是由于小波細(xì)節(jié)反映不同埋藏深度的場源分布，可用于揭示地殼三維密度擾動(dòng)的緣故.

3 重力場小波細(xì)節(jié)的等效層深度公式

區(qū)域重力場可表為二維單值實(shí)函數(shù)Δg（x，y），多尺度分解以后，小波細(xì)節(jié)為三維單值實(shí)函數(shù)ΔD（x，y，n）的取樣數(shù)據(jù)集，而地殼三維密度結(jié)構(gòu)模型表為密度擾動(dòng)源埋藏深度的函數(shù).必需附加約束性規(guī)律，把小波細(xì)節(jié)ΔD（x，y，n）轉(zhuǎn)換成ΔD（x，y，h）的取樣數(shù)據(jù)集，取得對應(yīng)不同埋藏深度場源等效層的重力異常子集.在最理想的情況下，我們希望分解為對應(yīng)上、中、下地殼和上地幔頂部起伏的四個(gè)重力異常子集.這時(shí)要附加布格重力場小波細(xì)節(jié)的尺度－深度轉(zhuǎn)換約束，位場頻率域解釋理論為這個(gè)問題的解答提供了基礎(chǔ).對位場而言，“頻率域”與“波數(shù)域”是同義的，指的都是異常隨距離變化的頻度.

位場頻率域解釋理論證明（Bhimasankaram et al.，1977；楊文采等，1978，1979；楊文采，1985，1986a），對單個(gè)重力場源，重力異常對數(shù)振幅譜的斜率與源體頂面或重心埋藏深度成反比.以水平圓柱體為例，對數(shù)功率譜為LnS（ω）＝ －Ln（2π2Gr2ρ）ωh，其中G為萬有引力常數(shù)，ρ為源體與圍巖密度差，ω為圓頻率（波數(shù)）；r為半徑，h為中心埋藏深度，對數(shù)功率譜與頻率（波數(shù)）成反比，比例系數(shù)β＝Ln（2π2Gr2ρ）.對區(qū)域眾多重力場的場源，區(qū)域重力異常對數(shù)功率譜的斜率與源體重心平均埋藏深度成反比.因此，如果地殼中重力場的場源集中在幾個(gè)不同的深度段，區(qū)域重力場對數(shù)功率譜對應(yīng)不同頻段將呈現(xiàn)有可分辨的幾個(gè)不同斜率的線段，埋藏深度越大的源體對應(yīng)斜率越陡的線段.換言之，如果區(qū)域重力場對數(shù)功率譜對應(yīng)有可分辨的幾個(gè)不同斜率的線段，則表明它可以按場源深度分解為幾個(gè)重力異常子集，它們就是區(qū)域重力場小波細(xì)節(jié)或細(xì)節(jié)組合，對應(yīng)不同埋藏深度的場源等效層.因此，布格重力場小波細(xì)節(jié)的等效層深度公式可近似表為

式中hn為n階小波細(xì)節(jié)對應(yīng)的場源等效層埋藏深度，為功率譜最大斜率；β為比例系數(shù)，隨源體形狀有一定變化.

小波細(xì)節(jié)合并之后得到了N個(gè)異常子集，要把它們分別作二維Fourier變換，求它們的平均功率譜，即各個(gè)方向功率譜曲線的平均.由公式（6）可以計(jì)算出異常子集對應(yīng)的等效層平均埋藏深度hn，然后把小波細(xì)節(jié)集ΔD（x，y，n）轉(zhuǎn)換成ΔD（x，y，h）的取樣數(shù)據(jù)集.必需指出，區(qū)域重力場能否按場源深度分解，取決于重力場的頻率域?qū)傩裕⒎撬械貐^(qū)都可以作有意義的分解.筆者多年來的研究表明，中國大陸的區(qū)域重力場可以按場源深度分解為多個(gè)異常子集，詳見圖2.圖2a—2d分別給出了中國大陸、青藏區(qū)、揚(yáng)子區(qū)和華北區(qū)的區(qū)域重力場對數(shù)功率譜，區(qū)域的范圍和數(shù)據(jù)比例尺見表2.由圖2可見，四個(gè)區(qū)域重力場對數(shù)功率譜不同頻段呈現(xiàn)有可分辨的幾個(gè)不同斜率的線段，說明四個(gè)地區(qū)重力場都可以作有意義的分解.

區(qū)域重力場的頻率分析得出兩個(gè)結(jié)果：一是重力場可否用作有意義的分解，二是如果可以重力場應(yīng)分解為幾個(gè)等效層的異常子集.圖2中劃出了中國大陸、青藏區(qū)、揚(yáng)子區(qū)和華北區(qū)重力場對數(shù)功率譜不同斜率的線段，除青藏區(qū)地殼構(gòu)造比較復(fù)雜可分為5～6層外，其他試驗(yàn)區(qū)確定都可分為4～5層.區(qū)域重力場范圍和數(shù)據(jù)長度都有限，尤其在高頻段，僅進(jìn)行Fourier變換并不能取得準(zhǔn)確的功率譜，有限離散的區(qū)域位場功率譜的計(jì)算需要特殊的技巧（楊文采，1985，1986a）.注意在圖2b—2c中青藏區(qū)和揚(yáng)子區(qū)重力場對數(shù)功率譜中選擇不同斜率的線段時(shí)，躲開了最高頻段，因?yàn)檫@里含有明顯噪聲，不符合頻率增高重力場對數(shù)功率譜斜率明顯減少的規(guī)律.

表2 中國試驗(yàn)區(qū)的范圍和數(shù)據(jù)比例尺Table 2 Number of equivalent layers in different test area

圖2 區(qū)域重力場對數(shù)功率譜（a）中國大陸；（b）青藏區(qū)；（c）揚(yáng)子區(qū)；（d）華北區(qū).Fig.2 Logarithmic power spectra in some test regions（a）Chinese continent；（b）Qinghai－Tibetan area；（c）Yangtze craton area；（d）Northern China.

進(jìn)行重力異常多尺度分析要注意正確選擇小波變換最高階次m.區(qū)域重力場分解時(shí)選擇小波變換最高階次m的原則是：m必須大于根據(jù)功率譜特征決定的等效層數(shù)N，等效層數(shù)見表2最右列.以青藏區(qū)為例，根據(jù)布格重力異常的功率譜特征決定重力異常分解成5個(gè)層次，根據(jù)m＞N的原則，可將m選為8.然后，探討將8階小波變換多尺度分析應(yīng)用于重力異常的分解，可計(jì)算出小波細(xì)節(jié)D1—D8和8階小波逼近S8.接下來的問題是如何組合m個(gè)小波細(xì)節(jié)形成對應(yīng)N個(gè)等效層場源的重力異常子集.根據(jù)小波變換多尺度分析理論可知，小波細(xì)節(jié)的尺度是隨2m次方增加的.因此，小波細(xì)節(jié)D2的尺度為D1的2倍，小波細(xì)節(jié)D3的尺度為D1的4倍.由于小波細(xì)節(jié)D1和D2的尺度相差不大，通?？砂眩―1＋D2）作為地殼最淺等效層的重力異常對待；它可能代表盆地上部沉積層，也可能代表上地殼結(jié)晶基底，還要通過計(jì)算場源等效層的平均深度來決定.另一方面，如果某個(gè)小波細(xì)節(jié)的異常幅度很小，也要通過與鄰近小波細(xì)節(jié)合并使等效層場源的重力異常子集相對均衡.

圖3 （a）青藏區(qū)布格重力異常圖；（b）最深等效層異常S8，等效層深度74.8kmFig.3 （a）The Bourger gravity anomalies in the Qinghai－Tibet area；（b）The wavelet approximation S8for the deepest equivalent layer

圖4 （a）最淺等效層（D1＋D2）異常圖，等效層深度3.03km；（b）最淺等效層脊形化系數(shù)圖Fig.4 （a）The wavelet detail（D1＋D2）of the Bourger gravity anomalies；refers to the equivalent layer of average depth 3.03km；（b）The corresponding ridge coefficient image of the layer

可以青藏區(qū)為例，根據(jù)布格重力異常的功率譜特征決定重力異常分解成5個(gè)層次，而m＝8，計(jì)算出小波細(xì)節(jié)D1—D8和S8.小波細(xì)節(jié)合并為（D1＋D2）、（D3＋D4）、（D5＋D6）、（D7＋D8）和S8五個(gè)異常子集.青藏高原主要數(shù)據(jù)來自地面實(shí)測，高原外圍的部分?jǐn)?shù)據(jù)來自衛(wèi)星觀測，數(shù)據(jù)融合后的青藏區(qū)布格重力場見圖3a.小波多尺度分解和細(xì)節(jié)合并之后得到各等效層的重力異常見圖4a—7a.由平均功率譜可計(jì)算出各層的場源平均埋藏深度，最淺等效層異常（D1＋D2）為3.03km，淺等效層異常（D3＋D4）為12.83km，中等效層異常（D5＋D6）為19.52km；深等效層異常（D7＋D8）為37.93km，最深等效層異常S8為74.8km.最淺等效層、淺等效層和中等效層位于青藏區(qū)上地殼頂部、中部和底部；深等效層位于青藏區(qū)中地殼底部.最深等效層位于青藏區(qū)上地幔頂部，可能主要反映莫霍面起伏，由于本文只研究地殼內(nèi)部密度成像不予討論.如果要專門研究下地殼密度擾動(dòng)，還可以單獨(dú)研究異常子集D8，或者把D7和D8分別當(dāng)成獨(dú)立的等效層.這時(shí)平均埋藏深度計(jì)算結(jié)果為D7對應(yīng)34km，D8對應(yīng)52km.

圖5 （a）淺等效層（D3＋D4）異常圖，等效層深度12.83km；（b）淺等效層脊形化系數(shù)圖Fig.5 （a）The 3rd and 4th order wavelet details（D3＋D4）refers to the equivalent layer of average depth 12.83km；（b）The corresponding ridge coefficient image of the layer

圖6 （a）中等效層（D5＋D6）異常圖，等效層深度19.52km；（b）中等效層脊形化系數(shù)圖Fig.6 （a）The wavelet details（D5＋D6）refers to the equivalent layer of average depth 19.52km；（b）The corresponding ridge coefficient image of the layer

4 等效層的密度擾動(dòng)反演

在完成重力場小波變換多尺度分解之后，要由小波細(xì)節(jié)轉(zhuǎn)換為地殼三維密度結(jié)構(gòu)模型，用廣義線性反演方法求取各等效層的密度擾動(dòng)（Yang and Du，1993；楊文采，1986b，1997，1993a，1993b，1995）.位場正演問題是求下述Poisson方程邊值問題的解：

其中v（x，y，z＝0）為地表位場值，u（x，y，z）為下半空間任意點(diǎn)位場函數(shù)；ρ（x，y，z）為場源項(xiàng).由數(shù)學(xué)物理方法知，在地表（7）式的解為（楊文采，1997）

在本文中v用小波細(xì)節(jié)ΔD（x，y，n＝1，2，3，…）代入；ρ用等效層體密度擾動(dòng)Δρ＝ρ－ρ′代入；（8）式可用廣義線性反演求解（楊文采，1986b）.

等效層的密度擾動(dòng)指平面上變化的場源密度與等效層平均密度的差.除小波多尺度分析的小波細(xì)節(jié)數(shù)據(jù)外，廣義線性反演還要求輸入各等效層的平均密度ρ′和場源平均埋藏深度z′.各等效層的平均密度主要根據(jù)巖石物理性質(zhì)測定結(jié)果給出.以青藏區(qū)為例，如最淺等效層平均密度取為2.60g·cm－3，上地殼中部淺等效層密度取為2.67g·cm－3，上地殼底部中等效層密度2.81g·cm－3；中地殼底部深等效層密度2.95g·cm－3；上地幔頂部最深等效層密度3.20g·cm－3.由于反演出的各等效層的橫向相對密度擾動(dòng)圖外觀與小波細(xì)節(jié)異常（圖4a—7a）差別不大，本文不再羅列.

5 重力場多尺度刻痕分析

在反演計(jì)算了重力多尺度密度擾動(dòng)圖件之后，我們就獲得了地殼三維密度擾動(dòng)的映像.這些映像包含有大量地殼構(gòu)造的信息，能否根據(jù)現(xiàn)代信息學(xué)理論用計(jì)算機(jī)把地殼構(gòu)造的信息直接提取出來？要回答這個(gè)問題，首先要問你想提取什么地殼構(gòu)造信息？我們首先想到的是現(xiàn)存的地殼變形帶分布信息.

從物理學(xué)的觀點(diǎn)看，地殼變形帶是動(dòng)力地質(zhì)作用在地殼中的刻痕，它造成地殼巖石的密度在局部呈條帶狀變化.因此，不同強(qiáng)度的地殼變形帶是巖石圈地質(zhì)作用對應(yīng)地殼中長條形的密度異常帶，它們又在重力場中造成“刻痕”.這種刻痕的特征參數(shù)包括：重力異常梯度急劇變化，異常在平面上各向異性較強(qiáng)，各向異性長軸走向比較穩(wěn)定，等等.把隱含在區(qū)域重力場中有關(guān)刻痕的所有特征參數(shù)都提取出來，可在綜合分析的基礎(chǔ)上識(shí)別不同時(shí)期形成的地殼變形帶.計(jì)算了重力等效層多尺度密度擾動(dòng)之后，從其中提取地殼變形帶信息要有新的數(shù)學(xué)工具，這就是表面刻痕分析（孫艷云和楊文采，2014）.面對尺度變化巨大的自然界，應(yīng)用數(shù)學(xué)中的局部分析和全局分析永遠(yuǎn)是準(zhǔn)確認(rèn)識(shí)自然的重要工具之一（Riley，1974）.表面刻痕分析就是局部分析的一個(gè)分支.

區(qū)域位場函數(shù)g＝g（x，y）是一個(gè)二維連續(xù)可微函數(shù)，可視為表面函數(shù)的一種特殊類型.把區(qū)域重力場網(wǎng)格化成一個(gè)個(gè)表面元，每個(gè)表面元的信息可用表面元的譜矩表征（Longust－Higgins，1962；Nayak，1973；Yanagi et al.，2001；黃逸云，1984，1985；楊叔子等，2007）.最簡單的譜矩為二階譜矩，它的三個(gè)元m20、m02和m11有各自的物理含義.m20、m02分別表示位場表面x和y方向的斜率的方差；m11是x和y方向的斜率的協(xié)方差.位場輪廓是指一個(gè)切平面與位場表面的交線.若已知表面的二階譜矩為m20、m02和m11，位場任意表面元的一條輪廓在θ（與x軸方向的夾角）方向上的二階譜矩m2（θ）由（9）式唯一確定

m2（θ）＝m20cos2θ＋2m11cosθsinθ＋m02sin2θ，（9）m2（θ）等于θ方向上輪廓的斜率的方差.若m2max和m2min是m2（θ）的最大值和最小值，則其所在的方向總是相互垂直的，這兩個(gè)方向稱為刻痕主方向.m2min的那條輪廓曲線與x軸的夾角稱為主方向角θp，表示表面刻痕的走向

位場表面譜矩依賴于坐標(biāo)系，如果坐標(biāo)系旋轉(zhuǎn)，盡管此時(shí)表面性質(zhì)不變，表面的譜矩將會(huì)改變.因此，需要考慮一種與坐標(biāo)系旋轉(zhuǎn)無關(guān)的統(tǒng)計(jì)不變量.位場表面的二階譜矩可以用統(tǒng)計(jì)不變量M2和Δ2來表達(dá)，它們是（Longust－Higgins，1962；Nayak，1973；Yanagi et al.，2001；黃逸云，1984，1985；楊叔子等，2007）

統(tǒng)計(jì)不變量M2表示位場表面元的二維斜率的方差，可以表征布格重力異常的刻痕強(qiáng)弱，定義為刻痕的強(qiáng)度系數(shù)；試驗(yàn)表明，M2著重表現(xiàn)幅度大的表面刻痕.統(tǒng)計(jì)不變量Δ2主要與位場表面元的各向異性有關(guān)，表示表面元刻痕脊形化的程度，且總有Δ2≥0.試驗(yàn)表明，脊形化反映的是場在局部區(qū)段各向異性的程度，Δ2可以較全面地反映表面刻痕隨方向變化的形貌.

圖7 （a）深等效層（D7＋D8）異常圖，等效層深度37.93km；（b）深等效層脊形化系數(shù)圖L—低密度擾動(dòng)中心；H—高密度擾動(dòng)中心.Fig.7 （a）The wavelet details（D7＋D8）refers to the equivalent layer of average depth 37.93km；（b）The corresponding ridge coefficient image of the layerLetter“L”denotes low－density units，and letter“H”denotes high－density units.

圖8 多尺度密度擾動(dòng)邊界刻痕系數(shù)圖像（a）為中等效層，等效層深度為19.52km；（b）為深等效層，等效層深度為37.63km.L—低密度擾動(dòng)中心；H—高密度擾動(dòng)中心.Fig.8 The ridge－boundary coefficient image of the studied area（a）The middle layer of depth 19.52km；（b）The deep layer of depth 37.63km.Letter“L”denotes low－density units，and letter“H”denotes high－density units.

把隱含在區(qū)域重力場中有關(guān)刻痕的特征參數(shù)M2和Δ2提取出來后，便可在綜合它們特征的基礎(chǔ)上識(shí)別局部結(jié)構(gòu)為脊索形的地殼變形帶.為了對位場表面元的刻痕準(zhǔn)確定位，我們綜合這兩個(gè)統(tǒng)計(jì)不變量定義刻痕的脊形化系數(shù)Λ為

Λ值變化區(qū)間為［0，1］，體現(xiàn)了表面形態(tài)的各向異性特征；與統(tǒng)計(jì)不變量M2成反比，可增加對位場表面弱刻痕的識(shí)別能力.脊形化系數(shù)Λ反映的是面元內(nèi)場源密度的各向異性度的變化，Λ＝0表明面元內(nèi)場源密度各向同性，Λ＝1表明面元內(nèi)場源密度呈線性高度各向異性.場源密度呈線性高度各向異性就是地殼變形帶的典型特征.

提取區(qū)域位場函數(shù)刻痕信息就是要由地殼三維密度結(jié)構(gòu)模型計(jì)算等效層密度擾動(dòng)的脊形化系數(shù)Λ（x，y，n＝1，2，3，…）.孫艷云和楊文采（2014）詳細(xì)闡述了提取脊形化系數(shù)的理論與方法，并將這種方法應(yīng)用于中國大陸重力場刻痕信息識(shí)別，對比了不同類型的重力場刻痕與已知的地殼變形帶的對應(yīng)關(guān)系.與之前不同的是，本文不是用區(qū)域重力場，而是用區(qū)域重力場多尺度分解之后取得的異常子集數(shù)據(jù)，或者多個(gè)等效層密度擾動(dòng)數(shù)據(jù)作為刻痕信息識(shí)別的輸入數(shù)據(jù)，所以叫多尺度刻痕分析.以青藏區(qū)為例，把圖4a—7a中布格重力等效層或等效層密度擾動(dòng)作為輸入，求得的脊形化系數(shù)Λ分別見圖4b—7b.

6 邊界刻痕系數(shù)和構(gòu)造邊界劃分

迄今為止，大陸構(gòu)造單元的劃分都是大地構(gòu)造學(xué)家各自進(jìn)行的，不同大地構(gòu)造學(xué)家劃分的大陸構(gòu)造單元各自不同，缺乏客觀的標(biāo)準(zhǔn).如果能夠用數(shù)學(xué)方法把構(gòu)造單元信息從地球物理場中識(shí)別和提取出來，用計(jì)算機(jī)自動(dòng)識(shí)別構(gòu)造單元就成為可能.在提取區(qū)域位場函數(shù)的刻痕信息，即計(jì)算等效層密度擾動(dòng)的脊形化系數(shù)Λ之后，我們還可以進(jìn)一步提取密度擾動(dòng)的刻痕邊界信息，因?yàn)樗鼈円彩谴箨憳?gòu)造單元的邊界.

在提取構(gòu)造單元邊界信息時(shí)，邊界越窄單元?jiǎng)澐值亩ㄎ辉綔?zhǔn)確.脊形化系數(shù)Λ＞0.5的條帶常常較寬，為了準(zhǔn)確地對邊界定位，要再定義一個(gè)參數(shù)以取得窄的邊界線，這個(gè)參數(shù)稱為邊界刻痕系數(shù)MΛ，它對應(yīng)脊形化系數(shù)Λ較大區(qū)域的邊界，可以從脊形化系數(shù)Λ作進(jìn)一步增強(qiáng)處理計(jì)算出，詳見參考文獻(xiàn)（Longust－Higgins，1962；孫艷云和楊文采，2014；Riley，1974；黃逸云，1984，1985；楊叔子等，2007）.邊界刻痕系數(shù)MΛ提取連續(xù)性好的脊形化系數(shù)分布區(qū)段，并將它的邊界銳化定位.試驗(yàn)表明，在結(jié)晶巖分布區(qū)，脊形化系數(shù)Λ圖像反映顯生宙地殼變形帶，邊界刻痕系數(shù)圖像反映構(gòu)造單元邊界位置（孫艷云和楊文采，2014）.以青藏區(qū)為例，把圖6—7中等效層密度擾動(dòng)和脊形化系數(shù)Λ作為輸入，求得的中、深等效層邊界刻痕系數(shù)分別見圖8a，8b.

提取等效層密度擾動(dòng)邊界刻痕系數(shù)MΛ之后，不僅可以取得自動(dòng)劃分的等效層構(gòu)造單元邊界位置，還可以取得等效層構(gòu)造單元的密度擾動(dòng)屬性.以青藏區(qū)為例，首先研究圖7a中的等效層密度擾動(dòng)極性正負(fù)，在圖中用字母“L”標(biāo)明等效層低密度擾動(dòng)中心位置；用“H”標(biāo)明等效層高密度擾動(dòng)中心位置.對比圖7a和8a可見，自動(dòng)劃分的等效層構(gòu)造單元都與標(biāo)明“L”或“H”的等效層密度擾動(dòng)中心對應(yīng)，因此可知等效層構(gòu)造單元屬于高密度地塊還是低密度地塊.

以上介紹了應(yīng)用區(qū)域重力場刻劃地殼三維密度結(jié)構(gòu)和構(gòu)造單元分區(qū)的數(shù)據(jù)處理、反演解釋和信息提取的方法研究的主要結(jié)果，青藏區(qū)的實(shí)際應(yīng)用僅作結(jié)果圖示，未對其地質(zhì)含義進(jìn)行深入討論.由于青藏區(qū)的重力場多尺度刻痕的地質(zhì)含義仍然要繼續(xù)研究，今后再另外擬文發(fā)表.

7 結(jié)論

（1）把小波多尺度分析、表面刻痕分析以及位場頻率域解釋理論和反演方法有機(jī)地結(jié)合起來，可形成一個(gè)數(shù)據(jù)處理、反演解釋和信息提取的方法系統(tǒng)，這一方法系統(tǒng)簡稱為區(qū)域重力場多尺度刻痕分析方法，應(yīng)用于刻畫地殼分層的三維密度結(jié)構(gòu)、地殼變形帶和構(gòu)造單元分區(qū).

（2）多尺度刻痕分析包含頻率域重力場按場源分層、小波變換多尺度分析、場源分層深度及密度擾動(dòng)反演、分層刻痕分析和構(gòu)造邊界定位四個(gè)子系統(tǒng).本文扼要地介紹其方法原理及應(yīng)用效果.

（3）多尺度刻痕分析除應(yīng)用小波變換和表面譜矩之外，還運(yùn)用了地球重力場本身的規(guī)律，即以前導(dǎo)出的重力場小波細(xì)節(jié)的等效層深度公式和本文導(dǎo)出的重力場小波變換的深度－尺度轉(zhuǎn)換律.

（4）從地球物理探測到大地構(gòu)造學(xué)發(fā)現(xiàn)，是一個(gè)多學(xué)科綜合研究的探索過程.要取得重大研究成果，必須研發(fā)和組合來自不同學(xué)科的多個(gè)新方法技術(shù)，使多學(xué)科綜合研究有寬厚的理論支撐.本文介紹的四個(gè)子系統(tǒng)組合的理論支撐分別來自應(yīng)用數(shù)學(xué)、地球物理學(xué)和信息科學(xué).

Bhimasankaram V L S， Nagendra R，Rao S V S.1977.Interpretation of gravity anomalies due to finite inclined dikes using Fourier transformation.Geophysics，42（1）：51－59.

Daubechies I.1990.The wavelet transform，time frequency localization and signal analysis.IEEETransactiononInformationTheory，36（5）：961－1005.

Hou Z Z，Yang W C.1997.Wavelet transform and multiscale analysis of gravity field in China.ChineseJ.Geophys.（in Chinese），40（1）：85－95.

Hou Z Z，Yang W C，Liu J Q.1998.Multi－scale inversion of density contrast within the crust of China.ChineseJ.Geophys.（in Chinese），41（5）：642－651.

Hou Z Z，Yang W C.2011.Multi－scale inversion of density structure from gravity anomalies in Tarim Basin.Sci.inChina（Ser.D），54（3）：399－409.

Hou Z Z，Yang W C.2012.Applications of Wavelet Multi－Scale Analysis（in Chinese）.Beijing：Science Press.

Huang Y Y.1984.The characterization of three－dimensional Radom surface topography.JournalofZhejiangUniversity（in Chinese），18（2）：138－148.

Huang Y Y.1985.Geometrical interpretation and graphical solution of second order spectrum moments and statistical invariants for random surface characterization.JournalofZhejiangUniversity（in Chinese），19（6）：143－153.

Longust－Higgins M S.1962.The statistical geometry of random surfaces in hydrodynamic instability.／／13th Sympo.Applied Mathem.Amer.Math.Soc.，105－143.

Mallat S G.1989.Multifrequency channel decompositions of images and wavelet models.IEEETransactionsonAcoustics，Speech，andSignalProcessing，37（12）：2091－2110.

Nayak P R.1973.Some aspects of surface roughness measurement.Wear，26（2）：165－174.

Riley K F.1974.Mathematical Methods for the Physical Sciences：An Informal Treatment for Students of Physics and Engineering.New York：Cambridge Univ.Press.

Sun Y Y，Yang W C.2014.Recognizing and extracting the information of crustal deformation belts from the gravity field.ChineseJ.Geophys.（in Chinese），57（5）：1578－1587，doi：10.6038／cjg20140521.

Yanagi K，Hara S，Endoh T.2001.Summit identification of anisotropic surface texture and directionality assessment based on asperity tip geometry.InternationalJournalofMachine ToolsandManufacture，41（13）：1863－1871.

Yang S Z，Wu Y，Xuan J P，et al.2007.Identification on Surface Topography，in Time Series Analysis in Engineering Application （in Chinese）.Wuhan：Huazhong University of Science ＆ Technology Press.

Yang W C，Guo A Y，Xie Y Q.1978.Theory and methods for interpretation of gravity and magnetic anomalies in the frequency domain.BullInst.Geophys.Geochem.Prospect（in Chinese），（2）：134－179.

Yang W C.1979.Theorem of equivalence between the twodimensional and three－dimensional potential fields in the frequency domain.ActaGeophys.Sin.（in Chinese），22（1）：89－96.

Yang W C，Guo A Y，Xie Y Q，et al.1979.Interpretation of gravity anomalies in frequency domain （A）.ComputingTechniquesfor GeophysicalandGeochemicalExploration（in Chinese），1（1）：1－16.Yang W C.1985.The power spectrum analysis of the field data（A）.ComputingTechniquesforGeophysicalandGeochemical Exploration（in Chinese），7（3）：188－196.

Yang W C.1986a.The power spectrum analysis of the field data（B）.ComputingTechniquesforGeophysicalandGeochemical Exploration（in Chinese），8（1）：14－22.

Yang W C.1986b.A generalized inversion technique for potential field data processing.ChineseJ.Geophys.（ActaGeophys.Sin.）（in Chinese），29（3）：283－291.

Yang W C.1993a.Nonlinear chaotic inversion of seismic tracesⅠ：Theory and numerical experiments（Ⅰ）.ChineseJ.Geophys.（in Chinese），36（2）：222－232.

Yang W C.1993b.Nonlinear chaotic inversion of seismic tracesⅡ：Lyapunov exponents and attractors.ChineseJ.Geophys.（in Chinese），36（2）：376－387.

Yang W C，Du J Y.1993.Approaches to Solve Nonlinear Problems of the Seismic Tomography.／／Wei Y，Gu B L eds.Acoustic Imaging Volume 20.New York：Plenum Press，591－603.

Yang W C.1995.BG－Inverse scattering method for inversion of seismic wavefield.ChineseJ.Geophys.（in Chinese），38（3）：358－366.

Yang W C.1997.Theory and Methods in Geophysical Inversion（in Chinese）.Beijing：Geological Publishing House.

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－541，doi：10.3321／j.issn：0001－5733.2001.04.012.

Yang W C.2009.Tectonophysics in Paleo－Tethys Domain （in Chinese）.Beijing：Petroleum Industry Press.

Yang W C，Wang J L，Zhong H Z，et al.2012.Analysis of regional magnetic field and source structure in Tarim Basin.ChineseJ.Geophys.（in Chinese），55（4）：1278－1287，doi：10.6038／j.issn.0001－5733.2012.04.023.

Zeng H L.2005.Gravity Field and Gravity Exploration （in Chinese）.Beijing：Geological Publishing House.

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

侯遵澤，楊文采.1997.中國重力異常的小波變換與多尺度分析.地球物理學(xué)報(bào)，40（1）：85－95.

侯遵澤，楊文采，劉家琦.1998.中國大陸地殼密度差異多尺度反演.地球物理學(xué)報(bào)，41（5）：642－651.

侯遵澤，楊文采.2011.塔里木盆地多尺度重力場反演與密度結(jié)構(gòu).中國科學(xué)（D輯），41（1）：29－39.

侯遵澤，楊文采.2012.小波多尺度分析應(yīng)用.北京：科學(xué)出版社.

黃逸云.1984.三維隨機(jī)表面形貌的識(shí)別.浙江大學(xué)學(xué)報(bào)，18（2）：138－148.

黃逸云.1985.隨機(jī)表面二階譜矩和統(tǒng)計(jì)不變量的集合解釋及其圖解法.浙江大學(xué)學(xué)報(bào)，19（6）：143－153.

孫艷云，楊文采.2014.從重力場識(shí)別與提取地殼變形帶信息的方法.地 球 物 理 學(xué) 報(bào)，57（5）：1578－1587，doi：10.6038／cjg20140521.

楊叔子，吳雅，軒建平等.2007.時(shí)間序列分析的工程應(yīng)用（下）.武漢：華中科技大學(xué)出版社.

楊文采，郭愛纓，謝玉清.1978.重磁異常頻率域解釋的理論與方法.物化探研究報(bào)導(dǎo)，（2）：134－179.

楊文采.1979.二維與三維位場在頻率域的等價(jià)定理.地球物理學(xué)報(bào)，22（1）：89－96.

楊文采，郭愛纓，謝玉清等.1979.重磁異常在頻率域的解釋方法（上）.物化探計(jì)算技術(shù)，1（1）：1－16.

楊文采.1985.位場數(shù)據(jù)的頻譜分析技術(shù)（上）.物化探計(jì)算技術(shù)，7（3）：188－196.

楊文采.1986a.位場數(shù)據(jù)的頻譜分析技術(shù)（下）.物化探計(jì)算技術(shù)，8（1）：14－22.

楊文采.1986b.用于位場數(shù)據(jù)處理的廣義反演技術(shù).地球物理學(xué)報(bào)，29（3）：283－291.

楊文采.1993a.地震道非線性混沌反演—Ⅰ：理論與數(shù)值試驗(yàn).地球物理學(xué)報(bào)，36（2）：222－232.

楊文采.1993b.地震道非線性混沌反演—Ⅱ：關(guān)于Lyapunov指數(shù)與吸引子.地球物理學(xué)報(bào)，36（3）：376－387.

楊文采.1995.地震波場反演的B—G逆散射方法.地球物理學(xué)報(bào)，38（3）：358－366.

楊文采.1997.地球物理反演的理論與方法.北京：地質(zhì)出版社.

楊文采，施志群，侯遵澤等.2001.離散小波變換與重力異常多重分解.地球物理學(xué)報(bào)，44（4）：534－541，doi：10.3321／j.issn：0001－5733.2001.04.012.

楊文采.2009.東亞古特提斯域大地構(gòu)造物理學(xué).北京：石油工業(yè)出版社.

楊文采，王家林，鐘慧智等.2012.塔里木盆地航磁場分析與磁源體結(jié)構(gòu).地球物理學(xué)報(bào)，55（4）：1278－1287，doi：10.6038／j.issn.0001－5733.2012.04.023.

曾華霖.2005.重力場與重力勘探.北京：地質(zhì)出版社.