龍門(mén)山地區(qū)介質(zhì)密度的衛(wèi)星重力數(shù)據(jù)反演

王毅鵬，張永志，韓 鳴，程 冬，槐巖柯

(長(zhǎng)安大學(xué) 地質(zhì)工程與測(cè)繪學(xué)院，西安 710054)

0 引言

長(zhǎng)期以來(lái)，密度界面的反演是重力學(xué)研究的主要內(nèi)容之一，研究龍門(mén)山地區(qū)地下介質(zhì)密度結(jié)構(gòu)對(duì)于分析該地區(qū)的地震活動(dòng)具有重要意義。目前國(guó)內(nèi)外研究地殼內(nèi)部密度結(jié)構(gòu)的方法主要有兩種：①通過(guò)分析密度與速度的關(guān)系，將速度轉(zhuǎn)化為密度[1]；②通過(guò)研究重力異常，反演地殼內(nèi)部的密度界面。李媛等[2]使用地面實(shí)測(cè)重力數(shù)據(jù)，反演汶川地震之前龍門(mén)山地區(qū)地殼密度隨深度的變化規(guī)律，經(jīng)過(guò)分析得出該地區(qū)的介質(zhì)密度隨時(shí)間變化曲線顯示為由突變到平緩，同時(shí)介質(zhì)密度隨地殼深度的增加變化特點(diǎn)明顯；明圓圓等[3]利用魚(yú)群算法進(jìn)行密度異常的反演，當(dāng)密度參數(shù)較少時(shí)，該算法效果明顯，當(dāng)密度參數(shù)較多時(shí)，則此算法使用受限；柯小平等[4]采用直立長(zhǎng)方體模型結(jié)合遺傳算法反演了青藏高原東緣的三維密度分布。此外，由于位場(chǎng)的可疊加性，導(dǎo)致重力反演結(jié)果具有不唯一性。為限制重力反演結(jié)果的隨機(jī)性，加入地震測(cè)深數(shù)據(jù)約束是有效方法之一。唐新功等[5]在深地震測(cè)深數(shù)據(jù)資料約束下，使用GeoSoft軟件反演龍門(mén)山地區(qū)的沉積層、康拉德界面和莫霍面的深度分布，研究表明龍門(mén)山地區(qū)在南北地震帶兩側(cè)具有明顯不同的地殼結(jié)構(gòu)，該地區(qū)地殼的突變和不均勻性是導(dǎo)致地震活動(dòng)強(qiáng)烈的主要原因之一。這些研究結(jié)果為分析龍門(mén)山地區(qū)近年來(lái)所發(fā)生地震的成因提供了一定的依據(jù)，然而由于地震研究的極端復(fù)雜性和以前研究中密度數(shù)據(jù)量不足、空間分辨率低等條件限制，以及存在計(jì)算方法復(fù)雜、速度慢的問(wèn)題。高瑋等[6]研究證明粒子群優(yōu)化算法是研究地質(zhì)力學(xué)參數(shù)反演的一種有效方法, 同時(shí)Parker[7]、Oldenburg(1974)[8]提出的Parker-Oldenburg頻率域算法具有計(jì)算速度快速以及有效的特點(diǎn)；盧鵬羽等[9]將Parker-Oldenburg模型應(yīng)用于重力張量數(shù)據(jù)的密度界面反演中，發(fā)現(xiàn)可以明顯提高所反演密度界面的分辨率。因此筆者使用EIGEN-6C2超高階衛(wèi)星重力場(chǎng)模型，在頻率域采用粒子群算法，反演龍門(mén)山地區(qū)的介質(zhì)密度分布，并用于分析其地球物理學(xué)含義和研究。

1 Parker-Oldenburg頻域反演方法

通常情況下，提高正演計(jì)算速度是改進(jìn)反演速度的關(guān)鍵問(wèn)題。Parker[7]提出了重力異常正演計(jì)算的頻率域快速計(jì)算公式。Oldenburg[8]在Parker公式的基礎(chǔ)上，采用FFT(快速傅里葉變換)提高正演的計(jì)算速度，同時(shí)使用低通濾波器來(lái)保證算法迭代的收斂性。由于FFT技術(shù)的使用，Parker-Oldenburg模型的計(jì)算速度得到了很大地提高，成為頻域計(jì)算的主要方法。Parker-Oldenburg模型假設(shè)在x-z直角坐標(biāo)系中，重力異常用Δg表示，場(chǎng)源層上部邊界為z=0，下部邊界為z=h(x)，此邊界顯示界面的起伏。

由EIGEN-6C2重力場(chǎng)模型所計(jì)算的重力異常經(jīng)布格改正后可以得到布格重力異常，該重力異常歸因于地表到地球深部所有因密度不均勻引起的重力效應(yīng)的疊加，為了對(duì)地殼的密度分布特征進(jìn)行研究，需要從布格重力異常中扣除地殼以下物質(zhì)所引起的重力異常，由于莫霍面是地殼與地幔的分界面，因此也是密度發(fā)生突變的分界面。利用Parker頻率域方法可以正演莫霍面起伏引起的重力異常，Parker頻率域莫霍面起伏與地表重力異常關(guān)系為：

(1)

(2)

圖1 異常場(chǎng)源坐標(biāo)示意圖

Parker-Oldenburg模型算法執(zhí)行流程如下：

假定已知地層與下部介質(zhì)間的密度差ρ，參考面深度z0已知或給定，則可以使用式(2)進(jìn)行如下迭代計(jì)算：

1)給定界面起伏h(x)的初值。

2)將h(x)的初值代入公式(2)，計(jì)算該公式右端項(xiàng)的傅里葉變換。

3)式(2)右端項(xiàng)的傅里葉逆變換即為改進(jìn)的界面起伏h(x)。

4)判斷計(jì)算結(jié)果是否滿足收斂標(biāo)準(zhǔn)，或達(dá)到給定的最大迭代次數(shù)。如果滿足，則算法停止；否則轉(zhuǎn)步驟2)，并以本次迭代結(jié)果作為初值，繼續(xù)迭代計(jì)算。因重力異常的高頻成分會(huì)導(dǎo)致反演結(jié)果不穩(wěn)定，為保證算法迭代的收斂性，故對(duì)重力異常波數(shù)域的傅里葉變換結(jié)果進(jìn)行低通濾波處理[9]。

(3)

其中:a、b為濾波參數(shù)；波數(shù)k=1/λ，波長(zhǎng)λ的單位為km。

2 試驗(yàn)數(shù)據(jù)及處理方法

2.1 EIGEN-6C2超高階衛(wèi)星重力場(chǎng)模型

EIGEN-6C2重力場(chǎng)模型采用由地面重力數(shù)據(jù)、衛(wèi)星測(cè)高數(shù)據(jù)、衛(wèi)星重力數(shù)據(jù)聯(lián)合解算，是GFZ(德國(guó)地學(xué)中心)在2012年發(fā)布的最高階數(shù)達(dá)1949階的重力場(chǎng)模型。李琦等[11]研究表明使用新的EIGEN-6C2重力場(chǎng)模型作為參考模型，最終解算所得的似大地水準(zhǔn)面精度明顯優(yōu)于EGM2008模型所得解算結(jié)果,其精度甚至可以提高到3 cm。

圖2 龍門(mén)山地區(qū)自由空間重力異常

圖3 龍門(mén)山地區(qū)布格重力異常

2.2 布格重力異常及處理

筆者利用EIGEN-6C2重力場(chǎng)模型計(jì)算得到龍門(mén)山區(qū)域的自由空間重力異常,然后依據(jù)Crust1.0所提供的中國(guó)大陸沉積層模型作沉積層校正和布格校正，則可得到布格重力異常，其中布格校正所用地形高程數(shù)據(jù)來(lái)自美國(guó)海洋研究所 (http://topex.ucsd.edu/cgi-bin/get_data.cgi)，空間分辨率為2 km。

我們根據(jù)已有的地球物理探測(cè)結(jié)果,在使用 Parker-Oldenburg模型反演時(shí)，假設(shè)初始界面的地殼厚度為0，基于Crust1.0地殼模型的計(jì)算結(jié)果，將龍門(mén)山地區(qū)及附近區(qū)域初始界面密度取為2.11 g/cm3，參考面深度取z0=45 km，密度差初值ρ=0.45 g/cm3，同時(shí)對(duì)重力異常數(shù)據(jù)使用邊界錐形余弦濾波方法進(jìn)行處理[10]，濾波參數(shù)值b=0.005,a=0.002 5，分別為濾波頻率的上、下限值，經(jīng)6次迭代計(jì)算后，MSE=0.275 4，最終得到龍門(mén)山地區(qū)的自由空間重力異常，如圖2、圖3所示。

由圖2可知，龍門(mén)山地區(qū)自由空間重力異常變化劇烈，其變化范圍為：-250 mGal～380 mGal，呈現(xiàn)出西北區(qū)域自由空間重力異常較高，東南區(qū)域自由空間重力異常較低，并自西向東逐漸減小，其原因可能與該地區(qū)地形的變化有密切聯(lián)系。

由圖3可知，該區(qū)域布格重力異常的變化范圍為：-370 mGal～10 mGal，總體變化值較自由空間重力異常略低。西北區(qū)域布格重力異常低，東南區(qū)域布格重力異常則相對(duì)較高，同時(shí)布格重力異常值呈現(xiàn)自西向東逐漸增加的趨勢(shì)；此外，由于川西高原的布格重力異常變化范圍為：-370 mGal～-200 mGal；龍門(mén)山地區(qū)布格重力異常變化范圍為：-200 mGal～-40 mGal；四川盆地布格重力異常變化范圍為：-40mGal～-10 mGal,因此在龍門(mén)山地區(qū)附近形成了布格重力異常的高梯度帶。

3 PSO算法

3.1 PSO算法原理

粒子群算法是人工智能群優(yōu)化算法的一種,該算法在搜索尋找最優(yōu)解的過(guò)程中可以將最優(yōu)解保存，然后將最優(yōu)解分配給下一次搜索的粒子，它通過(guò)追隨當(dāng)前的局部最優(yōu)解來(lái)最終尋找到全局最優(yōu)解；此外，粒子群算法由于具有收斂速度快、所需調(diào)整的參數(shù)少、計(jì)算復(fù)雜度低等優(yōu)點(diǎn)被廣泛應(yīng)用于眾多的科研領(lǐng)域[12]。

粒子群算法假設(shè)在D維搜索空間中，粒子的群體由n個(gè)粒子所組成，全部粒子在飛行中都具有確定的速度，Xi=(xi1,xi2,…,xiD)是第i個(gè)粒子的位置，Vi=(vi1,vi2,…,viD)是第i個(gè)粒子的飛行速度，pbesti=(pbesti1,pbesti2,…,pbestiD)是第i個(gè)粒子所能搜索到的局部最優(yōu)位置，gbest=(gbest1,gbest2,…,gbestD)是整個(gè)粒子群所能搜索到的全局最優(yōu)位置，同時(shí)在算法的每次迭代中，每個(gè)粒子使用式(3)和式(4)來(lái)更新其速度和位置：

(4)

(5)

(5)

其中:F是目標(biāo)函數(shù);Δg0為重力異常觀測(cè)值；Δg為由密度參數(shù)所計(jì)算的重力異常值；N為觀測(cè)點(diǎn)個(gè)數(shù)。粒子群優(yōu)化算法執(zhí)行的流程描述如下：

1)確定待反演參數(shù)的搜索區(qū)域，同時(shí)設(shè)置其上下邊界，以及在搜索區(qū)域內(nèi)隨機(jī)產(chǎn)生一定數(shù)量的粒子，并且使每個(gè)粒子都具有隨機(jī)的位置和速度。

2)使用所設(shè)計(jì)的目標(biāo)函數(shù)F，計(jì)算每個(gè)粒子的最小適應(yīng)值。

3)將每個(gè)粒子當(dāng)前的局部最優(yōu)適應(yīng)值與已有全局適應(yīng)值進(jìn)行比對(duì)，如滿足條件，則將其更新為當(dāng)前粒子的最優(yōu)適應(yīng)值

4)繼續(xù)將每個(gè)粒子的當(dāng)前最優(yōu)值與群體的最優(yōu)值進(jìn)行對(duì)比，并更新當(dāng)前群體的全局最優(yōu)適應(yīng)值。

5)使用式(4)和式(5)迭代更新粒子的速度和位置。

6)使用已經(jīng)更新過(guò)的粒子速度和位置繼續(xù)計(jì)算目標(biāo)函數(shù)，如果所得適應(yīng)值滿足算法的精度要求，則算法結(jié)束，否則轉(zhuǎn)步驟2)繼續(xù)執(zhí)行算法。

3.2 粒子群算法反演結(jié)果

為驗(yàn)證本方法的有效性，我們使用粒子群算法反演計(jì)算龍門(mén)山地區(qū)介質(zhì)的密度參數(shù)，反演區(qū)域?yàn)?02°E～106°E和30°N～34°N，由于反演計(jì)算區(qū)域較大，故采用長(zhǎng)方體剖分單元法，將研究區(qū)域劃分成12′×12′共20×20個(gè)直立長(zhǎng)方體進(jìn)行計(jì)算,同時(shí)將每個(gè)長(zhǎng)方體內(nèi)部同一水平層的剩余密度設(shè)為常數(shù)，并與長(zhǎng)方體高度(地殼厚度)成函數(shù)關(guān)系，各水平層的地殼密度可以由其剩余密度加上初始界面密度得到。由于計(jì)算數(shù)據(jù)量仍然較大，普通計(jì)算機(jī)無(wú)法完成運(yùn)算，因此本文中的計(jì)算是在課題組內(nèi)存為64GB的高性能服務(wù)器上進(jìn)行PSO迭代運(yùn)算，并對(duì)每次運(yùn)算結(jié)果進(jìn)行修正，以獲得最優(yōu)解，最終運(yùn)行結(jié)果如表1，2所示，從以下兩表可以看出，該方法能得到一個(gè)穩(wěn)定的結(jié)果并接近真值。

表1 粒子群算法反演介質(zhì)密度參數(shù)與真值

表2 粒子群算法反演參數(shù)及結(jié)果

表3 龍門(mén)山地區(qū)介質(zhì)密度與Crust1.0對(duì)比

馮銳等[12]利用重力和地震資料，研究得出龍門(mén)山地區(qū)的地殼各層初始密度分別為：沉積層為2.3g/cm3、上地殼為2.67g/cm3、中地殼為2.80 g/cm3、下地殼為2.90 g/cm3[13-14]。本文以此為基礎(chǔ)，并將反演結(jié)果與Crust1.0模型所得結(jié)果進(jìn)行對(duì)比，結(jié)果如表3所示，以評(píng)估反演結(jié)果的正確性。

由表3可知，本文反演計(jì)算所得的結(jié)果與CRUST1.0模型值相比整體一致性較好[14]，但是反演結(jié)果顯示上地殼、下地殼密度都略低于模型值，這可能是高頻濾波時(shí)邊界效應(yīng)所導(dǎo)致的結(jié)果。

4 龍門(mén)山地區(qū)重力異常反演結(jié)果

4.1 莫霍面引起的重力異常及剩余重力異常

由圖4可知，龍門(mén)山地區(qū)莫霍面起伏引起的重力異?？傮w變化值為：-420 mGal～160 mGal，重力異常自西向東逐漸增加，這與該區(qū)域的地殼厚度有較為密切的聯(lián)系。通常地殼厚度越厚，莫霍面起伏引起的重力異常較低；地殼厚度越薄，則莫霍面起伏引起的重力異常較高。川西高原重力異常值為：-420 mGal～-320 mGal，龍門(mén)山地區(qū)重力異常為：-320 mGal～-200 mGal，四川盆地重力異常為：-200 mGal～-160 mGal，因此龍門(mén)山地區(qū)位于莫霍面起伏所引起的重力異常的高梯度帶上，其兩側(cè)重力異常值之差達(dá)到了120 mGal。

圖4 莫霍面起伏引起的重力異常

圖5 剩余重力異常

圖6 龍門(mén)山地區(qū)密度分布

龍門(mén)山地區(qū)的剩余重力異常，可以通過(guò)從區(qū)域布格重力異常中扣除莫霍面起伏所引起的重力異常而得到。因此由圖5可知，該區(qū)域剩余重力異常自西向東逐漸增加，以龍門(mén)山地區(qū)為分界，西北區(qū)域剩余重力異常較低，說(shuō)明該區(qū)域介質(zhì)密度相對(duì)較小，這可能是青藏高原地下物質(zhì)和能量的“東流”的結(jié)果；東南區(qū)域則剩余重力異常較高，說(shuō)明該區(qū)域介質(zhì)密度較大，這可能和四川盆地的克拉通地殼有關(guān)[15-16]。

4.2 地下介質(zhì)密度

采用Parker-Oldenburg模型以及粒子群算法反演龍門(mén)山地區(qū)的地下介質(zhì)密度結(jié)構(gòu)，該區(qū)域地形起伏較大，包括川西高原、龍門(mén)山地區(qū)和四川盆地，反演所得介質(zhì)密度分布如圖6所示。

圖6顯示了龍門(mén)山地區(qū)密度分布隨地殼深度變化的特征：整體而言，不同深度的密度分布與區(qū)域剩余重力異常呈現(xiàn)對(duì)應(yīng)關(guān)系：高密度區(qū)對(duì)應(yīng)的區(qū)域剩余重力異常高；低密度區(qū)對(duì)應(yīng)的剩余重力異常低；不同深度層的介質(zhì)密度均表現(xiàn)為西北區(qū)域低、東南區(qū)域高，并且表現(xiàn)出自西向東逐漸增加的特點(diǎn)；顯示龍門(mén)山地區(qū)東西兩側(cè)密度差異較大；此外在圖6中，深度為10 km～20 km和30 km～40 km的橫向分辨率差異較小這可能是由于地下介質(zhì)密度之間中間層的連續(xù)性所導(dǎo)致的結(jié)果，反映出該區(qū)域劇烈的地下構(gòu)造運(yùn)動(dòng)特征[17]。龍門(mén)山地區(qū)地下介質(zhì)密度值隨深度變化對(duì)比關(guān)系如表4所示。

由表4分析可知，龍門(mén)山地區(qū)以及相鄰區(qū)域不同深度的地殼密度分布各異，介質(zhì)密度隨深度加深而增加；在同一深度，其密度大小分布也不相同，呈現(xiàn)出自西向東逐漸增加的特點(diǎn)。

表4 密度隨深度變化關(guān)系對(duì)比表

5 結(jié)論

利用GFZ提供的EIGEN-6C2超高階衛(wèi)星重力場(chǎng)模型，計(jì)算得到龍門(mén)山地區(qū)的自由空間重力異常，并結(jié)合界面反演確立龍門(mén)山地區(qū)莫霍面引起的區(qū)域異常，進(jìn)而獲取該區(qū)域的剩余重力異常，然后對(duì)該異常采用粒子群算法進(jìn)行反演，最終得到該區(qū)域的三維地殼密度分布。通過(guò)分析模型反演的結(jié)果可知，龍門(mén)山地區(qū)位于密度變化的梯度帶，東西兩側(cè)密度差值較大，屬于密度變化強(qiáng)烈的區(qū)域，地下不同深度層的密度均具有西北區(qū)域低，東南區(qū)域高，以及呈現(xiàn)介質(zhì)密度由西向東逐漸增加的特點(diǎn)。經(jīng)過(guò)將本文所得結(jié)果與Crust1.0模型提供的密度值進(jìn)行對(duì)比分析，結(jié)果表明本文所得密度分布與Crust1.0模型值的一致性較好,因此利用衛(wèi)星重力數(shù)據(jù)使用Parker-Oldenburg 模型結(jié)合粒子群算法,反演三維地殼密度分布具有快速、有效和方法相對(duì)簡(jiǎn)單的特點(diǎn)，對(duì)于深入研究地殼介質(zhì)密度有一定的參考價(jià)值。

致謝：

感謝編輯部的大力支持和匿名審稿專家的寶貴建議。