楊雄， 申重陽， 祝意青， 楊光亮， 孫凱， 談洪波， 王嘉沛

1 中國地震局地震研究所， 武漢 430071 2 中國地震局第二監(jiān)測中心， 西安 710054 3 太原理工大學(xué)， 太原 030024

0 引言

1976年7月28日河北唐山MS7.8級地震發(fā)生在燕山褶皺帶與華北裂陷盆地交接部位(圖1a)，亦即華北克拉通的東北部(圖1b)(朱日祥等，2012)，震中烈度達XI度，震源深度約11 km，對人民生命和財產(chǎn)造成巨大損害.地震發(fā)生后，為了探究唐山地震的成因機理，科學(xué)家們在唐山震區(qū)開展了大量的地震活動性、地質(zhì)構(gòu)造與深部地球物理研究(王椿鏞等，2017).

地震學(xué)結(jié)果顯示唐山MS7.8主震為近似直立的右旋走滑型破裂，破裂優(yōu)勢走向為北東30°，破裂方式為雙側(cè)不對稱破裂(Butler et al.，1979；張之立等，1980; 薛志照，1986).地表地質(zhì)調(diào)查表明，位于唐山塊陷中央的北東走向的唐山斷裂(圖1a)是這次MS7.8地震的發(fā)震構(gòu)造，唐山塊陷呈北東走向菱形狀，被北東向?qū)幒印钄嗔?、豐臺—野雞坨斷裂和北西向灤縣—樂亭斷裂、薊運河斷裂四大深斷裂所圍限(虢順民等，1977)；唐山斷裂具有高角度西傾的逆沖走滑性質(zhì)，晚更新世以來垂直錯距達15 m(尤惠川等，2002)；唐山斷裂是一條與褶皺相伴生的復(fù)雜斷裂構(gòu)造帶，主要由一系列相互平行的NNE-NE向的斷裂組成，自西向東分別為陡河斷裂、長山—巍山斷裂與古冶—唐山斷裂，這次大地震是晚第四紀曾發(fā)生多期位移事件的唐山斷裂帶的最新活動，地表破裂帶長約47 km，北段以右旋走滑為主，兼有西升東降垂直位移，南段以東升西降的傾滑活動為主，兼有右旋走滑活動(郭慧等，2011).

地震測深(深地震反射、深地震測深、密集地震臺陣等)獲取的殼幔速度結(jié)構(gòu)研究表明，唐山震區(qū)地殼可分為上、中和下三層結(jié)構(gòu)且殼內(nèi)界面波速差異較小(曾融生等，1988；段永紅等，2016；王椿鏞等，2017)，唐山菱形地塊下方上地幔頂部存在高度達10 km的異常隆起(劉啟元等，2007)，震源位于莫霍面隆起區(qū)的上方，震源下方的中下地殼存在明顯的低速體(劉昌銓和嘉世旭，1986；張先康等，1994；王未來等，2009；賴曉玲和孫譯，2013；段永紅等，2016)，認為低速體是深部巖漿沿裂隙侵入引起殼內(nèi)物質(zhì)熔融造成；唐山斷裂帶淺部構(gòu)造復(fù)雜表現(xiàn)為典型的花狀構(gòu)造，推測在中地殼合并為一條主斷裂延伸至莫霍面(劉保金等，2011)，與曾融生等(1988)推斷結(jié)果基本一致，但缺乏其他手段佐證.

大地電磁測深結(jié)果表明(劉國棟等，1983；秦馨菱等，1991；吳萍萍等，2021)，唐山震區(qū)電性呈三層結(jié)構(gòu)，淺表沉積層為低阻，上地殼為高阻層，下地殼為低阻層，唐山主震發(fā)生與高阻上地殼有關(guān).

重力探測比地震探測具有分辨率高、經(jīng)濟、覆蓋面廣等優(yōu)勢，被廣泛應(yīng)用于地球深部構(gòu)造探測等方面.方劍等(1999)利用華北地區(qū)的2°×2°網(wǎng)格的S波速度三維層析成像結(jié)果和15′×15′的網(wǎng)格平均布格重力異常，采用約束最小二乘方法反演得到了華北地區(qū)巖石圈內(nèi)6個層面上的密度分布結(jié)果；王新勝等(2012)利用華北地區(qū)5′×5′空氣重力異常資料和地震層析成像P波速度擾動資料，經(jīng)消除莫霍面和沉積層效應(yīng)，進一步反演給出了華北克拉通巖石圈三維密度結(jié)構(gòu)，反映出華北地區(qū)不同塊體間密度不均勻性，發(fā)現(xiàn)地震主要發(fā)生在殼內(nèi)低密度層附近.朱岳清等(1984)利用震前重復(fù)重力資料，探討了重力異常、莫霍面變形和唐山地震的孕育過程，提出了唐山地震地幔上隆孕震模式；Chen等(1979)研究了地震前后的重力變化，認為唐山地震與地殼和上地慢內(nèi)的質(zhì)量遷移有關(guān)；姜文亮和張景發(fā)(2012)利用布格重力異常反演了首都圈地區(qū)莫霍面，并對該區(qū)域的分層地殼結(jié)構(gòu)進行了深入分析；雷曉東等(2021)利用高精度區(qū)域重力資料識別研究了北京平原區(qū)淺部斷裂構(gòu)造.

這些研究成果大大加深了人們對唐山地震孕育發(fā)生背景的認識，但對唐山地震發(fā)震構(gòu)造及其深淺構(gòu)造關(guān)系缺乏足夠認識或更多證據(jù).已有的重力探測地下密度結(jié)構(gòu)成果較少，一般是較大區(qū)域(比如華北)或大尺度研究，缺乏震源區(qū)地下精細結(jié)構(gòu)的研究.此外，以往重力探測研究大多進行單一的密度體或密度界面反演研究，很少同時考慮兩者，這對重力資料解釋利用存在一定局限性.地震的孕育發(fā)生與地殼內(nèi)部的構(gòu)造運動有關(guān)，構(gòu)造運動伴隨著地殼形變及密度變化，從而引起重力場發(fā)生變化(申重陽，2005； Xing et al., 2021)，可以說重力場變化是地殼構(gòu)造運動的一種“化石”.為了較精細研究唐山震區(qū)地下密度結(jié)構(gòu)特征和進一步探討孕震發(fā)生機理，本文主要利用唐山地區(qū)較高分辨率的布格重力異常數(shù)據(jù)，經(jīng)有效分離和提取地下界面與局部密度體信息，反演研究震區(qū)地下密度三維結(jié)構(gòu)(密度體和密度界面)，并結(jié)合地質(zhì)與地球物理資料成果，從地下物質(zhì)遷移運動角度探討震區(qū)地震孕育發(fā)生的機理，為大震震源孕育的構(gòu)造條件識別和地震危險判定提供新的參考依據(jù).

1 布格重力異常數(shù)據(jù)

為了研究唐山地震的地下精細密度結(jié)構(gòu)和發(fā)震機制，中國地震局地震研究所在唐山地區(qū)開展了重力及GNSS(Global Navigation Satellite System)聯(lián)合測量，實測數(shù)據(jù)點共964個，點位分布如圖1a所示，點距約1～2 km、線距約5 km，重力觀測精度優(yōu)于0.01 mGal，GNSS觀測水平定位精度優(yōu)于0.04 m，垂直精度優(yōu)于0.06 m.圖2a為震區(qū)(39.3°N—40.7°N, 117°E—120°E)較精細的布格重力異常分布圖，其由實測數(shù)據(jù)及冀東地區(qū)1∶20萬布格重力異常圖(河北省地質(zhì)礦產(chǎn)局物探大隊1989年編制完成)經(jīng)統(tǒng)一、拼接而重新繪制.唐山震區(qū)的布格重力異?？傮w呈現(xiàn)出從北向南由負向正逐漸增加的變化趨勢，沿薊縣—遵化—撫寧—山海關(guān)一帶形成近東西向的重力梯度帶，北部的燕山隆起主要表現(xiàn)為負異常，南部的華北裂陷盆地以正異常為主，遠離梯度帶的兩側(cè)重力變化平緩，在興隆、青龍附近出現(xiàn)-40×10-5m·s-2左右的局部負異常，在唐山、豐潤、灤縣、山海關(guān)附近出現(xiàn)+30×10-5m·s-2左右的局部正異常區(qū)，唐山地震發(fā)生在重力局部高值區(qū)域，重力值高可能反映該區(qū)域地殼(或上地幔)存在莫霍面局部上隆(虢順民等，1977)或局部高密度體.

由于震區(qū)(圖1b虛線框)的空間范圍較小，后續(xù)研究需要利用較大區(qū)域(取36°N—42°N，112°E—121°E)重力場及其反映的深部構(gòu)造背景特征.為此，我們收集了國際公開的2′×2′全球重力場模型WGM2012(World Gravity Map 2012)，該模型數(shù)據(jù)源自EGM2008和DTU10，并從ETOPO1模型獲得的1′×1′分辨率地形校正，以球諧展開式計算得到的高精度地球重力場模型(Balmino et al., 2012).以震區(qū)地面重力數(shù)據(jù)為基準(圖2a)，將模型重力數(shù)據(jù)和震區(qū)地面重力數(shù)據(jù)進行擬合，去除系統(tǒng)差，再將模型重力數(shù)據(jù)和震區(qū)重力數(shù)據(jù)進行拼接處理，得到華北平原及鄰區(qū)布格重力異常(圖2b)，拼接后的震中區(qū)分辨率與實測數(shù)據(jù)一致，其他區(qū)域分辨率與WGM2012模型一致.從圖2b可以看出，從NW往SE區(qū)域重力異常變化總體呈現(xiàn)由負向正的變化趨勢，在中部存在一個NNE-SSW走向貫穿南北的重力梯度帶，即中國東部(大興安嶺—太行山—武陵山)重力梯度帶華北段(李安然等，1984)，重力梯度帶在石家莊附近發(fā)生偏轉(zhuǎn)，石家莊以南為NNE向，以北轉(zhuǎn)為NE向，重力梯度帶重力變化范圍為(-110～-20)×10-5m·s-2，北段梯度大，南段梯度小.重力梯度帶在北京附近出現(xiàn)分叉，一支繼續(xù)向北東方向延伸，另一支則轉(zhuǎn)向近東西向向東延伸(王椿鏞等，2017).太行山前重力梯度帶將華北平原及鄰區(qū)布格重力異常分為東西兩個異常變化區(qū)，東部重力變化平緩，變化范圍在(-20～50)×10-5m·s-2之間，在唐山震區(qū)附近存在局部圈閉，可能與唐山下方的莫霍面隆起或存在正密度體有關(guān)，在渤海附近表現(xiàn)為20×10-5m·s-2左右的局部正異常.西部主要為負異常變化，幅值在(-110～-180)×10-5m·s-2之間，在鄂爾多斯東緣達到最低值.

圖1 區(qū)域地形、活動斷裂與強震分布圖

圖2 布格重力異常

2 三維重力反演方法

重力異常是地下各層界面起伏和介質(zhì)密度不均勻的綜合反應(yīng)(楊文采等，2001；王新勝等，2012)，亦即重力異常是由許多具有一定密度的地下介質(zhì)離散體對測點單位質(zhì)量引力疊加所引起，這些地下介質(zhì)密度離散體可分為密度界面和局部密度體.密度界面一般指上下介質(zhì)明顯差異的近水平延伸界面(如莫霍面)，而局部密度體一般指內(nèi)外部密度差異明顯的局部圈閉體.設(shè)地表任一位置(x,y)觀測的重力異常Δg(x,y)可分解成

Δg(x,y)=Δgs(x,y)+Δgd(x,y),

(1)

其中Δgs(x,y)為密度界面效應(yīng)，Δgd(x,y)為密度體效應(yīng).重力反演的目的是最優(yōu)或合理確定這些密度界面與局部密度體分布，其關(guān)鍵問題是密度界面與局部密度體相關(guān)重力信號分離及其最優(yōu)反演.考慮到唐山地區(qū)地殼內(nèi)介質(zhì)不同界面之間波速差異較小、相應(yīng)的密度差異也較小(波速與密度一般成正比，圖3密度由波速-密度經(jīng)驗轉(zhuǎn)換關(guān)系計算)，后述研究只考慮莫霍面.基本研究思路為：利用布格重力異常數(shù)據(jù)，提取與莫霍面相關(guān)的區(qū)域趨勢異常(簡稱區(qū)域異常)；利用區(qū)域異常反演確定莫霍面起伏；利用扣除莫霍界面效應(yīng)的剩余布格異常，反演研究地下介質(zhì)三維密度體分布.

圖3 區(qū)域地殼介質(zhì)參數(shù)一維模型；實線為P波速度(黃雪源等，2021)，虛線為轉(zhuǎn)換密度(馮銳等，1986)

2.1 莫霍界面反演算法

已有研究表明，一定區(qū)域條件下，小波多尺度分解信號可較好反映不同構(gòu)造深度的密度體信號(侯遵澤等，2014； Xuan et al.，2016；雷曉東等，2021)，其中小波逼近信號與莫霍面起伏引起的重力異常具有較好一致性(侯遵澤和楊文采，1997；楊文采等，2001)，可用來反演莫霍面.因此，本文反演的基本思路是，先選取較大區(qū)域重力異常資料，利用小波多尺度分解提取與莫霍面相關(guān)的區(qū)域重力異常，在此基礎(chǔ)上采用空間域直接迭代法反演莫霍面.

小波多尺度分析通過尺度因子的變化可將信號分解為小波細節(jié)與逼近，其分別對應(yīng)高頻與低頻部分，對各尺度異常細節(jié)和逼近進行疊加重構(gòu)可得到局部異常和區(qū)域異常，實現(xiàn)重力場的二分，具體可表示為

(2)

其中Δg(x,y)為地表的觀測異常，AN(x,y)為N階小波逼近，Dj(x,y)為各尺度的小波細節(jié)，N為分解的最大階，階數(shù)N是人為選定的，小波變換的低階小波細節(jié)不變準則說明不管N如何選擇，小波變換出來的低階小波細節(jié)都是一樣的，不同的是小波細節(jié)的個數(shù)和N階逼近.

密度界面反演根據(jù)計算域的不同可分為空間域和頻率域，空間域方法主要包括直接迭代法，脊回歸法、壓縮質(zhì)面法和正則化方法等(馮娟等，2014；馮旭亮，2019)，頻率域主要是Parker-Oldenburg迭代方法(Parker，1972；Oldenburg，1974；馮娟等，2014).空間域方法中直接迭代法應(yīng)用較多，Bott(1960)首先提出用于沉積盆地基底反演，后經(jīng)不斷完善發(fā)展得到廣泛應(yīng)用，如：Silva等(2014)對迭代公式進行改進，提高了收斂速度；Santos等(2015)又引入正則化原理用于三維非光滑盆地基底反演.直接迭代算法基本原理是利用無限大平板引起的重力異常逐次逼近相關(guān)重力異常，可表示為

i=1,…,N,

(3)

其中，(xiyi)為第i個重力觀測點的坐標；mk-1(xiyi)和mk(xiyi)分別為(xiyi)點處第k-1次迭代和第k次迭代的密度界面深度反演結(jié)果；g0(xiyi)和g(xiyi)k-1分別為(xiyi)點處實測重力異常和第k-1次迭代反演模型正演計算的重力異常；Δρ為界面上下密度差.該方法用來確定地下界面的主要特點是原理簡單、計算速度快，只要較好分離提取與界面有關(guān)的重力異常信號，可達到理想效果.

2.2 局部密度體反演

Camacho等(2000)年提出的基于塊體生長的三維密度反演方法，在火山地區(qū)的應(yīng)用研究中取得較好效果(Camacho et al.，2000，2011).該算法先將研究區(qū)地下介質(zhì)體離散成一系列規(guī)則形狀(如棱柱體)單元，反演時模擬細胞生長的方式，即以規(guī)定的密度差Δρ填充某個單元，搜索使模型擬合度和平滑度達到最小的塊體單元，聚合生長這些單元以最優(yōu)擬合觀測到的重力異常值.該算法的優(yōu)點是需要的先驗信息少，且以自動客觀的模型空間探索方式構(gòu)建模型(Camacho et al.，2000；Berrino and Camacho，2008).

(4)

其中Aij表示單位密度的第j個棱柱體單元在第i個觀測點引起的重力異常(Nagy，1966；Camacho et al.，2000).

+py(yi-yM))=vii=1,…,n

(5)

(6)

塊體生長或反演過程中，目標函數(shù)采用基于模型擬合度和模型平滑度的混合最小化條件：

(7)

其中，v=(v1,…,vn)T是第j個模型殘差的列向量，m=(Δρl1,…,Δρlk,Δρj)T是第j個模型密度差的列向量，λ為平衡模型擬合度和模型平滑度的因子，QD是重力數(shù)據(jù)不確定性的先驗協(xié)方差矩陣，QM是模型參數(shù)不確定性的協(xié)方差矩陣.

生長過程中，對于每個可能聚合生長的新單元(以規(guī)定的正或負密度差填充)，比例因子f和趨勢參數(shù)p0,px,py在求解(5)式中同時被調(diào)整，再通過式(7)計算整體模型的殘差值.最后，選擇使模型殘差最小的第j個棱柱體進行聚合生長.對于每個連續(xù)步驟，調(diào)整后的比例值f減小，趨勢參數(shù)p0,px,py幾乎達到穩(wěn)定值(這里應(yīng)為零).當f接近1時，生長過程停止，生成完整的三維異常體(Camacho et al.，2000，2011).

3 反演計算與結(jié)果分析

3.1 華北及鄰區(qū)莫霍面的反演

由于震區(qū)(圖1a和圖2a)范圍較小，難以得到較好的莫霍面埋深信息，為此將研究區(qū)擴展到華北及鄰區(qū).

3.1.1 區(qū)域重力異常信號計算與分析

為了分離出由深部莫霍面起伏引起的重力異常，利用小波變換將華北及鄰區(qū)布格重力異常進行多尺度分解.不同尺度的小波細節(jié)和逼近對應(yīng)不同的場源深度(楊文采等，2001)，利用功率譜分析計算了華北平原及鄰區(qū)各階小波細節(jié)和小波逼近的近似場源深度，得到華北平原及鄰區(qū)1階、2階、3階和4階小波細節(jié)對應(yīng)的近似場源深度分別為3、11.5、18.5 km和28 km，4階小波逼近的對應(yīng)的近似場源深度為41 km.根據(jù)前人的研究結(jié)果，華北平原及鄰區(qū)的莫霍面深度一般在31～43 km(熊小松等，2011；段永紅等，2016)，這樣小波分解階數(shù)選取為4階即可.對圖2b的布格重力異常進行二維離散小波分解，取4階逼近為區(qū)域異常，得到華北平原及鄰區(qū)區(qū)域異常(圖4).由圖4可看出，區(qū)域異?？傮w呈現(xiàn)NW向SE由負向正逐漸增加的趨勢，變幅在(-170～40)×10-5m·s-2；在中部存在一條貫穿南北的北北東向重力梯度帶(即中國東部布格重力異常梯級帶)，梯度帶變化范圍(-120～-10)×10-5m·s-2，重力梯度帶在石家莊以南為NS向，石家莊-延慶走向NE，延慶以北區(qū)域又轉(zhuǎn)向NNE，遠離梯度帶的區(qū)域重力變化平緩.梯度帶以東主要呈現(xiàn)正異常區(qū)，變化范圍在(-10～40)×10-5m·s-1，梯度帶以西主要以負異常為主，變化范圍在(-120～-170)×10-5m·s-2.

圖4 華北平原及鄰區(qū)區(qū)域異常

3.1.2 莫霍面反演結(jié)果與分析

根據(jù)前人的研究結(jié)果(熊小松等，2011；段永紅等，2016；姜磊等，2021)，選取莫霍面的平均參考面深度為37 km、上下密度差為-0.4 g·cm-3.利用圖4的區(qū)域異常和上述界面反演方法得到華北平原及鄰區(qū)的莫霍面(圖5a).由圖5a可以看出，華北平原及鄰區(qū)的地殼厚度為31～45 km，總體表現(xiàn)為東薄西厚，盆地薄、山區(qū)厚的特征，東南部的華北盆地地殼厚度為31～35 km，唐山震區(qū)附近位于華北中部莫霍面陡變帶與東部莫霍面隆起區(qū)的過渡部位；南部的魯西隆起地殼厚度為34～35 km，從華北盆地到太行山隆起，地殼厚度從35 km急劇增加到41 km；北部的燕山隆起區(qū)地殼厚度由東邊渤海的34 km向西逐漸增加至張家口西部的43 km，地殼厚度最深的地方在西北部達到了44.5 km.莫霍面等值線的走向主要為北東向，在太行山隆起、燕山隆起與華北盆地交匯處莫霍面梯度最大，與重力異常梯度帶相對應(yīng).重力梯度帶和莫霍面陡變帶反映了太行山、燕山與華北盆地交匯區(qū)構(gòu)造運動強烈，太行山是華北克拉通內(nèi)部非常重要的構(gòu)造轉(zhuǎn)換帶(段永紅等，2016).

圖5b為華北平原及鄰區(qū)莫霍面(圖5a)與HBCrust 1.0(段永紅等，2016)的比較.總體來看，華北平原及鄰區(qū)莫霍面與HBCrust1.0結(jié)果基本一致，大部分區(qū)域偏差在2 km左右；唐山震區(qū)莫霍面偏差較小(小于2 km)；較大偏差在華北盆地的西南側(cè)，誤差達到4 km左右，可能與該區(qū)域地震測深剖面較少有關(guān).將本文反演結(jié)果和前人研究結(jié)果進行統(tǒng)計比較(表1)，可以看出研究區(qū)域莫霍面深度總體上具有一致性，但仍然存在一些差異，這可能與所采用觀測資料(如臺站數(shù)、剖面數(shù)及其分布)、反演方法的局限性或?qū)嶋H波速與密度之間存在的非線性有關(guān).本文利用的重力數(shù)據(jù)具有較好的均勻空間分辨，反演結(jié)果具有較好的可信度.

表1 華北平原及鄰區(qū)莫霍面與前人研究結(jié)果統(tǒng)計值

3.2 唐山震區(qū)密度體反演結(jié)果

利用公式(1)，并考慮到4階小波逼近得到的重力異?？赡苓€包含一些不是莫霍面(如密度體)引起的信號，故從震區(qū)布格重力異常(圖2a)中扣除震區(qū)莫霍面(圖5a)的重力效應(yīng)，得到震區(qū)剩余重力異常(圖6)，這種做法的優(yōu)點是密度反演結(jié)果的物理意義更加清晰些.剩余重力異常主要反映地殼內(nèi)部不均勻密度效應(yīng)，震區(qū)南部剩余重力異常以正異常為主，正異常主要分布在唐山菱形塊體及渤海灣附近，其中寶坻、豐潤及山海關(guān)附近存在多個局部重力高，武清附近存在局部重力低；震區(qū)北部以負異常為主，其中興隆、青龍和山海關(guān)北存在3個局部重力負異常區(qū).

圖5 華北平原及鄰區(qū)莫霍面

以震區(qū)剩余重力異常(圖6)為地表約束，將震區(qū)地下結(jié)構(gòu)根據(jù)重力點的空間分辨能力劃分為平均邊長5000 m的35000個棱柱體單元，為了避免邊界效應(yīng)，將震區(qū)范圍沿水平方向向四周延拓114 km.采用唐山地區(qū)的地殼分層速度結(jié)構(gòu)以及速度擾動信息(滕吉文等，1979；楊峰，2019)，利用馮銳等(1986)在華北地區(qū)的速度與密度轉(zhuǎn)換關(guān)系，可得到密度差取值范圍為(-160,160)kg·m-3.過高的λ值反演得到的模型簡單，擬合度差，過低的λ值反演得到一個非常復(fù)雜的模型且擬合度高，通過不斷的試錯得到一個λ=11.3的最優(yōu)值.由于不需要調(diào)整區(qū)域趨勢，所以在反演中將區(qū)域趨勢的系數(shù)全部設(shè)為0，利用得到的剩余重力異常、密度差及λ值，最終反演出地下密度結(jié)構(gòu)分布圖(圖7和圖8)，其中圖7為震區(qū)地下密度結(jié)構(gòu)水平切面圖，圖8為震區(qū)地下密度結(jié)構(gòu)垂直切面圖.

圖6 震區(qū)剩余重力異常；圖中黑線表示垂向剖面位置

圖7a為5 km深度的密度反演結(jié)果，總體呈現(xiàn)為正負相間分布的小局部異常體，但仍可看出大致以平谷—遷西—山海關(guān)為界，基本呈現(xiàn)南正北負的格局.北部以負密度異常體為主，在興隆、青龍及山海關(guān)北部出現(xiàn)多個局部負密度異常體，與燕山隆起區(qū)的興隆凹陷等多個局部凹陷區(qū)一致；南部以正密度異常體為主，在寶坻斷裂附近呈現(xiàn)顯著分異的北正南負局部密度異常體，分別對應(yīng)寶坻凸起和武清凹陷；潘莊西斷裂的東南，豐臺附近也出現(xiàn)北東走向的正密度異常體，對應(yīng)豐臺凸起；樂亭附近的負密度異常體和灤縣附近的正密度異常體分布在寧河—昌黎斷裂兩側(cè)，分別對應(yīng)樂亭凹陷和灤縣凸起，反映該斷裂的切割作用.

圖7b為10 km深度(與唐山地震的震源深度相當)密度反演結(jié)果，震區(qū)正負密度異常體位置與5 km深度幾乎一致，北負南正的格局更加明顯；北部負密度體更為突出，南部除東南北戴河局部正密度體外，正密度體呈現(xiàn)倒R型分布：唐山斷裂附近出現(xiàn)明顯的正密度異常體，且與西南豐臺凸起的正密度異常體相連；寶坻—玉田、遷西、遷安和灤縣附近均出現(xiàn)明顯的正密度異常體..

圖7 地下密度結(jié)構(gòu)水平剖面圖；現(xiàn)代地震(1970—2021)

圖7c為15 km深度密度反演結(jié)果，與10 km結(jié)果基本一致，繼續(xù)維持北負南正的格局，南部正密度體仍然呈倒R型.隨著深度逐漸增加，密度體開始聚攏，反映淺部構(gòu)造的一些小局部異常體消失，例如平谷、薊縣、撫寧、昌黎等密度體；震區(qū)南部沿唐山斷裂分布的正密度異常體匯聚到震中附近，遷西的正密度異常體向南(震中方向)偏移，寶坻的正密度體收縮，倒R型稍微減弱；遷安、樂亭附近正負密度異常體幾乎消失，山海關(guān)隆起的正密度異常體更加顯著.

圖7d為20 km深度密度反演結(jié)果，繼續(xù)維持北負南正的格局，南部正密度體倒R型分布近乎消失.北部的負密度異常體收縮，遷西附近的正密度異常體增大，與震中的正密度異常體相連，寶坻和灤縣附近的正密度異常體消失，山海關(guān)的正密度異常體變小，武清附近負密度異常體收縮.

圖7e為25 km深度密度反演結(jié)果，北負南正的格局明顯減弱.北部負密度異?；鞠?；南部正密度體倒R型分布消失，收縮至遷西—豐南一帶，遷西的正密度異常體向震中方向遷移、收縮，反映了震中區(qū)域深部構(gòu)造的活躍性.圖7f為30 km 深度密度反演結(jié)果，北負南正的格局消失，震區(qū)內(nèi)負密度異常體消失，正密度異常體繼續(xù)向震中附近(豐潤—唐山)匯聚，在豐潤附近呈現(xiàn)正密度異常體.

為了研究唐山斷裂與密度體的關(guān)系，分別沿平行和垂直唐山斷裂(參見圖6)切取垂向密度結(jié)構(gòu)剖面圖.平行唐山斷裂剖面圖(圖8a)顯示唐山震中下方正密度體從7 km左右向下延伸至30 km左右與豐潤、遷西下方的正密度重合，繼續(xù)向下延伸至50 km左右，形似Y型，反映唐山下方的正密度體與遷西、豐潤下方的正密度體屬于同一構(gòu)造; 垂直唐山斷裂剖面圖(圖8b)顯示密度體呈棍柱狀垂直向下切穿莫霍面，密度異常體在地殼內(nèi)部的走向與唐山斷裂及推斷的地殼深斷裂(劉保金等，2011)基本一致；唐山斷裂下方的正密度體形態(tài)明顯反映了唐山斷裂活動時的物質(zhì)運移狀態(tài).

圖8 地下密度結(jié)構(gòu)垂直剖面圖，黃色虛線表示莫霍面深度，(b)中的黑線為深地震反射剖面結(jié)果(劉保金等，2011)

4 討論與結(jié)論

4.1 孕震構(gòu)造環(huán)境條件分析

從大區(qū)域來看，唐山地震所處的華北克拉通及鄰區(qū)在太平洋板塊西向俯沖、西伯利亞板塊南南東擠壓和印度板塊北東向碰撞或匯聚作用下，中新生代以來主要經(jīng)歷了燕山運動和克拉通破壞，巖石圈發(fā)生伸展減薄并伴隨著大量的巖漿活動(朱日祥等，2012；王瑜等，2018；董樹文等，2019).華北平原及鄰區(qū)布格重力異常(圖2b)和莫霍面起伏(圖5a)大體反映了區(qū)域地質(zhì)構(gòu)造長期活動引起的物質(zhì)質(zhì)量分布基本構(gòu)架：從西往東，重力異常逐步增大，莫霍面逐步變淺，可分為西部、中部和東部三個小區(qū).西部包括鄂爾多斯東緣和燕山隆起西緣，重力異常較低(小于-120×10-5m·s-2)、莫霍面較深(大于41 km)，主要受到青藏地塊東向推擠和西太平洋板塊南西西向反向擠壓作用；中部包括山西地塹、太行山隆起和燕山隆起，重力異常呈現(xiàn)寬緩的重力梯度帶(介于-120×10-5～-10×10-5m·s-2之間)、莫霍面呈現(xiàn)往東抬升的陡坡(介于35～41 km之間)，中部重力梯度帶兩側(cè)地殼厚度差異反映了華北地區(qū)東西兩側(cè)巖石圈經(jīng)歷了不同程度的改造過程，為華北克拉通內(nèi)部非常重要的構(gòu)造轉(zhuǎn)換帶，其北部燕山隆起附近可能受到西伯利亞板塊南南東擠壓遠程作用和華北平原板塊的向北伸展作用，也曾出現(xiàn)大量巖漿活動(王瑜等，2018)；東部包括華北平原，重力異常較高(大于-10×10-5m·s-2)、莫霍面較淺(小于35 km)，可能與華北裂陷盆地內(nèi)巖漿上涌作用有關(guān).這種構(gòu)架的形成主要與西部青藏地塊東向推擠、東部太平洋板塊西向擠壓、北部西伯利亞板塊和南部華南地塊阻擋長期聯(lián)合作用有關(guān)，也與華北克拉通破壞密切相關(guān).華北地區(qū)地震構(gòu)造應(yīng)力場的平均主壓應(yīng)力場(P軸)為NEE-SWW向，反映深部各向異性的SKS波優(yōu)勢方向為NWW-NW向(王椿鏞等，2017)，與研究區(qū)內(nèi)莫霍面等值線的NE優(yōu)勢走向基本一致，說明板塊之間的水平作用起到十分重要作用.在布格重力異常與莫霍面起伏線性趨勢明顯的地帶往往具有較強的地震活動，如：中部梯度帶南段的西側(cè)和東側(cè)分別存在與梯度帶走向基本平行的汾渭地震帶和華北平原地震帶(圖1b)；中部梯度帶拐彎部位被張渤地震帶垂直穿過，唐山地震就發(fā)生在重力異常梯度帶及莫霍面陡變帶東側(cè)的隆起區(qū)，反映震區(qū)上地幔存在垂直向上的推力.

從唐山震區(qū)來看，震區(qū)的布格重力異常(圖2a)南高北低特征展示了燕山隆起與華北平原過渡區(qū)之間的差異性，震中附近顯著正異常可能反映地下高質(zhì)量物質(zhì)活動，震區(qū)地殼密度結(jié)構(gòu)(圖7和圖8)給予了進一步印證.密度體結(jié)構(gòu)圖像(圖7和圖8)的顯著特點有二：一是南正北負，南部與北部之間介質(zhì)密度的這種分異性一直延深至25 km(上、中地殼)；二是鄰近唐山斷裂的顯著高密度體三維特征與斷裂帶產(chǎn)狀相吻合，從地下7 km一直延伸至40 km(上地幔).從整體來看，地殼密度的南正北負分布反映了華北盆地與燕山隆起在中國東部構(gòu)造應(yīng)力場作用下的差異構(gòu)造運動，震區(qū)北部的低密度體可能反映了燕山隆起過程中，由于重力均衡作用調(diào)整使地殼處于物質(zhì)虧損狀態(tài)；震區(qū)南部地殼高密度體特征則反映出華北伸展盆地受到高質(zhì)量物質(zhì)的大量流入或侵入(伸展作用下地殼密度應(yīng)該降低)，從高密度體均呈現(xiàn)小局部特征說明高質(zhì)量物質(zhì)侵入的可能性更大.從震區(qū)密度結(jié)構(gòu)水平截面圖(圖7)可看出，淺部分散高密度體隨深度增加逐步聚合特征：10、15 km深局部高密度體呈倒R分布(寶坻—遷西—遷安—灤縣—唐山)；20、25 km深局部高密度體主要分布在遷西—唐山，至30 km收縮至豐潤—唐山，這種特征似反映了高密度體源于更深處(上地幔)巖漿的上侵流動，在地殼內(nèi)部形成相互貫通的巖漿囊體(高密度)，這與燕山運動、華北克拉通破壞過程中大量巖漿活動有關(guān).從震區(qū)密度結(jié)構(gòu)垂直截面圖(圖8)可看出，鄰近唐山斷裂的顯著高密度體貫穿莫霍面至上地幔，平行唐山斷裂走向(圖8a)高密度體在莫霍面以上明顯發(fā)散，而垂直唐山斷裂(圖8b)高密度體一直呈現(xiàn)條柱狀，且密度體走向與劉保金等(2011)唐山斷裂帶下方存在錯斷中下地殼的深斷裂一致，這更進一步說明地幔巖漿易沿斷裂薄弱部位上侵活動的特點.唐山斷裂處深部高密度體特征說明：唐山斷裂切穿莫霍面至上地幔，上地幔巖漿沿斷裂向上活動一方面造成莫霍面局部隆起，同時巖漿沿先存唐山斷裂面向上侵入，在地殼中形成巖漿囊，直至地下7 km左右受到地表沉積成巖作用而受到阻擋，熱的巖漿充填了唐山斷裂深部破裂面.從地震和大地電磁測深展示的中地殼和下地殼存在P波、S波低速體和高導(dǎo)層特征來看，這些侵入的巖漿仍然處于高溫狀態(tài)，這可能為唐山地震提供了能量積累的重要源泉.

4.2 發(fā)震力學(xué)模式分析

從已有研究成果和上述分析來看，巖漿沿斷裂面的侵入作用對唐山地震的形成和發(fā)展具有重要作用(朱岳清，1984；李紹柄，1986).一般認為唐山地震孕育發(fā)生的力源主要包括兩個方面：水平向力源和垂直向力源.華北及震區(qū)重力異常、莫霍面、以及密度體特征優(yōu)勢方向與華北區(qū)域構(gòu)造應(yīng)力場方向和SKS波各向異常結(jié)果大體一致，說明唐山地震震源處水平作用力的東西分量主要源于太平洋板塊的俯沖和青藏塊體東向流動，水平力的南北分量主要源于華北平原往北伸展作用和西伯利亞板塊的阻擋，進而引起的燕山隆升均衡作用(密度的南高北低).從震區(qū)密度體結(jié)構(gòu)來看，莫霍面往東翹起的上隆作用(可能與克拉通破壞和地幔對流有關(guān))和唐山斷裂下方的高密度體均說明上地幔巖漿的垂直上侵作用是唐山地震的主要垂直力源，莫霍面的上隆引起唐山斷裂深部發(fā)育有張性裂縫，上地幔巖漿沿張性裂縫上侵，在上地殼下部、中下地殼斷裂面附近形成巖漿囊，由于巖漿的上下對流維持高溫狀態(tài).這也是唐山地震前重力異常增加、垂直形變增加、地溫升高等前兆現(xiàn)象的原因(朱岳清，1984；李紹柄，1986).復(fù)雜的水平力與垂直力的聯(lián)合作用造就了唐山地震的孕育和發(fā)生，也引起唐山地震復(fù)雜的地表構(gòu)造變形形式(郭慧等，2011)，以及唐山斷裂淺部的花狀構(gòu)造(劉保金等，2011).綜上所述，我們提出唐山地震深部巖漿上侵與斷層聯(lián)合作用的發(fā)震模式如圖9所示.

圖9 唐山地震發(fā)震構(gòu)造環(huán)境模式示意圖(參考朱岳清，1984繪制)；地殼分層參考(段永紅等，2016)

4.3 結(jié)論

綜上所述，本文主要從重力異常與密度結(jié)構(gòu)角度探討了唐山地震的孕震構(gòu)造背景，主要取得以下認識：

(1) 利用華北及鄰區(qū)重力資料反演得到的莫霍面(圖5a)總體呈現(xiàn)由北西向南東逐漸減薄(隆起)的趨勢，可大體反映華北地區(qū)、唐山震區(qū)的地殼厚度分布，其與西太平洋板塊的俯沖、西伯利亞板塊南南東擠壓和印度板塊北東向碰撞或匯聚作用有關(guān)，是燕山運動、克拉通破壞等深淺構(gòu)造運動的總體體現(xiàn)，與地幔物質(zhì)向上運移和侵入有關(guān).

(2) 唐山震區(qū)三維精細密度結(jié)構(gòu)(圖7和圖8)顯示南正北負密度特征反映出華北盆地與燕山隆起構(gòu)造分區(qū)之間過渡的差異聯(lián)系，這種差異性反映出華北裂陷盆地與燕山褶皺的構(gòu)造應(yīng)力的不同，也是唐山地震形成的重要構(gòu)造動力背景或條件.

(3) 與唐山地震發(fā)震構(gòu)造唐山斷裂相關(guān)的高密度體特征(圖7和圖8)進一步凸顯斷裂與上地幔巖漿的侵入活動之間的聯(lián)系，結(jié)合層析成像與深地震測深結(jié)果分析，可能是板塊的俯沖引起上地幔巖漿物質(zhì)往上運移，沿深大斷裂上侵，使上地殼發(fā)震層熔融并弱化，引發(fā)地震，這可能是唐山地震形成的主因.該結(jié)果也為地震測深結(jié)果(曾融生等，1988；劉保金等，2011)提供了重力學(xué)證據(jù).

(4) 歷史地震主要分布在深大斷裂圍限的唐山菱形塊體內(nèi)部和中上地殼正密度異常區(qū)附近，且淺部的密度結(jié)構(gòu)顯示沿寶坻斷裂(F7)、潘莊西斷裂(F6)和寧河—昌黎斷裂(F1)兩側(cè)形成一正一負的局部隆起與拗陷.10～15 km深(圖7)高密度體呈倒R型，其應(yīng)與菱形斷裂分布有一定聯(lián)系(需要進一步研究)并說明深淺構(gòu)造存在較大差異，對唐山地震的形成具有控制作用，由于未看到與菱形斷裂有關(guān)的較深密度擾動，說明其作用有限.

致謝感謝馮旭亮老師和Camacho A G教授提供的界面反演及密度體反演程序，感謝野外數(shù)據(jù)采集者的辛勤付出.同時十分感謝匿名審稿專家的寶貴修改意見.