劉代芹 陳 石 王曉強(qiáng) 張 貝 李 杰 吳傳勇 盧紅艷

1)中國(guó)科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院，合肥 230026 2)新疆帕米爾陸內(nèi)俯沖國(guó)家野外科學(xué)觀測(cè)研究站，烏魯木齊 830011 3)中國(guó)地震局烏魯木齊中亞地震研究所，烏魯木齊 830011 4)新疆維吾爾自治區(qū)地震局，烏魯木齊 830011 5)中國(guó)地震局地球物理研究所，北京 100081 6)北京白家疃地球科學(xué)國(guó)家野外觀測(cè)科學(xué)研究站，北京 100095 7)中山大學(xué)地球科學(xué)與工程學(xué)院，廣東省地球動(dòng)力作用與地質(zhì)災(zāi)害重點(diǎn)實(shí)驗(yàn)室，廣州 510275

0 引言

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

本文以 “以場(chǎng)求源、場(chǎng)源結(jié)合”的思想為指導(dǎo)，基于地表重復(fù)重力觀測(cè)數(shù)據(jù)，引入貝葉斯重力平差方法，將儀器的漂移率、氣壓導(dǎo)納、潮汐因子、儀器的格值系數(shù)等作為未知量，通過(guò)貝葉斯原理以及ABIC評(píng)價(jià)準(zhǔn)則選取最優(yōu)的參數(shù)值，在得到平差值的同時(shí)可以獲取每臺(tái)重力儀的漂移特性及格值系數(shù)等，進(jìn)而得到時(shí)變重力點(diǎn)值序列；在此基礎(chǔ)上，采用時(shí)變重力信號(hào)的等效場(chǎng)源反演方法，最小化局部性淺源高頻干擾，并獲得研究區(qū)地殼內(nèi)部近10a的區(qū)域性重力場(chǎng)源動(dòng)態(tài)變化特征。具體而言，可分為以下幾個(gè)步驟： 首先，對(duì)塔里木盆地西緣的重力觀測(cè)網(wǎng)絡(luò)系統(tǒng)、等效源反演方法的原理和研究區(qū)的場(chǎng)源監(jiān)測(cè)能力進(jìn)行介紹；之后，對(duì)實(shí)際重復(fù)重力觀測(cè)數(shù)據(jù)進(jìn)行了測(cè)試，并反演得到重力時(shí)空?qǐng)鲈醋兓卣?；最后，結(jié)合2020年伽師6.4級(jí)地震的發(fā)震構(gòu)造區(qū)重力場(chǎng)變化過(guò)程進(jìn)行分析與討論。

本文較全面地獲得了2020年1月19日伽師6.4級(jí)地震(39.83°N，77.21°E)前后柯坪斷裂及其周邊測(cè)網(wǎng)的重力場(chǎng)變化特征，為深入研究該區(qū)域深部的地殼物質(zhì)運(yùn)移過(guò)程和中強(qiáng)地震的孕育機(jī)理提供了地球物理依據(jù)。所得研究成果對(duì)于認(rèn)識(shí)孕震構(gòu)造背景和斷裂帶周邊構(gòu)造運(yùn)動(dòng)模式具有重要意義。

1 塔里木盆地西緣的重力測(cè)網(wǎng)

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

圖1 塔里木盆地西緣及周邊的地震重力觀測(cè)網(wǎng)絡(luò)Fig.1 Seismic gravity observation network in the western margin of Tarim Basin and its surrounding areas.白色四邊形為陸地重力觀測(cè)點(diǎn)，藍(lán)色五角星為絕對(duì)重力控制點(diǎn)，黃色圓圈為地震震中，紅色實(shí)線為活動(dòng)斷裂構(gòu)造

利用 2 臺(tái)高精度石英彈簧相對(duì)重力儀對(duì)監(jiān)測(cè)網(wǎng)進(jìn)行定期常規(guī)復(fù)測(cè)。測(cè)網(wǎng)內(nèi)由庫(kù)車、塔什庫(kù)爾干和烏什3個(gè)絕對(duì)點(diǎn)提供時(shí)空重力基準(zhǔn)約束(圖1 中的藍(lán)色五角星)。從塔里木盆地西緣的流動(dòng)重力分布點(diǎn)位(圖1)可以發(fā)現(xiàn)，測(cè)點(diǎn)間距相對(duì)不均，測(cè)網(wǎng)西北部分布的測(cè)點(diǎn)較密。但塔里木盆地內(nèi)部測(cè)點(diǎn)—和田地區(qū)一帶僅有1條測(cè)線通過(guò)，監(jiān)測(cè)能力較弱，存在局部監(jiān)控盲區(qū)。

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

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

2 場(chǎng)源反演方法

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

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

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

2.1 基本原理

對(duì)于時(shí)變重力場(chǎng)的等效源反演問(wèn)題，目標(biāo)函數(shù)可以采用式(1)(Lietal.，1998)描述：

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

(1)

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

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

2.1.1 空間約束

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

(2)

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

2.1.2 時(shí)間約束

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

(3)

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

在反演時(shí)，需要確定式(1)—(3)中的W0—W4共5個(gè)超參數(shù)后才可求解。為了合理地確定這些超參數(shù)，引入貝葉斯方法中的ABIC準(zhǔn)則來(lái)實(shí)現(xiàn)超參數(shù)的優(yōu)化計(jì)算(Akaike，1980)。具體而言，是通過(guò)ABIC最小化實(shí)現(xiàn)：

ABIC=-2log(L)+2N

(4)

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

2.2 檢測(cè)板實(shí)驗(yàn)

在反演實(shí)際重力場(chǎng)觀測(cè)數(shù)據(jù)之前，需要對(duì)不均勻分布的重力測(cè)網(wǎng)進(jìn)行場(chǎng)源分辨能力測(cè)試。本文采用地震層析成像中常用的檢測(cè)板方法開(kāi)展場(chǎng)源分辨率測(cè)試。由于單獨(dú)的重力異常不具有垂向分辨能力，因此，本文只進(jìn)行水平分辨率的檢測(cè)板實(shí)驗(yàn)。具體方法是設(shè)計(jì)一組正負(fù)相間的同尺度場(chǎng)源體，通過(guò)重力正演方法得到實(shí)際測(cè)點(diǎn)位置的理論重力異常，然后利用反演方法獲得場(chǎng)源模型參數(shù)，并與已知場(chǎng)源參數(shù)進(jìn)行對(duì)比以評(píng)價(jià)由于實(shí)際測(cè)網(wǎng)分布不均勻?qū)е碌目臻g非均勻場(chǎng)源分辨能力的差異性。通常地震重復(fù)重力測(cè)量的空間覆蓋范圍可達(dá)幾百km，地球的曲率變化不能忽略，故采用球面六面體(Tesseroid單元)作為等效場(chǎng)源的離散化單元。

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

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

圖2 研究區(qū)重力測(cè)網(wǎng)分辨率檢測(cè)板Fig.2 Gravity network resolution test board of the study area.

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

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

圖3 測(cè)網(wǎng)反演的檢測(cè)板圖像Fig.3 Image of test board of the survey network obtained by inversion.a 無(wú)噪聲模型；b 含噪聲模型。白色四邊形為重力測(cè)點(diǎn)的位置；黑色實(shí)線為斷裂構(gòu)造；黑色虛線為有場(chǎng)源監(jiān)測(cè)能力范圍；黃色實(shí)心圓圈為地震震中

3 實(shí)際數(shù)據(jù)測(cè)試

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

3.1 貝葉斯平差

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

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

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

3.2 視密度反演

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

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

對(duì)比圖5d和5f在柯坪斷裂區(qū)域范圍內(nèi)的場(chǎng)源視密度變化可以看出，地震前后出現(xiàn)了明顯的與構(gòu)造走向NEE較一致的區(qū)域性場(chǎng)源視密度變化過(guò)程。那么，從反演結(jié)果與構(gòu)造形態(tài)相似性的角度，在一定程度上可以認(rèn)為本文得到的反演結(jié)果反映了與斷裂構(gòu)造相關(guān)的深部重力場(chǎng)源信號(hào)。在地震前后(圖5d，f)，場(chǎng)源的視密度增加過(guò)程由正轉(zhuǎn)負(fù)，而且在地震后的視密度變化特征與柯坪構(gòu)造高度相關(guān)，說(shuō)明這些由重力監(jiān)測(cè)獲得場(chǎng)源變化信號(hào)與此次地震事件及構(gòu)造控制的場(chǎng)源環(huán)境變化相關(guān)。

對(duì)比圖5a—c可以看出等效源視密度變化由弱轉(zhuǎn)強(qiáng)。首先，在阿圖什及周邊區(qū)域出現(xiàn)近EW向的場(chǎng)源視密度變化特征，自圖5d開(kāi)始變化方向與柯坪斷裂系統(tǒng)高度相關(guān)，在圖5d和5e中阿圖什—和田區(qū)域范圍內(nèi)出現(xiàn)了1個(gè)小幅度的視密度減小區(qū)域，這個(gè)區(qū)域逐漸變大，發(fā)震后擴(kuò)展至整個(gè)柯坪斷裂系范圍(圖5f)。

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

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

3.3 殘差分析

按照前述的平差和反演過(guò)程，圖5 所示的場(chǎng)源視密度反演結(jié)果是基于圖4 的重力點(diǎn)值平差結(jié)果計(jì)算得到的，反演的意義在于可以對(duì)一些局部高頻的場(chǎng)源干擾和周期性的時(shí)變場(chǎng)源信號(hào)進(jìn)行壓制。兩者之間的殘差，即通過(guò)視密度模型正演得到的重力點(diǎn)值與反演輸入的平差重力點(diǎn)值之差，就是反演模型剔除掉的局部場(chǎng)源干擾信號(hào)，該信號(hào)的直方圖特征如圖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期參與反演的重力平差點(diǎn)值變化范圍為±120μGal，而反演視密度模型解釋的重力點(diǎn)值變化量約為±80μGal(圖6b)，其殘差部分的±40μGal重力變化是反演模型剔除掉的局部場(chǎng)源干擾和周期性高頻變化。對(duì)于本文研究的塔里木盆地西緣地區(qū)多年期的累積重力場(chǎng)變化信號(hào)而言，這個(gè)量級(jí)的干擾基本合理。

4 結(jié)論與討論

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

(1)塔里木盆地西緣及周邊地區(qū)現(xiàn)有的地震重力測(cè)網(wǎng)的場(chǎng)源監(jiān)測(cè)能力對(duì)于柯坪斷裂系統(tǒng)及周邊范圍具有1°×1°的分辨能力，但對(duì)和田西側(cè)構(gòu)造體系及至阿克蘇之間的塔里木盆地內(nèi)部區(qū)域的監(jiān)測(cè)能力相對(duì)較弱。

(2)伽師6.4地震前呈現(xiàn)的區(qū)域性視密度增加量級(jí)約為2kg/cm3，等效源模型的視密度變化約為正常地殼密度的0.74‰，該量級(jí)的場(chǎng)源多年累積的重力變化可以很好地被地表重力測(cè)網(wǎng)監(jiān)測(cè)到。

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

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

綜上所述，地殼內(nèi)部地質(zhì)構(gòu)造形變和地殼密度變化是地殼運(yùn)動(dòng)的2種基本形式，且兩者相互耦合。地殼介質(zhì)因形變引起的密度變化產(chǎn)生了重力效應(yīng)(陳運(yùn)泰等，1979；申重陽(yáng)等，2009；陳石等，2014)，地殼物質(zhì)密度變化是地殼運(yùn)動(dòng)的基本形式之一，探測(cè)地殼物質(zhì)質(zhì)量的重新分布(或物質(zhì)運(yùn)移)引起的密度變化，對(duì)了解深部物質(zhì)的分異、調(diào)整和運(yùn)動(dòng)以及構(gòu)造分區(qū)深層過(guò)程的動(dòng)力學(xué)機(jī)制具有重要意義。

本文圍繞伽師6.4地震前后的時(shí)空重力變化特征和等效場(chǎng)源視密度參數(shù)開(kāi)展了一系列分析，所得結(jié)果可為新疆地區(qū)強(qiáng)震孕育、發(fā)生環(huán)境的研究提供地球物理場(chǎng)參考。此外，本文研究過(guò)程中采用的場(chǎng)源檢測(cè)板實(shí)驗(yàn)方法可以為將來(lái)的區(qū)域重力測(cè)網(wǎng)設(shè)計(jì)、改進(jìn)和優(yōu)化布局提供指導(dǎo)。

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