小波多尺度分析在青藏高原東緣重力場的應用

劉 威, 王緒本, 王向鵬, 張 翔

(成都理工大學 地球勘探與信息技術教育部重點實驗室，成都 610059)

0 引言

青藏高原東緣是青藏高原與揚子板塊接壤的過渡帶,受印度板塊北向地應力,青藏高原內物質在應力作用下向川滇區(qū)塊運移[1]。該區(qū)域發(fā)生強烈變形,發(fā)育多條大型南北向深大斷裂帶,是我國大陸構造運動最為強烈的區(qū)域之一。近代以來，該區(qū)發(fā)生多起特大地震事件，是國內、外學者研究的熱點區(qū)域。學者們先后利用不同方法對該區(qū)域構造背景成像及物質運移動力學進行相應的研究，丁志峰等[1]、王椿鏞等[2]、白志明等[3]利用地震反射剖面資料對青藏高原東緣及鄰近區(qū)域上地幔地殼反射結構及P波波速結構進行分析研究，探討了地殼厚度變化規(guī)律及深部構造運移背景；孫潔等[4]、王緒本等[5]、白登海等[6]利用大地電磁剖面針對青藏高原東緣各構造地塊電性結構特征進行研究，揭示了深部電性結構及提出了下地殼流模型，為動力學研究提供了科學依據；蔣福珍等[7]、孟小紅等[8]、楊文采等[9]利用重磁資料對東緣及鄰近區(qū)域重磁異常進行分析，基于重磁場研究了橫向及縱向密度變化特征，并對鄰區(qū)內主要構造背景進行綜合分析。前人在此也取得了豐富的研究成果以及科學認知，但由于方法及技術的應用不同，產生了不同的看法，對此筆者利用小波多尺度分析方法對衛(wèi)星重力數據進行位場分離處理，結合徑向平均功率譜對不同界面的重力異常以及莫霍面起伏進行分析。重力異常是地球內部不均勻密度物質引起的不同深度重力場的疊加場效應，通過實現有效的位場分離，提取不同深度、不同介質因密度不均而引起的重力場，楊文采等[10-14]利用小波多尺度分析進行了有效的位場分離分析研究，并取得了較好地應用。

本文數據來源于WGM2012地球重力場模型(地形和重力資料來自美國加州大學圣迭戈分校斯克里普斯海洋研究所http://topex.ucsd.edu/)，利用多尺度小波分析提取不同尺度的重力異常，采用徑向功率譜分析提取異常的位置，進行莫霍面界面起伏反演，從而可以分析深部構造[6-9]。

1 地質構造及重力場數據

在四川盆地的阻擋作用下,青藏高原東向運移的過程中,在東緣形成高原的強邊界,在川滇地區(qū)形成弱邊界[15],使得青藏高原構造背景復雜,地震活動頻繁，研究區(qū)范圍為98°E～104°E，26°N～34°N，區(qū)域地形起伏較大，向東地勢急劇下降，構造復雜，自由空氣重力異常與地形起伏相關。研究區(qū)內分布眾多斷裂帶，分別為龍門山斷裂，鮮水河斷裂，哀牢山-紅河斷裂，瀾滄江斷裂，以及對青藏高原主要隆升起控制作用的金沙江縫合帶，班公湖-怒江縫合帶，雅魯藏布江縫合帶，該研究區(qū)地殼厚度顯著增大，具有較明顯的重力梯度帶。

這里使用的WGM2012地球重力場模型，是基于EGM2008重力場模型和ETOPO1模型基礎的，網格數據為2’×2’，球諧系數擴展至2159階。根據位場理論，利用引力位來計算重力場，地面重力異常(自由空氣重力異常)與擾動位為一階偏導的關系，自由空氣異常表示為式(1)。

(1)

圖1 自由空氣重力異常圖

Δg中間層=-2πfMρΔh=-0.1118Δh

(2)

(3)

Δg布格=Δg自由+Δg中間層+Δg地形

(4)

式中：ρ為地表巖石密度，2.67g/cm3; Δx、Δy為方格網中節(jié)點網格距;rij為測點與節(jié)點的間距，hij為測點與節(jié)點的高程差，cij為數值積分系數，改正區(qū)內點系數cij為1，改正區(qū)內頂點系數cij為0.75，改正區(qū)邊緣點系數cij為0.5，改正區(qū)外頂點系數cij為0.25。進行正常場改正和地形改正后的重力布格重力異常，可以較為精細地揭示研究區(qū)的深部結構，重力布格異常圖如圖2所示。

圖2 布格重力異常圖

經自由層校正后的重力布格異常，去除了地表起伏帶來的異常影響，能較好地反映地下深部構造，研究區(qū)布格重力異常的異常值范圍為-100 mGal～-680 mGal,重力異常值呈現自西北向西南逐漸增大的趨勢，且異常呈串珠狀分布，與研究區(qū)內斷裂及地體分布有一定關系。其負低異常分布在揚子地塊，川滇菱形塊體，騰沖地塊；高負異常分布在青藏高原的巴顏喀拉塊體，松潘甘孜塊體以及川滇北部部分區(qū)塊，且存在條帶狀負異常。

2 方法原理

2.1 二維小波多尺度分析

地球重力場受地下各種密度不均勻介質體的影響，布格重力異常中包含了不同深度，不同尺度的各種重力場源。在研究不同尺度的重力場時，要對這種復雜的場進行場分離，分別得到不同研究尺度的場。小波多尺度分析可以將信號按各種不同頻率進行分解，可以在任意細節(jié)上體現不同細節(jié)的地球物理意義[10,16]。為了研究不同地質體引起的對應重力場特征，可以利用小波多尺度分析方法分離出不同尺度的地質體在橫向和縱向分布的重力異常場[8,17]。

根據小波分析理論，假設二維重力異常場為：

Δg(x,y)=f(x,y)

(5)

重力分解表達式可簡化為：

Δg(x,y)=ANG+D1G+D2G+…DNG

例如：在《黃鶴樓送孟浩然之廣陵》相關內容教學中，小學生由于自身經歷的局限性，難以真正理解文章中描述的內容，從而難以融入自身的情感，朗讀質量也不高。針對這種情況，教師可以利用多媒體工具，為學生盡可能還原當時的情景，通過視頻、音頻、圖片等直觀的刺激，讓學生有更加直觀的認識，更好地融入到教學情境中進行朗讀，自身的情感也能更好地融入，朗讀興趣與積極性也會隨之提升。

(6)

式中:D1～DNG為異常的1階到N階小波細節(jié)；ANG為異常N階小波逼近。圖3為小波分解示意圖。

圖3 小波分解示意圖

2.2 徑向平均功率譜

功率譜分析是由 Spector等[18]最早提出來的一種重磁場解析處理方法，用來估測地下異常體平均深度。經傅里葉變換后的重力異常對其頻譜進行對數功率譜分析，在極坐標系下重磁異常的徑向平均對數功率譜曲線圖可以反映異常場源的平均深度，定量估計重磁異常的場源深度，徑向功率譜及場源埋深深度表達式為：

E(r1)=A2(r1)

(7)

(8)

式中：E(r1)、lnE(r1)為異常波譜的功率譜與徑向對數功率譜；A(r1)為頻譜幅值；r為圓波數；h為場源埋深。徑向對數功率譜值與徑向頻率線性相關，線性直線斜率常用來推斷估計不同尺度重磁異常的場源平均深度。

2.3 Parker迭代反演

Parker-Oldenburg位場迭代公式常用于重磁反演計算中，由于其采用了快速傅里葉變化和反演迭代，能計算物性橫向變化的連續(xù)界面，計算的速度較快。馮銳等[25]對Parker-Oldenburg迭代反演方法進行詳細的敘述，通過反演迭代式(10)計算迭代修正量Δhij，利用Parker正演式(9)計算下界面引起的重力異常，當前后兩次迭代結果差值小于閾值時，截止迭代，輸出界面深度。

(9)

(10)

(11)

式中:F[h(r)]、F[Δg(r)]分別是介質界面深度h(r)和重力異常Δg的傅里葉變換;z0為界面平均參考深度;ρ為界面密度差;r為位置矢量;k為波數;G為萬有引力常數;MN為總計算點;l為迭代次數;∈為設置閾值。

3 重力模型試算

3.1 建立正演模型

筆者建立了由不同尺寸長方體組合的三層地質異常模型，用5個同深度不同大小的小長方體模擬淺層地質體；用2個中等尺寸的長方體模擬中部穿插構造單元，用兩個相接不同埋深的大長方體模擬深部區(qū)域背景。模擬背景區(qū)域為32 km×32 km，模擬示意圖及模型參數分別如圖4和表1所示。

圖4 模型示意圖

表1 正演模型參數

3.2 模型多尺度分析

利用不同尺度長方體組合模擬地下結構分布，組合的模型正演重力異常如圖5所示。在進行多尺度小波分解的過程中，小波基函數的選擇對小波分解的效果有一定的影響，通常選取的小波基函數不僅要滿足一維分解重構的準確性，而且還要能體現二維數據成像的優(yōu)良性[19],正交小波可以保證低階細節(jié)不變形，對稱性影響小波分解后的重構。筆者選取了具有正交性和對稱性的coif3函數，即能夠突出一維分解重構的準確性，還能體現二維數據成像的優(yōu)良性[20-21]。圖6為coif3的尺度函數和小波函數圖像，對模型組合重力異常進行六階多尺度分析，多尺度分析結果如圖7所示。

圖5 模型正演重力異常

圖6 coif3小波基函數圖

通過比較正演模型與小波多尺度分析結果，小波多尺度分析的各階細節(jié)能清晰反應不同深度尺寸的重力異常體分布。如圖7所示，1階～5階小波細節(jié)組合可以清晰看到5個不同位置的異常體，形態(tài)與正演形態(tài)有差別，但位置較一致，較好地體現了局部異常；小波6階細節(jié)(圖7(c))體現出了中部穿插異常體的位置，而小波6階逼近較好地分離出了深部背景場，將1階～6階小波細節(jié)進行結合，可以分離出中淺部異常體的分布。因此小波多尺度分析既可以有效對區(qū)域場進行異常分離，還能針對不同深度異常體分布進行有效分離，有助于分析地下介質的重力場源情況。

圖7 模型小波分析

3.3 徑向功率譜分析

經小波多尺度分析后的各不同尺度場源異常，利用徑向平均對數功率譜對其進行深度估計。對計算后的徑向對數功率譜數據的最低波數段進行線性擬合，由式(8)可知，擬合直線二分之一斜率為場源的平均深度值。1階～5階小波細節(jié)異常圖的波譜圖8(a)擬合直線斜率為4.86，反應了深度約為2.43 km的淺層長方體；6階小波細節(jié)異常圖的波譜圖8(b)擬合直線斜率為13.88，反應了深度約為6.94 km的中部長方體；6階小波逼近異常圖的波譜圖8(c)擬合直線斜率為36.34，反應了深度約為18.17 km的深部背景長方體。與模型數據對比，徑向對數功率譜估計的平均深度較準確。

圖8 模型徑向功率譜

4 多尺度分析結果

基于式(6)對青藏高原東緣的重力布格異常進行多尺度小波分解，小波細節(jié)(D1G～D5G)和小波逼近(A5G)結果如圖9所示。對圖9(a)～圖9(e)小波細節(jié)分別進行徑向平均功率譜計算分析，分別如圖10(a)～圖10(e)所示。徑向平均功率譜計算結果表明，1階和2階細節(jié)的場源深度分別為4 km和11 km，反映了上地殼沉積物質密度不均勻體，3階和4階的場源深度分別為26 km和37 km，反映了中地殼和下地殼物質密度不均勻性，5階小波的場源深度為54 km，反映了莫霍面的起伏。

圖10 研究區(qū)多尺度小波分析重力異常徑向平均功率譜圖

小波多尺度分解的1階小波細節(jié)(圖9(a))和2階小波細節(jié)(圖9(b))在松潘甘孜地塊、川滇地塊、騰沖地塊異常比較復雜且劇烈，表明淺部存在高密度體,橫向上受構造影響具有不均勻性。在三江區(qū)域內顯示多條南北向條帶狀負異常，條帶狀負異常的極大值與該處的哀牢山-紅河斷裂，瀾滄江斷裂等南北向斷裂吻合性較好，鮮水河斷裂在1階和2階小波細節(jié)中的負異常值吻合性也比較好，揭示了該區(qū)域沉積物質高頻異常的特征，表明了該區(qū)復雜的構造背景。

隨著小波分解階數的增大，分解的重力異常的場源深度也隨之變大，重力異常逐漸變得緩和與平滑。地殼內部不同物質密度分布不均會導致重力場出現大小不一的正負異常，通常情況下，正異常與正密度異常相關，負異常與負密度異常相關。3階小波細節(jié)(圖9(c))和4階小波細節(jié)(圖9(d))主要反映青藏高原東緣區(qū)域中地殼和下地殼物質密度的不均一性，如深大斷裂引起的密度差異導致的異常，龍門山斷裂附近存在正負異常錯綜復雜的條帶，對研究區(qū)重力場的分區(qū)特征有所反映，可揭示青藏高原東緣深部存在不同性質的塊體之間的碰撞與融合。

5階小波細節(jié)(圖9(e))的場源深度與研究區(qū)內莫霍面的平均深度接近，異常值從西至東呈現“高低高”現象，在局部區(qū)域出現團狀負異常和團狀正異常，正異常主要分布在松潘甘孜地塊，四川盆地所在的揚子地塊，巴彥克拉西側，表現了剛性塊體特征。負異常主要分布于川滇地塊。正異常形成的原因有可能是上地幔高密度熔融物質沿斷裂通道上涌，負異常可能是深部軟、熱物質向東運移所產生的。通過分析不同異常產生的原因，可以幫助我們更好地去理解深部塊體間的物質運移以及深部大尺度的異常源。

圖9 研究區(qū)多尺度小波分析圖

小波5階小波逼近(圖5(g))可反映深部莫霍面起伏引起的重力異常，為研究區(qū)內區(qū)域場，1階～5階小波細節(jié)之和(圖9(f))可作為場分離后的局部場。相比布格異常圖，消除了區(qū)域背景場后的局部異常，異常值分布有所變化，變化呈現中間低兩邊高的異常趨勢，相比重力布格異常的梯度變化而言，局部異常更趨向于片區(qū)狀以及串珠狀分布，在橫向上變化明顯，體現了莫霍面以上物質分布不均而引起的重力場變化，在川滇地塊與揚子地塊，甘孜地塊與揚子地塊的接連處重力異常皆有明顯變化，且存在與川西地區(qū)地震事件具有緊密聯系的鮮水河斷裂和龍門山斷裂，有可能是因為青藏高原深部熱物質向東逃逸過程中延斷裂通道向上運移且在不同深度環(huán)境下發(fā)生聚集，形成局部高密度異常體，在研究區(qū)域西南邊出現一系列重力異常的過渡帶、梯度帶、突變，呈現局部狹窄的北西向異常帶，極有可能是此處俯沖斷裂所引起的成巖物質所導致的。

5 莫霍面結構分析

莫霍面為地球內部與地幔的分界面,是引起區(qū)域重力異常的深部主要分界面,通過對莫霍面的確定,可以為認識地殼與地幔的構造發(fā)育以及構造運動提供參考。5階小波逼近見圖5(g)，利用徑向功率譜計算的平均場源深度為52 km，與測深地震識別的莫霍面深度相似[22]，異常的趨勢與莫霍面的趨勢具有相似性，可以反映莫霍面起伏形態(tài)。地球內部物質密度分布影響著地球外部重力場，地殼與地幔間的密度差為重力資料確定莫霍面的位置提供了良好的基礎。綜合前人在該區(qū)域及附近區(qū)域的研究成果及地質背景資料[9,23]，利用5階小波逼近(圖9(g))采用Parker-Oldenburg界面反演方法，對研究區(qū)莫霍面深度進行迭代反演，選取重力基本場為5階小波逼近，參考深度z0=50 km，殼幔密度差的基本參數，反演結果與CRUST1.0地殼模型一致性較好[24]，莫霍面界面埋深范圍約為37 km～59 km，莫霍面深度從西北至東南逐漸抬升，地殼厚度逐漸減薄(圖11)。

圖11 研究區(qū)莫霍面深度圖

基于前人研究的深部地震以及大地電磁探測資料背景，將研究區(qū)域內的地殼結構分為沉積層、上地殼、中地殼、下地殼4 層，地殼各層及上地幔的初始密度分別為：沉積層為2.3 g/cm3、上地殼為2.60 g/cm3、中地殼為2.70 g/cm3、下地殼為2.80 g/cm3，上地幔為3.30 g/cm3[25]。使用Oasis Montaj平臺中的GM-SYS模塊對東西走向，橫跨川滇地塊和楊子地塊的剖面AB進行人機交互重力密度擬合反演，調整界面深度以及模型密度參數，擬合效果較好(圖12)。松潘甘孜，巴彥克拉地塊的莫霍面深度在52 km以上，川滇地塊莫霍面深度主要分布在47 km左右，與地震測深資料獲得的地殼厚度圖形態(tài)且厚度相一致，展現為西北至東南逐漸變淺的趨勢，且青藏塊體與川滇地塊和揚子地塊深度變化較大，形成了環(huán)青藏高原莫霍深淺變化的梯級帶狀。龍門山斷裂附近莫霍面深度差較大，在該區(qū)域發(fā)生多起特大地震災害，青藏高原深部物質向東逃逸過程中在剛性四川盆地的阻擋作用下，物質運移能量發(fā)生聚集且延龍門山構造進行釋放，從而消耗掉了板塊之間巨大的碰撞應力。川滇塊體南部莫霍深度帶較為輕緩地向東南漸變，與下地殼流方向基本一致，軟弱的下地殼物質延南向、東南向流動，導致地殼厚度大于剛性體的四川盆地地殼厚度。青藏高原東緣地殼的東西部厚度差異，可以合理地解釋青藏高原物質東流的channel flow模式，四川盆地剛性塊體阻擋了高原物質的東流，高原物質在此堆積擠壓，形成了龍門山局部擠壓推覆構造帶，形成的推覆構造帶吸收了來自西部青藏高原的能量，這也是青藏高原東緣地震能量的主要來源。

圖12 AB剖面密度模型擬合

6 結論

1)使用WGM2012地球重力場模型，解譯計算后的重力數據基于多尺度小波分析應用于青藏高原東緣，結合徑向平均功率譜計算不同尺度小波細節(jié)所揭示的地下場源平均深度，參考前人研究成果以及各種地質資料，分析不同尺度的異常與地質背景間的對應。

2)分解后的1階小波細節(jié)反映了上部沉積層物質密度橫向不均勻變化，3階～4階小波細節(jié)揭示了地殼內部物質在深部動力學影響下，物質縱向變化及橫向擴展而導致的不均勻性；5階小波細節(jié)及5階逼近反映了莫霍面起伏引起的區(qū)域場異常變化，能夠揭示地殼厚度總體性變化。

3)采用Parker-Oldenburg界面反演進行重力莫霍面反演，并結合深部地震研究進行重力剖面擬合驗證，反演的莫霍面深度與測深地震所探測的地殼厚度趨勢一致，呈梯度向東南減薄，間接說明在推覆擠壓的過程中深部物質向東逃逸，與剛性塊體接觸為地應力集中及釋放提供場所，導致龍門山等地活動構造激活，發(fā)生多起重大自然災害。