龍門山地區介質密度的衛星重力數據反演

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

(長安大學 地質工程與測繪學院，西安 710054)

0 引言

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

1 Parker-Oldenburg頻域反演方法

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

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

(1)

(2)

圖1 異常場源坐標示意圖

Parker-Oldenburg模型算法執行流程如下：

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

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

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

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

4)判斷計算結果是否滿足收斂標準，或達到給定的最大迭代次數。如果滿足，則算法停止；否則轉步驟2)，并以本次迭代結果作為初值，繼續迭代計算。因重力異常的高頻成分會導致反演結果不穩定，為保證算法迭代的收斂性，故對重力異常波數域的傅里葉變換結果進行低通濾波處理[9]。

(3)

其中:a、b為濾波參數；波數k=1/λ，波長λ的單位為km。

2 試驗數據及處理方法

2.1 EIGEN-6C2超高階衛星重力場模型

EIGEN-6C2重力場模型采用由地面重力數據、衛星測高數據、衛星重力數據聯合解算，是GFZ(德國地學中心)在2012年發布的最高階數達1949階的重力場模型。李琦等[11]研究表明使用新的EIGEN-6C2重力場模型作為參考模型，最終解算所得的似大地水準面精度明顯優于EGM2008模型所得解算結果,其精度甚至可以提高到3 cm。

圖2 龍門山地區自由空間重力異常

圖3 龍門山地區布格重力異常

2.2 布格重力異常及處理

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

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

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

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

3 PSO算法

3.1 PSO算法原理

粒子群算法是人工智能群優化算法的一種,該算法在搜索尋找最優解的過程中可以將最優解保存，然后將最優解分配給下一次搜索的粒子，它通過追隨當前的局部最優解來最終尋找到全局最優解；此外，粒子群算法由于具有收斂速度快、所需調整的參數少、計算復雜度低等優點被廣泛應用于眾多的科研領域[12]。

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

(4)

(5)

(5)

其中:F是目標函數;Δg0為重力異常觀測值；Δg為由密度參數所計算的重力異常值；N為觀測點個數。粒子群優化算法執行的流程描述如下：

1)確定待反演參數的搜索區域，同時設置其上下邊界，以及在搜索區域內隨機產生一定數量的粒子，并且使每個粒子都具有隨機的位置和速度。

2)使用所設計的目標函數F，計算每個粒子的最小適應值。

3)將每個粒子當前的局部最優適應值與已有全局適應值進行比對，如滿足條件，則將其更新為當前粒子的最優適應值

4)繼續將每個粒子的當前最優值與群體的最優值進行對比，并更新當前群體的全局最優適應值。

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

6)使用已經更新過的粒子速度和位置繼續計算目標函數，如果所得適應值滿足算法的精度要求，則算法結束，否則轉步驟2)繼續執行算法。

3.2 粒子群算法反演結果

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

表1 粒子群算法反演介質密度參數與真值

表2 粒子群算法反演參數及結果

表3 龍門山地區介質密度與Crust1.0對比

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

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

4 龍門山地區重力異常反演結果

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

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

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

圖5 剩余重力異常

圖6 龍門山地區密度分布

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

4.2 地下介質密度

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

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

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

表4 密度隨深度變化關系對比表

5 結論

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

致謝：

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