伽師6.4級地震前后震源區視密度變化及其構造意義

劉代芹 陳 石 王曉強 張 貝 李 杰 吳傳勇 盧紅艷

1)中國科學技術大學地球和空間科學學院，合肥 230026 2)新疆帕米爾陸內俯沖國家野外科學觀測研究站，烏魯木齊 830011 3)中國地震局烏魯木齊中亞地震研究所，烏魯木齊 830011 4)新疆維吾爾自治區地震局，烏魯木齊 830011 5)中國地震局地球物理研究所，北京 100081 6)北京白家疃地球科學國家野外觀測科學研究站，北京 100095 7)中山大學地球科學與工程學院，廣東省地球動力作用與地質災害重點實驗室，廣州 510275

0 引言

通過研究重力場的靜態異常特征(布格重力異常、自由空氣重力異常、剩余均衡重力異常)和動態時間變化(隨時變化的內部物質交換、質量遷移等)，可獲取地殼內部結構、形態和特性的物理場變化信息(傅容珊，1988；滕吉文等，2006；Crossleyetal.，2013)。地震孕育和發生的整個過程是地殼內部、殼幔物質變動乃至核幔相互作用的結果，可使地殼出現物質運移、應力聚集和釋放等活動，這些運動特征將導致地震斷裂帶內部出現因質量遷移引起的不均衡變化，使得某些活動斷裂帶具備孕育中強地震的能力及背景。因此，認識地震的孕育及發生的全過程，離不開對地殼構造運動及物質遷移的系統性研究。地下介質結構變化的異常信息可以通過地表重力觀測進行捕捉。地球內部由于質量分布不均衡導致的狀態變化，常會引起地表觀測的重力異常發生變化。研究表明，地震活動區在震前較長的時間范圍內將會出現區域重力異常變化信息(祝意青等，2008；申重陽等，2009)。時變重力場異常是與地震孕育發生相關的前兆異常之一。以往的震例，如1997年3月伊豆半島地震群(Yoshidaetal.，1999)、2001年昆侖山口8.1級地震(祝意青等，2003)、2004年12月印尼9.1級大地震(Hanetal.，2006)、2008年汶川8.0級地震(Zhuetal.，2010；Zhangetal.，2020)、2013年蘆山7.0級地震(祝意青等，2013，2015)以及2015年尼泊爾MW7.8地震(Chenetal.，2016)等均表明，地殼內部的物質運移與地震孕育環境的變化有關。在強震發生前數年，地殼能量的匯集或驅散將在一定程度上導致地殼內部區域性的物質運移，特別是對于6級以上地震的孕震過程而言，在地表通?？捎^測到μGal量級的重力場變化過程。因此，通過開展高精度重力及地球物理探測工作研究地殼深部特征構造的物性結構、監測震前地球物理場異常變化動態，對深入認識地殼深部的孕震環境和識別前兆異常信息具有重要的科學研究價值和現實意義。

本文以 “以場求源、場源結合”的思想為指導，基于地表重復重力觀測數據，引入貝葉斯重力平差方法，將儀器的漂移率、氣壓導納、潮汐因子、儀器的格值系數等作為未知量，通過貝葉斯原理以及ABIC評價準則選取最優的參數值，在得到平差值的同時可以獲取每臺重力儀的漂移特性及格值系數等，進而得到時變重力點值序列；在此基礎上，采用時變重力信號的等效場源反演方法，最小化局部性淺源高頻干擾，并獲得研究區地殼內部近10a的區域性重力場源動態變化特征。具體而言，可分為以下幾個步驟： 首先，對塔里木盆地西緣的重力觀測網絡系統、等效源反演方法的原理和研究區的場源監測能力進行介紹；之后，對實際重復重力觀測數據進行了測試，并反演得到重力時空場源變化特征；最后，結合2020年伽師6.4級地震的發震構造區重力場變化過程進行分析與討論。

本文較全面地獲得了2020年1月19日伽師6.4級地震(39.83°N，77.21°E)前后柯坪斷裂及其周邊測網的重力場變化特征，為深入研究該區域深部的地殼物質運移過程和中強地震的孕育機理提供了地球物理依據。所得研究成果對于認識孕震構造背景和斷裂帶周邊構造運動模式具有重要意義。

1 塔里木盆地西緣的重力測網

2005年，新疆維吾爾自治區地震局在西南天山地區初步建立了由40個流動監測點構成的流動重力測網(喀伽網)，為了準確捕捉南天山地區的地震異常信息，分別于2012年、2013年對該重力網進行了優化改造，在原有的喀什-伽師監測網的基礎上，將監測范圍擴展到西昆侖、塔里木盆地、和田以及庫車一帶，形成了覆蓋阿克蘇、庫車、喀什、烏恰、塔什庫爾干、和田地區，共包括90余個重力監測點的陸地監測網(圖1)。

圖1 塔里木盆地西緣及周邊的地震重力觀測網絡Fig.1 Seismic gravity observation network in the western margin of Tarim Basin and its surrounding areas.白色四邊形為陸地重力觀測點，藍色五角星為絕對重力控制點，黃色圓圈為地震震中，紅色實線為活動斷裂構造

利用 2 臺高精度石英彈簧相對重力儀對監測網進行定期常規復測。測網內由庫車、塔什庫爾干和烏什3個絕對點提供時空重力基準約束(圖1 中的藍色五角星)。從塔里木盆地西緣的流動重力分布點位(圖1)可以發現，測點間距相對不均，測網西北部分布的測點較密。但塔里木盆地內部測點—和田地區一帶僅有1條測線通過，監測能力較弱，存在局部監控盲區。

據中國地震臺網測定，2020年1月19日伽師MS6.4地震的震源深度為16km，發震構造屬于柯坪塔格褶皺-逆斷裂帶的前緣斷裂(柯坪斷裂)(溫少妍等，2020)。該震區的主體位于柯坪逆沖推覆構造帶，由5、6排褶皺-逆斷裂帶及其之間的新生代盆地或谷地構成(馮先岳，1991；Burchfieletal.，1999；鄧起東等，2000；Fuetal.，2003；曲國勝等，2003；楊曉平等，2006；冉勇康等，2006；管樹巍等，2007；李安等，2016)。GPS數據研究表明南天山山前-盆地結合帶的變形明顯大于天山內部，在變形梯度最大的部位歷史上曾發生多次7～8級大地震(楊少敏等，2008；喬學軍等，2011；李杰等，2015；劉代芹等，2016)。天山地區的逆沖斷裂活動引發的地震破壞性極強(張培震等，1996；徐錫偉等，2006)，因此南天山山前盆-山結合部位的構造活動應該是我們需要重點關注的目標(陳杰等，2001；程建武等，2006；冉勇康等，2006)?？缕耗鏇_推覆構造系位于塔里木盆地的西北緣(馮先岳，1991；Burchfieletal.，1999；鄧起東等，2000；陳杰等，2006)，其北以喀拉鐵克斷裂與南天山晚古生代層系為界，南以柯坪塔格南緣逆斷裂與塔里木盆地內NWW向的塔里木西南坳陷、巴楚斷隆和阿瓦提坳陷為界，西界位于八盤水磨一帶，與阿圖什-八盤水磨逆沖構造帶相接，東端位于阿克蘇附近，與庫車逆沖構造系逐漸過渡(鄧起東等，2000；何文淵等，2002；李安，2013；吳傳勇等，2014)?？缕核衲鏇_推覆系自1990年以來共發生6級以上地震16次，2020年1月19日發生的伽師MS6.4地震是繼2003年巴楚-伽師6.8級地震后再次發生的6級以上地震，該地震打破了該地區長達近20a的6級地震平靜期(王筱榮，2001；高國英等，2007；溫少妍等，2020)。

另外，2015年7月3日在研究區重力監測的有效時間范圍內還發生了皮山MS6.5地震(37.6°N，78.2°E)，震源深度為10km。該地震位于青藏高原西北端和塔里木盆地的交界帶，震中距澤普斷裂較近，地震發生在推覆體前緣的皮山背斜上(李金等，2016；吳傳勇等，2017)。但該地震周邊的有效重力監測點較少，僅有1條測線通過(圖1)。

2 場源反演方法

基于陸地時變重力觀測研究地殼深部場源視密度結構特征時，通常需要考慮測量誤差和場源干擾2種不確定性。由于重力場隨時間變化的量級很小，通常在幾十μGal的水平，測量誤差主要體現在相對重力儀的非線性漂移和格值系數誤差2方面。本文基于貝葉斯平差方法(Chenetal.，2019)重新處理了2013年第2期(2013-C2)—2020年第1期(2020-C1)原始重力觀測資料(共14期)，獲得了重力點值變化結果。

此外，要應用這些點值變化結果研究地殼內部場源的視密度變化特征，還需要分離局部淺源干擾，常見的方法有低通濾波方法、小波分離方法、最小二乘配置、遺傳算法和神經網絡等優化插值類方法(Chackoetal.，1980；侯遵澤等，1997；Raoetal.，1999；樓海等，2005；柯小平等，2009；陳國雄等，2014；李長波等，2014；許闖，2014)，然而以上方法對于單期資料的處理效果較為明顯。

本文基于經典的重力位場反演方法理論(Lietal.，1998)，引入時變重力點值序列光滑的正則約束，以實現多期重力場變化的優化模型計算。采用 “等效源”近似描述時變重力場的場源視密度特征，這種等效類似于衛星重力數據反演的 “等效水厚度”(Schramaetal.，2007)概念。本文反演方法的優點在于，可以通過控制等效源的單元尺度壓制高頻干擾，通過場源的時空平滑正則約束項剔除局部性、脈沖式的場源干擾因素，進而更有效地獲得地殼深部的場源物性特征。

2.1 基本原理

對于時變重力場的等效源反演問題，目標函數可以采用式(1)(Lietal.，1998)描述：

Φ=W0Gm-g(x，y，t)2+Φs(m)+ΦT(m)

(1)

其中，等式右端第1項為觀測部分，第2和第3項分別為空間光滑約束和時間光滑約束，具體形式與模型的先驗假設有關。式(1)中，函數g表示不同期次的重力觀測時空變化點值；x、y、t分別為經度、緯度、時間；G為核函數，與等效源幾何參數和觀測點位置相關，可預先計算得出；m為待反演的等效源視密度參數；W0為權參數，與觀測數據的噪聲方差成反比。

下面將具體說明式(1)右端約束部分的構成。

2.1.1 空間約束

式(1)右端第2項Φs可以表示為

(2)

其中，W1、W2和W3是與模型在空間各方向光滑程度相關的權系數(超參數)，其具體的物理含義是期望待求的場源模型在X、XY和Y各個水平方向上都具有二階光滑特征，具體光滑程度與W1、W2和W3這3個待優化的未知超參數有關。

2.1.2 時間約束

式(1)右端第3項ΦT可以表示為

(3)

其中，W4是與模型在時間上的光滑程度相關的權系數(超參數)，其物理含義是期望待求的場源模型在時間上具有二階光滑特征，光滑程度由超參數W4控制。

在反演時，需要確定式(1)—(3)中的W0—W4共5個超參數后才可求解。為了合理地確定這些超參數，引入貝葉斯方法中的ABIC準則來實現超參數的優化計算(Akaike，1980)。具體而言，是通過ABIC最小化實現：

ABIC=-2log(L)+2N

(4)

式中，L是模型的似然函數，N是超參數個數。ABIC最小化問題是非線性問題，可以采用單純形方法或者牛頓法等具有多參數非線性優化能力的方法求解。

2.2 檢測板實驗

在反演實際重力場觀測數據之前，需要對不均勻分布的重力測網進行場源分辨能力測試。本文采用地震層析成像中常用的檢測板方法開展場源分辨率測試。由于單獨的重力異常不具有垂向分辨能力，因此，本文只進行水平分辨率的檢測板實驗。具體方法是設計一組正負相間的同尺度場源體，通過重力正演方法得到實際測點位置的理論重力異常，然后利用反演方法獲得場源模型參數，并與已知場源參數進行對比以評價由于實際測網分布不均勻導致的空間非均勻場源分辨能力的差異性。通常地震重復重力測量的空間覆蓋范圍可達幾百km，地球的曲率變化不能忽略，故采用球面六面體(Tesseroid單元)作為等效場源的離散化單元。

首先，采用0.5°×0.5°的場源體模型進行離散化，設定正負相間的場源體大小為1°×1°，每個場源體的視密度為±1×10-3g/cm3，場源體埋深為10km，每個場源模型的等效厚度為1km。通過理論計算得到每個等效場源體在地表可觀測到的最大重力異常范圍約為±32μGal。

圖2 是1°×1°分辨率的檢測板在地表可觀測的理論重力異常與地表實際測點的空間分布情況。由于地表實際的重力測點分布不均勻，對于已知的檢測板模型，實際僅能利用有限的陸地重力測網獲得觀測結果。本文在理論檢測板的重力異常值(圖2)中再加入5%的噪聲，采用上節所述的反演方法分別進行無噪聲和有噪聲干擾條件下的等效源視密度特征反演，得到的反演結果如圖3 所示。

圖2 研究區重力測網分辨率檢測板Fig.2 Gravity network resolution test board of the study area.

圖3 是基于圖2 所示的地表實際測網采集的重力異常反演后的場源參數結果。研究區測點間的平均距離為40～60km，因此，每個球坐標系下的六面等效場源體的大小設計為0.5°×0.5°。依據采樣定理，采樣頻率應大于信號最高頻率的2倍，因此，選擇1°×1°分辨率的檢測板進行實驗。本文采用的1°×1°實驗結果是相對合理的，圖3 的相關結果主要用于解釋不同區域的重力變化是否可靠提供參考。在圖3a中，距離測點較遠的區域無法恢復場源參數，在塔里木盆地西端的和田—阿克蘇之間也存在大面積的監測空白區。由于每個場源體模型的大小為0.5°×0.5°，而檢測板分辨率為1°×1°，故在例如和田西側的無測點覆蓋區內可以看到相鄰場源體之間的密度值存在一定的 “漸變”特征，這就是空間平滑正則化約束的結果。在圖3b的含噪聲反演結果中采用了與圖3a一致的顏色標注場源視密度特征，與圖3a進行對比可以看出整體上兩者區別不大，但在一些局部區域存在一定差異，主要表現為視密度值偏小。產生這種結果的原因是由于加入了隨機誤差，使得式(1)中與數據相關的估計權重參數W0減小。為了便于之后的分析，在圖3b中用黑色虛線標示出具備1°場源分辨能力的區域位置，這對于后續實際數據反演結果的分析和解釋具有指導意義。

另外在圖3b中可以看出，對于塔里木盆地南部的皮山MS6.5地震，震中周邊的場源檢測板模型無法恢復；而對于伽師MS6.4地震，地表重力測網對震中周邊的場源參數具有較好的分辨能力。因此，本文將重點分析北部的伽師MS6.4地震震中周邊區域的場源參數變化特征。雖然本研究資料的時間范圍也覆蓋2015年的皮山MS6.5地震的時段，但是由于測網對該地震震中及周邊監測能力較弱，故不再過多地嘗試提取該地震的震前異常特征。

圖3 測網反演的檢測板圖像Fig.3 Image of test board of the survey network obtained by inversion.a 無噪聲模型；b 含噪聲模型。白色四邊形為重力測點的位置；黑色實線為斷裂構造；黑色虛線為有場源監測能力范圍；黃色實心圓圈為地震震中

3 實際數據測試

在針對實際重力測網進行理論場源模型的分辨率測試基礎上，本文進一步對實際觀測到的重力數據進行平差和反演。陸地時變重力場一般每年觀測2次，分別在3—4月和8—9月進行，本文將上半年和下半年的觀測分別定義為C1和C2。

3.1 貝葉斯平差

本文共處理了2013-C2—2020-C1共14期測量數據，通過貝葉斯重力平差方法得到了每期的重力點值和誤差估計。首先，以2013-C2期的結果為基準，計算了全部14期的結果中具備同點位重復觀測重力點的累積變化圖像，選擇2015-C1—2020-C1共6期重力時間序列變化數據繪制了圖像，如圖4 所示。

圖4 貝葉斯平差后的重力場變化圖像(2015—2020年)Fig.4 Image of gravity field change after Bayesian adjustment(2015—2020).

在圖4a中，我們將圖3b所示的有效場源監測范圍用黑色虛線進行標記。圖4a—e為2020年伽師MS6.4地震前的重力場變化圖像，圖4f為震后的變化圖像。對比圖4a—e可以看出，在地震前的2017—2018年(圖4c，d)出現了較為明顯的重力場增加變化，且從2015—2018年的演化過程中可以看出，重力場的時空變化逐年增加，能量不斷積累，構造應力不斷增強，特別是在伽師6.4級地震之前，2017—2018年(圖4c，d)該地區重力值呈快速增加趨勢，表明震中附近的區域構造應力突然增強，物質可能不斷向該區域運移。但在2019年下半年，重力場開始趨于逐漸減弱的過程，在重力場正負交替階段即發生了伽師6.4級地震。地震發生之后3個月針對該重力網再次進行了復測，獲得了震后該地區的重力數據。從震后重力場變化圖像(圖4f)可以看出，重力場的時空演變特征繼續延續伽師6.4級地震之前(圖4e)的變化趨勢，即重力場在該網NW側出現大面積負值變化區域，特別是在震中出現了負值高值區，說明該地區在地震之后能量得到釋放，從而導致地殼物質遷移，構造應力處于逐漸減小或地殼密度不斷虧損的趨勢。從圖4 整體可以看出，2020年地震之后的重力場相對于2015年震前的重力場變化結果呈現區域性的減小態勢，這與2017—2018年震前區域重力場突增的趨勢特征明顯不同，這表明伽師6.4級地震發生之前地殼構造應力的積累較為集中，該特征與本次地震的孕育過程密切相關。

3.2 視密度反演

在上節的平差結果的基礎上，我們進一步應用第2部分的場源反演方法對全部期次的重力點值平差結果進行了模型反演。一般認為，重力變化與測點周邊的環境密切相關，如測點的高程變化、地下水位變化、土壤含水、河流、降雨、蒸發、徑流等。如果將這些不確定性統一歸類為場源因素，則這些場源因素非常復雜，若要定量地扣除每一種場源干擾，實施難度又極大。因此，本研究將采用等效源反演技術，通過在反演模型引入時空平滑先驗條件以壓制局部高頻干擾、周期性起伏波動變化等特征信號，進而實現深部時變重力場異常信號的提取和場源參數估計。

將與圖4 中每個期次對應的反演結果繪制于圖5 中，同樣在圖5a中用黑色虛線標注測網分辨能力較好的位置，且統一采用圖5a中的顏色標尺來可視化圖5b—f的結果。整體上，圖5 所示的重力場源視密度變化特征與圖4 的重力點值累積變化的過程相似，圖5 中的場源視密度變化范圍為±3kg/m3(地殼的平均密度約為2i700kg/m3)，反演得到的等效源層的密度變化約為正常地殼密度的1‰。從圖5 的結果可以看出，由于反演模型中的時空正則化約束，反演得到的場源模型時空演化過程更加平滑。

對比圖5d和5f在柯坪斷裂區域范圍內的場源視密度變化可以看出，地震前后出現了明顯的與構造走向NEE較一致的區域性場源視密度變化過程。那么，從反演結果與構造形態相似性的角度，在一定程度上可以認為本文得到的反演結果反映了與斷裂構造相關的深部重力場源信號。在地震前后(圖5d，f)，場源的視密度增加過程由正轉負，而且在地震后的視密度變化特征與柯坪構造高度相關，說明這些由重力監測獲得場源變化信號與此次地震事件及構造控制的場源環境變化相關。

對比圖5a—c可以看出等效源視密度變化由弱轉強。首先，在阿圖什及周邊區域出現近EW向的場源視密度變化特征，自圖5d開始變化方向與柯坪斷裂系統高度相關，在圖5d和5e中阿圖什—和田區域范圍內出現了1個小幅度的視密度減小區域，這個區域逐漸變大，發震后擴展至整個柯坪斷裂系范圍(圖5f)。

圖5 反演得到的地震前后震源區的視密度變化特征(2015—2020年)Fig.5 Variation characteristics of apparent density before and after earthquakes obtained by inversion(2015—2020).

綜上所述，MS6.4伽師地震前，場源視密度的變化過程具有繼承性發展特征，如果進一步從反演得到的等效源視密度變化的量級上估算，則本文研究區10km深度范圍內、1km厚度等效層內約1‰的密度變化可以被陸地時變微重力測量很好地觀測到。對于柯坪塊體而言，此次伽師地震前后的場源視密度演化特征，未來可作為典型震例輔助判斷強震風險源的位置。

3.3 殘差分析

按照前述的平差和反演過程，圖5 所示的場源視密度反演結果是基于圖4 的重力點值平差結果計算得到的，反演的意義在于可以對一些局部高頻的場源干擾和周期性的時變場源信號進行壓制。兩者之間的殘差，即通過視密度模型正演得到的重力點值與反演輸入的平差重力點值之差，就是反演模型剔除掉的局部場源干擾信號，該信號的直方圖特征如圖6 所示。

圖6 平差重力變化(a)、反演模型重力變化(b)和殘差重力變化直方圖(c)(2015—2020年)Fig.6 The histogram of adjustment gravity change(a)，gravity change obtained by model inversion(b)and residual gravity change(c)(2015—2020).

從圖6a中可以看出，全部14期參與反演的重力平差點值變化范圍為±120μGal，而反演視密度模型解釋的重力點值變化量約為±80μGal(圖6b)，其殘差部分的±40μGal重力變化是反演模型剔除掉的局部場源干擾和周期性高頻變化。對于本文研究的塔里木盆地西緣地區多年期的累積重力場變化信號而言，這個量級的干擾基本合理。

4 結論與討論

本文應用最新的時變重力貝葉斯平差方法和時變重力場視密度反演方法，系統性地跟蹤分析了2020年伽師6.4級地震前后的重力場變化特征，并反演獲得了場源視密度演化特征，得到的主要研究結論如下：

(1)塔里木盆地西緣及周邊地區現有的地震重力測網的場源監測能力對于柯坪斷裂系統及周邊范圍具有1°×1°的分辨能力，但對和田西側構造體系及至阿克蘇之間的塔里木盆地內部區域的監測能力相對較弱。

(2)伽師6.4地震前呈現的區域性視密度增加量級約為2kg/cm3，等效源模型的視密度變化約為正常地殼密度的0.74‰，該量級的場源多年累積的重力變化可以很好地被地表重力測網監測到。

(3)伽師6.4地震前的顯著重力變化自2017年開始，視密度變化整體上呈現區域性增加的趨勢，而在形態上先呈現EW向形態，后逐漸轉為NEE向，與柯坪斷裂系的構造方向趨于一致。2019年視密度變化趨勢減弱，發震后出現NEE向的視密度減小。

(4)伽師6.4地震后出現了與柯坪構造走向較一致的場源視密度減小特征，其主要視密度較小的位置位于地震震中—阿圖什一帶，這可能與震后斷層附近的構造應力得以迅速均衡調整，進而導致的地殼內部流體物質的再分布過程相關。另外，在伽師6.4級地震發生后3個月對塔里木盆地西緣及周邊進行重力復測，所得結果表明，由于與震后間隔的時間較短，震中附近區域場源的視密度變化仍延續震前減小的態勢，說明2020年4月觀測到的流動重力數據有可能包含著此次地震的同震效應信息。

綜上所述，地殼內部地質構造形變和地殼密度變化是地殼運動的2種基本形式，且兩者相互耦合。地殼介質因形變引起的密度變化產生了重力效應(陳運泰等，1979；申重陽等，2009；陳石等，2014)，地殼物質密度變化是地殼運動的基本形式之一，探測地殼物質質量的重新分布(或物質運移)引起的密度變化，對了解深部物質的分異、調整和運動以及構造分區深層過程的動力學機制具有重要意義。

本文圍繞伽師6.4地震前后的時空重力變化特征和等效場源視密度參數開展了一系列分析，所得結果可為新疆地區強震孕育、發生環境的研究提供地球物理場參考。此外，本文研究過程中采用的場源檢測板實驗方法可以為將來的區域重力測網設計、改進和優化布局提供指導。

致謝本研究基于所有參與天山地區流動重力觀測的同事們的辛勤付出才得以開展；中國地震局重力觀測技術管理部提供了相關基礎資料；GEOIST開源Python軟件包(1)https: ∥cea2020.gitee.io/geoistdoc為本文的模型測試和反演提供了支持；審稿專家對本文提出了建設性意見和建議。在此一并表示感謝!