翟麗娜，賈曉東，焦明若，孔祥瑞，王松陽(yáng)

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

0 引言

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

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

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

1 基本原理

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

(1)

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

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

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

2 參考模型計(jì)算

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

根據(jù)上述對(duì)歐拉反褶積方法的認(rèn)識(shí)和應(yīng)用需求，結(jié)合遼寧地區(qū)重力觀測(cè)網(wǎng)絡(luò)形態(tài)和規(guī)模，設(shè)計(jì)符合研究區(qū)觀測(cè)規(guī)模的理論模型，該模型覆蓋上地殼、中地殼和下地殼不同深度位置場(chǎng)源的3個(gè)直面正六面體模型。該模型如圖1所示，在圖中矩形框?yàn)?個(gè)場(chǎng)源模型在地表觀測(cè)網(wǎng)上的水平投影，理論異常范圍在0～120×10-8m/s2，測(cè)網(wǎng)點(diǎn)間距設(shè)為5 km，測(cè)網(wǎng)長(zhǎ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個(gè)場(chǎng)源模型深度不同，具體參數(shù)信息如表1所示。

圖1 不同深度場(chǎng)源組合模型及正演異常

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

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

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

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

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

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

圖3 遼寧地區(qū)測(cè)點(diǎn)及斷裂分布圖

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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