基于歐拉反褶積方法計算遼寧地區重力變化場源特征*

翟麗娜，賈曉東，焦明若，孔祥瑞，王松陽

(遼寧省地震局，遼寧 沈陽 110034)

0 引言

地震是自然界最具破壞力的自然地質災害之一，是地球構造運動的一種表現形式。通常強震和大震是由于受到局部應力場以及周邊活動斷裂的影響而形成應力高度積累，從而在應力高度積累部位及其附近地區形成的明顯構造差異運動，往往還伴有重力場顯著變化(滕吉文等，2008；Chenetal，1979；Kuoetal，1999；Zhuetal，2010)。在地震地球物理觀測與分析的學科領域中，一般會采用定期流動重力觀測和連續重力重復觀測的測量技術，來捕捉整個設定測網覆蓋區內與場源變化相關的震前重力場變化信息(陳運泰等，1980；祝意青等，2009a；陳石等，2011)。地下介質變化過程在地表場兆反映中的重力場表現較為明顯，且其物理意義明確，在構造運動中是一種較好的表現形式(劉芳等，2016；高倩，陳石，2015)。祝意青等(2001，2003，2008，2009b，2013，2015)、陳石等(2011，2014)和Chen等(2016)將前人研究成果應用于我國川滇地區、新疆地區以及青藏高原東北部的中長期地震危險分析研究中，并捕捉到震前重力變化高梯度變化特征，為地震危險區的判定提供一定參考。

對于所獲取的地表觀測重力信號，如何開展其異常的定量化研究，一直以來都是研究人員的研究重點。焦新華和吳燕岡(2009)引入地球物理反演方法對重力場變化的場源特征進行研究，基于“以場溯源，場源結合”的分析思路對已知場源變化信息進行定量化解釋。在反演過程中需要精確分離地球物理變化過程中蘊含的復雜信息，來獲取有效觀測和精確的先驗信息以及參數，這一直以來也是重力場變化反演和解釋的難點(高倩，陳石，2015；劉芳等，2016)。

大地測量學的分析方法是地震重力場分析的常用方法，本文嘗試采用勘探地球物理學中的歐拉反褶積方法，對2012—2018年遼寧地區觀測的流動重力實測數據進行差分變化分析。結合研究區實際地質情況合理設計理論模擬的地表觀測重力場模型，對觀測到的流動重力場變化進行正演計算，選用歐拉反褶積方法反演場源參數，反復計算選取合適的構造指數和參數后，對2012—2018年實測數據進行反演和解釋。

1 基本原理

歐拉反褶積方法是由英國地球物理學家Reid等(1990)結合Thompson(1982)對歐拉齊次方程的研究總結后提出的方法，該方法可對重力位場數據進行快速反演，在先驗信息較少的情況下，在自動、半自動掃描后確定場源位置，找出構造范圍并推算其異常源信息(范美寧，2006)，計算公式為：

(1)

式中：(x，y，z)為觀測點位置；(x0，y0，z0)為場源位置；T為位場異常；N為構造指數；B為背景場。

構造指數N反映的是位場異常伴隨場源深度變化的衰減率，N參數的選取是場源位置確定的關鍵，一般與場源幾何結構有關。結合前人研究成果分析重力異常發現，N一般選擇0～2較為合適(范美寧，2006；劉芳等，2016；高倩，陳石，2015；魯寶亮等，2009)。

結合勘探地球物理研究中基于歐拉反褶積所建立的超定方程組，計算出最佳構造指數N，選取該構造指數進行歐拉反演計算成圖。具體步驟如下：①對離散數據進行網格化處理；②計算在x，y，z3個方向的梯度的位場異常；③根據異常體性質選擇合適的N和滑動窗口的大小且窗口滑動大小，最小為3×3以上；④解方程組求出異常場源參數等，并反復調整到最佳反演結果。

2 參考模型計算

歐拉反褶積方法一直以來在勘探學中應用較廣，它是在場源信息未知時就可對場源參數進行反演和解釋。計算要點主要是在反演時，最佳構造指數N與在x，y，z3個方向的導數組成線性方程組，計算求解即可。N在重力反演計算中與場源體幾何形狀相關，實際應用時要根據復雜的地質情況和場源形態多次試算出合適的N。再根據測點間的間距和異常形態特征等因素選取合理的滑動窗口，窗口大小能夠覆蓋場源異常，反復調整參數直至得到最佳反演結果(劉芳等，2016；高倩，陳石，2015)。

根據上述對歐拉反褶積方法的認識和應用需求，結合遼寧地區重力觀測網絡形態和規模，設計符合研究區觀測規模的理論模型，該模型覆蓋上地殼、中地殼和下地殼不同深度位置場源的3個直面正六面體模型。該模型如圖1所示，在圖中矩形框為3個場源模型在地表觀測網上的水平投影，理論異常范圍在0～120×10-8m/s2，測網點間距設為5 km，測網長度在X，Y兩方向上為180 km。其中M1模型位于中地殼，參數設為：X(50，60)，Y(100，110)，Z(10，13)，中心埋深11.5 km；M2模型位于下地殼，參數設為：X(80，110)，Y(40，70)，Z(20，25)，中心埋深22.5 km；M3模型位于位于上地殼，參數設為：X(80，110)，Y(150，155)，Z(5，6)，中心埋深5.5 km，密度變化選取為3.0×10-3g/cm3。3個場源模型深度不同，具體參數信息如表1所示。

圖1 不同深度場源組合模型及正演異常

表1 理論場源模型參數

通過上述設定模型計算分析后發現，當選取構造指數N=1、滑動窗口選5～8個測點距離時，反演誤差深度在15%以內，效果較好，并得到如圖2所示的最佳反演結果。圖中共計216組解，實際模型位置用矩形框表示，反演模型的結果用黑色圓點表示，計算結果聚集程度較好，將不同深度場源特征明顯區分出來。M2和M3模型在水平位置上誤差較小，結果較為理想，而M1誤差較大。由于上述模型試算結果較好，因此采用該參數對遼寧流動重力測網實測重力變化資料進行反演計算。

圖2 基于歐拉反褶積的反演結果三維空間分布

3 遼寧地區重力場變化的歐拉反褶積反演計算

遼寧地區處于太平洋板塊和歐亞大陸板塊交匯的構造復雜部位，屬華北地臺北部、赤峰—開原深大斷裂(內蒙地軸)以南的地臺區。著名的郯廬地震帶呈NW向、沿沈陽—海城—遼東灣—渤海東部—萊州灣—山東濰坊一線貫穿整個研究區。華北地臺、蒙吉—興安造山帶、郯廬斷裂帶基本控制著研究區的構造背景。遼寧地區的華北地臺部分自東向西依次分布有膠遼臺隆、下遼河斷陷、燕山臺褶帶。遼寧地區流動重力監測是我國東北部地區對郯廬地震帶中北段前兆監測的常規地球物理方法之一。

2012年以來，遼寧地區重力綜合測網每年觀測2期，觀測時間為每年5月和9月，承接101個重力測點、111個重力測段。根據學科發展的規劃，對遼寧地區地震重力區域測網進行調整、優化和改造，2018年以來的2期測量包括69個重力測點和84個重力觀測測段，測點分布見圖3。

圖3 遼寧地區測點及斷裂分布圖

利用2012年9月至2018年9月遼寧地區共計13期流動重力觀測數據進行分析，圖4給出了1年尺度重力變化差分結果。由于重力測點分布不均勻，筆者對其進行了空間差值。

(1)2012年9月至2013年9月重力分布圖顯示(圖4a)，遼寧地區整體以重力正變化分布為主，在遼寧南部岫巖—海城老震區和遼寧西部朝陽—北票地區出現明顯重力四象限現象，NW向和NE向有2個明顯的重力高梯度帶展布，遼寧東部地區尤其在沈陽—海城附近重力有自西向東由負轉正的趨勢變化。

(2)2013年9月至2014年9月的重力變化顯示(圖4b)，與上一年度相比全區發生正負重力異常反

向變化，遼寧西部地區朝陽一帶重力異常由正變負尤為明顯，該區朝陽—北漂斷裂、阜新—錦州斷裂梯度變化較大；遼東地區重力負值面積擴大。

(3)2014年9月至2015年9月重力場分布顯示(圖4c)，遼寧地區整體出現大范圍負變化，最大負值區出現在朝陽—北票斷裂附近，可達-40×10-8m/s2重力異常區；遼東地區由負轉正，發生正負變化較大說明斷裂斷層活躍有所增強。

(4)2015年9月至2016年9月，遼寧北部沈陽及其周邊重力場由負變正，遼西地區出現明顯重力四象限現象，變化最大處在朝陽—北票斷裂附近，達30×10-8m/s2，遼南海城—岫巖老震區亦出現重力四象限現象(圖4d)，且其周邊地震活動增加，表明老震區附近斷裂帶有所活動。

圖4 2012—2018年遼寧地區1年尺度重力變化

(5)2016年9月至2017年9月，與上一期相比重力負值區增大(圖4e)，總體重力變化較為平穩，說明斷裂活動穩定。

(6)2017年9月至2018年9月，在遼寧地區測網整體改造后，測點間距變大等值線稀疏，重力整體變化平穩，遼東出現重力高值變化區，測區重力變化等值線走向大體與上期一致(圖4f)。

本文選用歐拉反褶積方法對遼寧地區7期年尺度重力變化做反演，計算得到如圖5所示的重力變化的場源形態特征。圖中圓圈表示反演的場源位置，圓圈大小與場源深度一一對應，彩色實心圓圈直徑從小到大依次分別代表場源深度0～5 km，5～10 km，10～15 km，空心彩色圓圈代表場源深度15～25 km。選取紅、粉、黃、粉紅、藍、紫這6種顏色表示年尺度變化反演解，在反演計算時，以20 km為網格間距，以8倍間距為窗長，對2012—2018年共計7期年尺度流動重力數據計算，得到反演后位于地殼內部的場源深度為10～50 km。2012—2018年研究區內共發生5級以上地震1次，4級以上地震7次，3級以上地震48次。從反演場源位置分布整體特征上看，引起重力場變化的場源體位置落在6次差分結果上，其中重合度較好的在朝陽—北票斷裂和錦州—阜新斷裂帶區域、場源深度在20～40 km，和營口—海城斷裂(郯廬斷裂帶中北段)和莊河斷裂帶的條狀區域、其場源深度集中在20～50 km。

圖5 遼寧地區基于2012—2018年重力變化反演的場源特征空間分布

4 結論

通過歐拉反褶積方法的計算，在理論模型分析基礎上，選用最優解參數對遼寧地區流動重力數據重力場變化進行反演計算，得到與場源位置相關的幾何參數，獲得如下結論：

(1)歐拉反褶積方法在位場資料解釋中應用較廣，經過試算和分析發現該方法可反演重力場異常的場源三維空間位置，對多期重復觀測重力資料反演，發現持續聚集的位置區域與區域構造活動性增強區域一致。

(2)根據構建的理論模型和實驗結果，歐拉反褶積可以選擇合適的構造指數和滑動窗口大小，不需要給出準確場源的物性參數即可以得到較合理的反演結果。借助歐拉反褶積反演，得到的場源解位置集中的區域，往往更加接近真實解區域位置。通過計算遼寧地區7期年尺度重力場差分變化發現，重合度較好的地區集中在朝陽—北票斷裂附近NE向、海城—岫巖斷裂周圍NE向區域，場源深度集中在20～40 km，這個位置與近年來遼南及沿海地區和遼蒙交界地震活動頻繁活動地區相對應。

(3)遼寧臨海海域內并無實測點，在數據差值計算和反演時會產生虛假解，因此對于海域內反演結果深度及位置一致性較差。由于場源變化具有隨機性，進而部分地區集中度較低，雖然短期內存在一定的場源運動(物質運移)，但連續性較差。

綜上認為，采用歐拉反褶積方法對遼寧地區重力場變化進行分析，所識別的場源特征對于區域潛在危險性分析是有意義的。