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

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

0 引言

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

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

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

1 基本原理

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

(1)

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

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

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

2 參考模型計算

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

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

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

表1 理論場源模型參數(shù)

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

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

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

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

2012年以來，遼寧地區(qū)重力綜合測網(wǎng)每年觀測2期，觀測時間為每年5月和9月，承接101個重力測點、111個重力測段。根據(jù)學(xué)科發(fā)展的規(guī)劃，對遼寧地區(qū)地震重力區(qū)域測網(wǎng)進(jìn)行調(diào)整、優(yōu)化和改造，2018年以來的2期測量包括69個重力測點和84個重力觀測測段，測點分布見圖3。

圖3 遼寧地區(qū)測點及斷裂分布圖

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

通過歐拉反褶積方法的計算，在理論模型分析基礎(chǔ)上，選用最優(yōu)解參數(shù)對遼寧地區(qū)流動重力數(shù)據(jù)重力場變化進(jìn)行反演計算，得到與場源位置相關(guān)的幾何參數(shù)，獲得如下結(jié)論：

(1)歐拉反褶積方法在位場資料解釋中應(yīng)用較廣，經(jīng)過試算和分析發(fā)現(xiàn)該方法可反演重力場異常的場源三維空間位置，對多期重復(fù)觀測重力資料反演，發(fā)現(xiàn)持續(xù)聚集的位置區(qū)域與區(qū)域構(gòu)造活動性增強(qiáng)區(qū)域一致。

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

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

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