吳萍萍， 丁志峰* ， 譚捍東， 楊歧焱， 王興臣

1 中國地震局地球物理研究所， 北京 100081 2 中國地質(zhì)大學(北京)地球物理與信息技術學院， 北京 100083 3 河北地質(zhì)大學地球科學學院， 石家莊 050031

0 引言

華北克拉通中新生代以來經(jīng)歷了獨特的巖石圈改造和破壞過程，并伴隨著廣泛的巖漿活動和大規(guī)模的構造伸展盆地(翟明國等, 2005; 段永紅等, 2016)，是中國大陸內(nèi)部地震活動強烈地區(qū)之一(王椿鏞等, 2017).地震活動分布圖(圖1)反映華北地區(qū)的地震活動具有明顯的分帶性，主要有四條重要的地震帶：張家口—渤海地震帶(簡稱張渤地震帶)、河北平原地震帶、汾渭地震帶和鄂爾多斯西緣地震帶.其中張渤地震帶位于華北平原的北部，整體呈現(xiàn)北西向的地震活動帶特征，該區(qū)域分布著一系列雁行排列的北西至北西西向斷裂(圖2)，該斷裂帶與多條NNE-NE向斷裂帶交匯，M6.0以上的強震主要群集于NE向斷裂交匯地段.王椿鏞等(2017)統(tǒng)計該區(qū)域的發(fā)震深度主要集中在8～25 km范圍，在深度大于25km的地震數(shù)量急劇減少.深部介質(zhì)的電阻率對地球內(nèi)部溫度、流體、熔融/半熔融等異常反映靈敏，已有研究表明地震孕育和發(fā)生與地下高導體(例如火山高導巖漿囊、殼內(nèi)滑脫層、含水高導層等區(qū)域)和高阻體的接觸邊界息息相關(趙國澤等, 2001; 魏文博等, 2010; Zhang et al., 2016; Ye et al., 2018; 葉濤等, 2018).因此通過獲取張渤地震帶殼內(nèi)高導和高阻異常體的分布特征，可為認識研究區(qū)的物質(zhì)狀態(tài)、地殼運動過程、深部孕震環(huán)境等科學問題提供深部電性結構依據(jù).

大地電磁法以天然電磁場為場源，具有探測深度大、對低阻體反映靈敏、不受高阻屏蔽、設備輕便、成本低等優(yōu)勢，成為研究地球深部電性結構的主要地球物理方法之一.華北克拉通也開展了大量的大地電磁探測研究工作(魏文博等, 2006, 2008, 2010; Ouro-Djobo等, 2018; 詹艷等, 2011; 徐光晶等, 2015)，但因在張渤地震帶大地電磁測點數(shù)量少，深部電阻率結構成像分辨率有限.同時大地電磁法因受限于方法本身的局限性，采集易受電磁干擾、體積效應、縱向分辨隨深度增加而降低等因素影響，反演獲取的電阻率模型在深部電阻率經(jīng)常會出現(xiàn)明顯的拖尾現(xiàn)象，對于深部異常體的底界面形態(tài)的分辨率不足(金勝等, 2010).

已有研究表明基于交叉梯度結構耦合約束(Gallardo and Meju, 2003, 2004)的大地電磁和地震資料的聯(lián)合反演，可以實現(xiàn)兩種方法的優(yōu)勢互補，克服單方法的局限性，提高反演精度(彭淼等, 2013; Bennington et al., 2015; Ogunbo et al., 2018; Peng et al., 2019; Wu et al., 2020).但地震數(shù)據(jù)和大地電磁數(shù)據(jù)的聯(lián)合反演需要考慮兩種不同數(shù)據(jù)的耦合性，在聯(lián)合反演過程中需要權衡多種權重參數(shù)(例如，交叉梯度項權重因子、正則化因子、數(shù)據(jù)權重因子等)，在實際操作中存在困難.

目前華北地區(qū)已有很多地震資料的研究成果，包括人工地震深反射地殼結構(段永紅等，2016)、天然地震體波層析成像(黃金莉和趙大鵬, 2005; 楊婷等, 2012; 楊歧焱等, 2018)、地震面波和噪聲面波成像(魯來玉等, 2009; 何正勤等, 2009; Fang et al., 2010; 唐有彩等, 2011; 潘佳鐵等, 2011; 房立華等, 2013)、地殼厚度和泊松比(嵇少丞等, 2009; Lei, 2012)等.這些研究成果為認識該區(qū)域深部速度結構提供很好的素材.地下介質(zhì)速度和電阻率模型在一些特殊構造環(huán)境(例如局部熔融半熔融異常體、含水的裂隙或斷裂、殼內(nèi)滑脫層等)在三維結構分布上具有相似性(Newman et al., 2008; Berdichevsky et al., 2008)，根據(jù)這一特征可以相互驗證各自反演結果的可靠性，這也為本文基于地震學模型結構耦合約束的大地電磁二維反演算法開發(fā)提供思路.

彈性參數(shù)VP、VS、VP/VS、泊松比等屬性的空間分布特征對揭示地球內(nèi)部介質(zhì)的非均勻性有重要的意義，影響彈性參數(shù)的因素主要有礦物組成、溫度、壓力和空隙流體等.根據(jù)已有的實驗資料(陳颙等, 2009)可知，不同的地震波屬性對介質(zhì)溫度和成分變化的反應靈敏程度不同，例如當溫度發(fā)生變化時，S波速度對溫度異常的響應能力比P波速度強.VP/VS波速比模型(泊松比可以通過波速比換算得到)同時兼顧了介質(zhì)VP和VS屬性特點，可以提供比單純P波和S波速度參數(shù)更為豐富的信息(李志偉等, 2009).在熔融/半熔融體中，波速比值隨著熔融程度的增加而增大(劉瓊林等, 2011)，在這個地質(zhì)環(huán)境下電導率通常也表現(xiàn)出類似的特征.雖然很難定量分析波速比與電導率的物性關系式，但是可以根據(jù)這兩種物性在這些特殊環(huán)境中表現(xiàn)出來的結構相似性進行耦合約束反演，減小反演的多解性.

因此，為了避免地震資料和大地電磁資料聯(lián)合反演過程中多參數(shù)的權衡，同時充分利用已有的地震學模型，本文提出將VP/VS波速比模型作為固定模型約束，通過交叉梯度結構耦合函數(shù)帶入大地電磁二維反演算法中，實現(xiàn)基于地震學模型約束的大地電磁二維反演法，該方法即克服了聯(lián)合反演的局限，同時又可以有效克服大地電磁法深部電阻率分辨率不足的缺點.文中首先介紹大地電磁二維OCCAM(Constable et al.,1987)反演算法中加入VP/VS波速比結構耦合約束策略，通過理論模型合成數(shù)據(jù)檢驗算法的可靠性.然后將該方法應用于張渤地震帶大地電磁測深剖面，對比分析大地電磁法無約束和有約束反演結果，檢驗約束反演算法在張渤地震帶應用的優(yōu)越性和實用性.根據(jù)張渤地震帶反演結果設計一套理論模型，通過理論模型合成數(shù)據(jù)試算，檢驗算法在該研究區(qū)應用的有效性.最后結合該區(qū)域已有的地質(zhì)-地球物理資料，探討張渤地震帶深部電性結構、孕震環(huán)境及其動力學機制等問題.

1 地質(zhì)構造背景

圖1為華北地區(qū)大地構造背景及地震分布圖，華北克拉通主要由東、中和西三個主要構造區(qū)組成.西部為穩(wěn)定的鄂爾多斯高原，地震活動性主要集中在地塊的邊緣地區(qū)，地塊內(nèi)部表現(xiàn)為穩(wěn)定性特征，在塊體西緣(鄂爾多斯地震帶，圖1中標注為D的地震帶)曾在1739年發(fā)生寧夏平羅—銀川M8.0地震.中部地區(qū)為太行山造山帶，位于鄂爾多斯高原和華北平原之間，該區(qū)域構造復雜，地震頻發(fā)，曾發(fā)生M8.0以上地震3次(1303年山西洪洞地震、1556年陜西華縣地震、1695年山西臨汾地震).東部地區(qū)華北平原是研究華北克拉通巖石圈活化的熱點地區(qū)(Zhu et al., 2012; 翟明國, 2019)，該區(qū)域內(nèi)地震分布呈現(xiàn)兩條共軛地震條帶(圖1)，為北東向活動斷裂的河北平原地震帶(圖1中標注為B的地震帶)和北西向活動斷裂的張渤地震帶(圖1中標注為A的地震帶).本文研究區(qū)為張渤地震帶.

圖1 華北地區(qū)大地構造背景及地震分布圖A 張渤地震帶；B 河北平原地震帶；C 汾渭地震帶；D 鄂爾多斯地震帶. 紅色虛線表示主要斷裂分布.Fig.1 Geological background of the North China and the earthquakes distributionA Zhangbo seismic belt; B Hebei plain seismic belt; C Fenwei seismic belt; D Ordos seismic belt. Red lines donate the distribution of the main faults.

張渤地震帶西起太行山，橫貫燕山和華北盆地，東入渤海，是中國東部地區(qū)一條重要的北西向地震活動帶，主要發(fā)育NNE和NE向斷裂，部分斷層走向為NWW向.通過GPS資料分析，張渤地震帶以左旋走滑為主，兼有擠壓運動.華北平原塊體和燕山塊體相對運動時張渤帶左旋走滑的直接動力源(陳長云, 2016).

本文大地電磁測線沿著該地震帶布設(圖2)，自東向西經(jīng)過唐山、三河、北京、延慶到張家口地區(qū)，主要穿過9條斷裂帶，依次為寧河—昌黎斷裂(F1)、唐山斷裂(F2)、豐田—野雞坨斷裂(F3)、夏墊斷裂(F4)、黃莊—高麗營斷裂(F5)、南口斷裂(F6)、延慶—礬山盆地北緣斷裂(F7)、懷(來)—涿(鹿)盆地北緣斷裂(F8)和張家口斷裂(F9).其中在黃莊—高麗營斷裂(F5)是一條京西山前大斷裂，是華北平原和京西隆起的分界線，走向NE20°～50°，總長約110 km，斷裂面傾向南東，傾角55°～75°.

圖2 張渤地震帶構造背景和MT測點分布圖F1 寧河—昌黎斷裂; F2 唐山斷裂; F3 豐臺—野雞坨斷裂; F4 夏墊斷裂; F5 黃莊—高麗營斷裂; F6 南口斷裂; F7 延慶—礬山盆地北緣斷裂; F8 懷(來)—涿(鹿)盆地北緣斷裂; F9 張家口斷裂. 紅色五角星: 三河平谷M 8.0地震; 藍色圓M 7.0～M 8.0之間地震; 黃色圓M 6.0～M 7.0之間地震; 綠色圓表示M 5.0～M6.0之間地震. 01～31數(shù)字為大地電磁測點編號.Fig.2 Geological settings of Zhangbo seismic belt and the locations of the MT sitesF1 Ninghe-Changli fault; F2 Tangshan fault; F3 Fengtai-Yejituo fault; F4 Xiadian fault; F5 Huangzhuang-Gaoliying fault; F6 Nankou fault; F7 North margin fault of Yanqing-Fanshan basin; F8 North margin fault of Huailai-Zhoulu basin; F9 Zhangjiakou fault. Red pentagram: Shanhe-Pinggu M8.0 earthquake; blue circles: M 7.0～M 8.0 earthquakes; yellow circles: M6.0～M7.0 earthquakes; green circles: M5.0 ～M6.0 eaqthquakes. The numbers from 01 to 31 are the number of the MT sites.

測線最東邊為由寧河—昌黎斷裂(F1)、唐山斷裂(F2)、豐臺—野雞坨斷裂(F3)組成的唐山強震區(qū)斷裂構造.唐山斷裂帶是1976年唐山7.8級地震的發(fā)震構造，人工地震資料展示唐山斷裂帶走向NE30°，在地殼伴生復雜斷裂構造帶，表現(xiàn)為花狀構造斷裂特征，主干斷裂沿伸到Moho面(劉保金等, 2011a)，屬走滑斷層性質(zhì).夏墊斷裂(F4)是一條NNE向巖石圈尺度的區(qū)域性深斷裂帶，全長28 km，總體走向45°，傾向SE，傾角70～80°，是北京平原歷史上唯一發(fā)生過8級地震的活動斷裂.

測線最西側為張家口斷裂帶(F9)，屬于首都圈地震重點監(jiān)視防御區(qū)，為張家口-渤海地震構造帶與山西地震構造帶的交匯區(qū)(趙博等, 2011)，由北西西向和近東西向兩組、 多段斷裂組成， 北西西向斷層為斷裂主體，傾向南東，傾角60°左右，全長約30 km，屬正斷層性質(zhì).懷涿盆地北緣斷裂(F8)控制著懷來—涿鹿盆地的發(fā)育，總體走向北東，傾向南東，傾角50～70°，是重要的控盆活動正斷裂之一.延慶—礬山盆地北緣斷裂(F7)是一條控制延慶—礬山盆地北緣的主邊界正斷裂帶，總體走向NE，傾向SE，傾角50～80°，全長約102 km.

2 基于VP/VS波速比模型約束的二維大地電磁反演方法

2.1 交叉梯度項

Gallardo和Meju(2003, 2004)通過定義兩種不同物性參數(shù)的梯度的叉乘來定量描述兩種物性的結構相似性，這是本文實現(xiàn)基于VP/VS波速比模型約束的大地電磁二維反演算法的關鍵.根據(jù)定義可將電阻率和波速比模型的交叉梯度函數(shù)表達式寫成：

(1)

mr表示電阻率模型，ms表示VP/VS模型.根據(jù)求解旋度公式可將交叉梯度項分解為

(2)

(3)

(4)

ty=Wcg,s·mr.

(5)

Wcg,s為與網(wǎng)格和ms相關的常數(shù)矩陣. 當電阻率反演初始模型mr0為均勻半空間模型時，則Wcg,s·mr0=0，那么公式(5)可以改寫成：

ty=Wcg,s·(mr-mr0).

(6)

取公式(6)的平方和構建VP/VS波速比模型的結構耦合約束項，則公式可寫成：

(7)

2.2 基于VP/VS波速比模型約束的二維大地電磁反演目標函數(shù)

在大地電磁無約束反演目標函數(shù)中，通過權重因子η將VP/VS波速比模型的結構耦合約束項加入到無約束反演目標函數(shù)中，形成基于VP/VS波速比結構耦合約束的二維大地電磁反演目標函數(shù)，具體表示如下：

(8)

f(·)是二維大地電磁正演響應函數(shù)，λ為拉格朗日因子，Cd為數(shù)據(jù)協(xié)方差矩陣，Wm為模型光滑矩陣，dobs為觀測數(shù)據(jù)矩陣.對目標函數(shù)φ求梯度，則其梯度表達式可以寫成：

(9)

構建

(10)

根據(jù)OCCAM反演算法(Constable et al., 1987)的迭代公式，即可獲得基于VP/VS波速比模型約束反演第k次模型更新公式：

(11)

公式中λ因子的選取采用自動搜索算法實現(xiàn)，每一次迭代都要求根據(jù)RMS值選取最優(yōu)λ值，具體算法可參考Constable等(1987).圖3為本文開發(fā)的基于VP/VS波速比模型約束的大地電磁二維反演算法流程圖.

圖3 基于VP/VS波速比模型約束的大地電磁二維反演流程圖Fig.3 The flow chart of 2-D MT inversion algorithm with the VP/VS wave velocity ratio model constraints

3 理論模型合成數(shù)據(jù)反演試算

3.1 理論模型設計

為了檢驗基于VP/VS波速比模型約束的大地電磁反演算法的可靠性，本文設計一套復雜電阻率模型(圖4a)，其中包含有A、B、C1、C2、F電阻率塊體、近垂直構造體D和上涌構造體E.圖4a中黑色三角形為大地電磁測點，基于有限單元法(Constable et al.,1987)正演計算10-3Hz～1103Hz 區(qū)間42個頻點的TE和TM模式的視電阻率和相位，分別加入5%的高斯隨機誤差，作為觀測數(shù)據(jù).

圖4b為采用的VP/VS波速比固定模型，波速比值范圍為1.5～1.9.對比電阻率模型，波速比模型中設計了與電阻率模型結構相似(A、B、C1、F、D、E)和結構不相似(C2)的異常塊體或構造體.將VP/VS波速比模型(圖4b)作為固定模型約束，分別采用無約束和有約束大地電磁反演算法獲得電阻率結構，對比分析這兩種方法的反演結果，檢驗約束反演算法的可靠性和優(yōu)越性.

圖4 理論電阻率模型(a)和VP/VS波速比模型(b)(a) 理論電阻率模型; (b) 固定約束的VP/VS波速比地震學模型. A、B、C1、C2、F表示不同的物性塊體；D為近垂直構造體；E為上涌構造體.Fig.4 The resistivity model (a) and the VP/VS wave velocity ratio model (b)(a) Synthetic resistivity model; (b) The fixed VP/VS wave ratio model. A, B, C1, C2 and F represent the different physical blocks; D represents the near vertical structure; E represents the upwelling structure.

3.2 理論模型反演結果

圖5a、b分別為無約束和有約束反演結果.從圖5a中可以發(fā)現(xiàn)大地電磁法在無約束反演情況下，淺部電阻率結構具有較高的分辨率，例如塊體A的異常值和異常邊界都恢復得與真實模型相似、C1和C2塊體的上邊界接近于真實模型邊界、塊體B的電阻率值在淺部恢復得與真實值相近.但無約束反演結果在深部分辨率明顯不足，具體表現(xiàn)在：異常體C1和C2底邊界缺失，形成了與F貫通的兩個上涌構造體；近垂直構造體D在淺部有明顯的邊界信息，隨著深度增加，構造體的邊界信息減弱；上涌構造體E幾乎沒有體現(xiàn)出來.

圖5 理論模型無約束和有約束的大地電磁二維反演結果(a) 無約束大地電磁反演結果; (b) 有約束大地電磁反演結果; (c) 無約束反演電阻率模型與VP/VS波速比模型的交叉梯度值空間分布; (d) 有約束反演電阻率模型與VP/VS波速比模型的交叉梯度值空間分布. RMS為數(shù)據(jù)擬合差.Fig.5 The results derived from the unconstrained and constrained inversion method(a) The resistivity model derived from the unconstrained inversion method; (b) The resistivity model derived from the constrained inversion method; (c) The spatial distribution of the cross-gradient value between the VP/VS ratio model and the unconstrained inverted resistivity model; (d) The spatial distribution of the cross-gradient value between the VP/VS ratio model and the constrained inverted resistivity model. RMS represents the data residual.

圖5b為大地電磁法約束反演結果，從圖中可以看出電阻率反演結果在深度分辨率有明顯提高.在電阻率模型和VP/VS波速比模型結構相似處電阻率邊界信息有明顯改善，具體表現(xiàn)在：C1和C2塊體由原來的連通狀態(tài)恢復成兩個獨立圈閉的塊體，其中C1的邊界信息與VP/VS波速比模型相似，在C1的異常值和邊界都恢復得與真實模型相似，而在C2位置處，因在VP/VS波速比模型中沒有存在C2異常，雖然C2異常值有改善，但邊界信息與真實情況仍有差距；上涌構造E和近垂直構造D的邊界信息恢復得與真實模型相近.圖5c和d分別為無約束反演結果、有約束反演結果與VP/VS波速比模型的交叉梯度值空間分布圖，從圖中可以發(fā)現(xiàn)在塊體C1周邊、近垂直構造D、上涌構造E的交叉梯度值有明顯降低，這說明經(jīng)過約束反演后的模型，在這幾個區(qū)域的結構更相似.從反演結果中不難發(fā)現(xiàn)對于淺部電阻率值和形態(tài)(例如塊體A)在無約束和有約束反演都可以很好地恢復，這說明基于地震學模型約束的大地電磁反演算法可以在滿足數(shù)據(jù)擬合差的前提下，獲取電阻率模型和地震學模型結構上更為相似的模型，有效提高大地電磁法深部電阻率分辨率.

4 張渤地震帶大地電磁測深剖面

4.1 大地電磁數(shù)據(jù)采集和處理

圖2為張渤帶大地電磁測點布設圖，測線跨越張家口、北京、唐山等城市，全長400 km，測線方向為北偏西54°.測點布設采用“+”字形正南北東西向布極方式，采集加拿大鳳凰公司的MTU-5A寬頻帶大地電磁觀測設備.從該測線中挑選出數(shù)據(jù)質(zhì)量較高的31個MT測點，分別對每個測點數(shù)據(jù)進行濾波、重采樣、頻譜分析等預處理，基于WinGLink平臺獲取0.001 Hz到320.0 Hz的TE模式和TM模式的視電阻率和相位曲線.同時采用Berdichevskiy和Rhoplus方法檢驗視電阻率和相位曲線的一致性，將處理后的數(shù)據(jù)作為反演數(shù)據(jù).

圖6中黑色曲線為TE模式的視電阻率和相位，紅色曲線為TM模式視電阻率和相位.根據(jù)地質(zhì)構造特征，將測線上MT測點數(shù)據(jù)分成四段，分別為唐山斷裂帶區(qū)域(圖6a，紅色測點)、三河—平谷斷裂帶區(qū)域(圖6a，黃色測點)、華北平原—京西丘陵過渡區(qū)(圖6a，藍色測點)和京西丘陵區(qū)域(圖6a，黑色測點).從圖中可以看出除了華北平原—京西丘陵過渡帶的視電阻率和相位曲線形態(tài)差異比較大之外，其他三個區(qū)域的數(shù)據(jù)形態(tài)較為一致.其中唐山斷裂帶區(qū)域視電阻率曲線形態(tài)和京西丘陵區(qū)域曲線形態(tài)相似，但京西丘陵段的最高值出現(xiàn)在0.15 Hz左右，而唐山斷裂帶區(qū)域最高值出現(xiàn)在0.06 Hz左右，這兩段的相位曲線差異比較大，說明地下電性結構山區(qū)和平原存在差別.三河—平谷斷裂帶視電阻率值比較低，最高值約300 Ωm，該區(qū)域的相位曲線形態(tài)與唐山斷裂帶相似.而華北平原—京西丘陵過渡帶的視電阻率和相位曲線形態(tài)復雜，沒有明顯的規(guī)律，這也間接說明該區(qū)域地下電性結構復雜的特征.

圖6 實測大地電磁資料視電阻率和相位曲線(a) 紅色三角形為Segment 1測點位置; 黃色三角形為Segment 2測點位置; 藍色三角形為Segment 3測點位置; 黑色三角形為Segment 4測點位置. 黑色曲線為TE模式的視電阻率和相位; 紅色曲線為TM模式的視電阻率和相位.Fig.6 Measured MT apparent resistivity and phase data(a) Red triangles represent the MT sites corresponding to the data of segment 1; Yellow triangles represent the MT sites corresponding to the data of segment 2; Blue triangles represent the MT sites corresponding to the data of segment 3; Black triangles represent the MT sites corresponding to the data of segment 4. Black curves represent the apparent resistivity and phase data of TE mode; red curves represent the apparent resistivity and phase data of TM mode.

4.2 VP/VS波速比模型

楊歧焱等(2018)基于華北地區(qū)數(shù)字地震臺網(wǎng)176個固定臺站，采用近震體波層析成像法獲取了張渤地震帶VP/VS波速比地震學模型.本文提取了張渤地震帶對應大地電磁測線剖面的VP/VS波速比模型(圖7)，對該剖面進行平滑，去除模型的突變點，然后將該模型作為固定約束代入公式(5—11)中的Wcg,s，實現(xiàn)基于VP/VS波速比模型約束的大地電磁法二維反演.

圖7 張渤地震帶VP/VS波速比剖面圖TS 唐山； SH 三河； BJ 北京； ZJK 張家口.Fig.7 The VP/VS wave velocity ratio profile along Zhangbo seismic beltTS Tangshan; SH Shanhe; BJ Beijing; ZJK Zhangjiakou.

4.3 反演網(wǎng)格

二維大地電磁反演網(wǎng)格剖分中，用均勻和非均勻剖分法分別對目標研究區(qū)域和邊界區(qū)域進行網(wǎng)格化處理(圖8).在目標研究區(qū)采用均勻網(wǎng)格剖分，X軸網(wǎng)格間距為3 km，Z軸網(wǎng)格間距為2 km.在邊界網(wǎng)格采用的是指數(shù)遞增網(wǎng)格間距.圖中8紅色框為VP/VS波速比模型約束區(qū)，在反演過程中只在紅色框區(qū)域內(nèi)進行結構耦合約束.

圖8 大地電磁正、反演網(wǎng)格剖分黑色線為網(wǎng)格；黑色方框為約束反演區(qū)域.Fig.8 Mesh discretization for MT forward and inversionBlack lines represent the mesh; black rectangle shows the constrained inversion region.

4.4 張渤地震帶反演結果

圖9a為無約束和有約束的大地電磁法反演迭代RMS曲線圖，從圖中可知這兩種方法都具有較好的收斂性.從RMS迭代變化趨勢可以看出，當加入VP/VS波速比模型約束項時，RMS值迭代下降的速度變慢，但經(jīng)過多次迭代后，有約束和無約束反演結果的RMS都收斂到相同水平上.

圖9 RMS迭代曲線和約束項權重因子L-Curve曲線(a) 無約束和有約束反演RMS迭代曲線; 空心圈曲線為無約束反演(SP)迭代RMS曲線; 實心點曲線為有約束反演(CG)迭代RMS曲線. (b) 約束項權重因子L-Curve曲線; 實心圈曲線為不同權重因子的RMS值; 黑色箭頭標注為本文選取的權重因子.Fig.9 Plots of RMS iteration curves and L-curve analysis for selecting constrained weight factors(a) RMS iteration curves for no-constrained and constrained inversion. The hollow circle curve represents the RMS iteration curve for unconstrained inversion. The solid dot curve represents the RMS iteration curve for constrained inversion. (b) L-curve anasys for selcting constrained weight factors; The solid dot curve represents the RMS value for different weighting factors; the black arrow represnets the optimal weighting value.

圖9b為公式中VP/VS波速比模型約束項權重因子(公式8中η)的選擇.本文采用“L-Curve”方法，選取了多個權重因子進行二維大地電磁法約束反演(圖9b)，圖中黑色箭頭為本文最終確定的權重因子.圖10展示了6個測點無約束和有約束反演模型響應曲線，從圖中可以發(fā)現(xiàn)無約束和有約束反演模型響應曲線基本上重合，與觀測數(shù)據(jù)有較好的一致性，說明這兩個模型都可以很好地擬合實測數(shù)據(jù).

圖10 無約束和有約束反演模型在測點01、09、16、19、26和30處的正演響應紅色點: 實測數(shù)據(jù)(real); 黑色圈: 無約束反演結果正演響應(SP); 藍色叉: 有約束反演結果正演響應(CG).Fig.10 The synthetic data sets of the unconstrained and constrained inverted model at the sites of 01, 09, 16, 26, and 30.Red dots represent the measured data (real); black circles represent the synthetic data sets of the unconstrained inverted model (SP); blue crosses represent the synthetic data sets of the constrained inverted model (CG).

圖11展示了大地電磁法二維無約束和有約束反演結果，圖11a為測線對應的地形曲線，可以發(fā)現(xiàn)測線的西部測點主要位于京西丘陵區(qū)，測線的東部主要位于華北平原地區(qū)，測線上由東向西穿越了唐山斷裂帶、三河—平谷斷裂帶、華北平原—京西丘陵過渡帶、延懷盆地和張家口斷裂帶.

圖11 大地電磁法無約束和有約束反演結果(a) 地形剖面圖; (b) 無約束反演電阻率模型; (c) 有約束反演電阻率模型. H1—H5為高阻異常體編號; L1—L3為低阻異常體編號.Fig.11 The resistivity models derived from the unconstrained and constrained inversion method(a) Topographic profile; (b) The resistivity model derived from the unconstrained inversion method; (c) The resistivity model derived from the constrained inversion method. H1—H5 represent the numbered high resistivity anomalies; L1—L3 represent the numbered low resistivity anomalies.

為了方便對比這兩種方法反演結果，將圖中5個高阻體標注為H1、H2、H3、H4和H5，3個低阻區(qū)域標注為L1、L2和L3.從圖中可以看出，基于VP/VS波速比模型約束的反演結果(圖11c)的橫向分辨率和縱向分辨率都有明顯的提高，主要表現(xiàn)有：在唐山強震區(qū)的高阻異常體H5沿伸到下地殼，經(jīng)過約束反演后，高阻異常的拖尾現(xiàn)象有明顯的改善，底部邊界約在20 km左右；從夏墊斷裂帶到京西山前斷裂帶的高阻體異常體H4近垂直地沿伸到地幔，很難識別出該異常體在地殼內(nèi)的異常形態(tài)，經(jīng)約束反演后，深部的分辨率有明顯的提升，高阻異常體H4在地殼中表現(xiàn)為圈閉異常形態(tài)；在京西丘陵區(qū)域的地殼高阻異常體H1、H2和H3在無約束反演中表現(xiàn)為整體高阻異常，深度延伸到下地殼，經(jīng)約束反演后，可以明顯看出H1、H2和H3這三個高阻異常體是獨立分開的三個高阻異常，同時這三個異常體的底界面大約在20 km；經(jīng)過約束反演后，發(fā)現(xiàn)深部高導區(qū)域L1和L2從30～40 km上涌，這個現(xiàn)象在無約束反演結果中沒有發(fā)現(xiàn)；在夏墊斷裂帶區(qū)域無約束和有約束反演都可以發(fā)現(xiàn)明顯的低阻異常體L3，但低阻異常的形態(tài)、規(guī)模有很大的差異.對比VP/VS波速比模型(圖7)和大地電磁法約束反演結果(圖11c)，在20 km范圍內(nèi)，電阻率模型與VP/VS波速比模型的異常體表現(xiàn)出明顯的相似性.而在大于20 km深度，這兩個模型存在比較大的差異，表現(xiàn)在約30 km處波速比模型出現(xiàn)低波速比異常層，而電阻率模型上未發(fā)現(xiàn)該特征.

圖12為大地電磁法無約束和有約束反演的電阻率模型與VP/VS波速比模型的交叉梯度值空間分布圖，mean|ty|為約束反演區(qū)域(圖8黑色框)中每個網(wǎng)格的交叉梯度值的絕對值的平均值.從圖中可以發(fā)現(xiàn)有約束反演結果與VP/VS波速比模型的交叉梯度值比無約束反演結果的交叉梯度值有明顯下降(mean|ty|由2.94降到1.5)，這說明基于VP/VS波速比模型約束的大地電磁反演可以在滿足數(shù)據(jù)擬合差的前提下，可以獲取與地震學模型結構上更為耦合的電阻率模型.

圖12 無約束、有約束反演電阻率模型與VP/VS波速比模型的交叉梯度值分布圖(a) 無約束反演電阻率模型與VP/VS波速比模型的交叉梯度值空間分布; (b) 有約束反演電阻率模型與VP/VS波速比模型的交叉梯度值空間分布. mean|ty|為約束反演區(qū)域(圖8紅色框)中每個網(wǎng)格的交叉梯度值的絕對值的平均值.Fig.12 The spatial distribution of the cross-gradient value between the VP/VS ratio model and the non-constrained or the constrained inverted resistivity model.(a) The spatial distribution of the cross-gradient value between the VP/VS ratio model and the unconstrained inverted resistivity model; (b) The spatial distribution of the cross-gradient value between the VP/VS ratio model and the constrained inverted resistivity model. mean |ty| represents the average of the absolute value of the cross-gradient values in the common region (red box in Fig.8).

根據(jù)上述現(xiàn)象，不難看出要使得基于地震學模型約束反演算法能有效改善大地電磁反演結果的前提條件是地震學模型和電阻率模型存在結構相似性.當兩種物性結構相似時，基于交叉梯度的約束反演可以有效地加強兩種物性結構耦合關系，提高結構相似處的分辨率.在結構不相似的區(qū)域，交叉梯度項幾乎不起作用，保留單方法反演的優(yōu)勢.Wu 等(2020)、Peng 等(2019)通過大量的理論模型算例充分地論證這一特征.

4.5 張渤地震帶約束反演結果可靠性和有效性檢驗

通過對比無約束、有約束的大地電磁反演結果，可知在地震學模型的約束下，張渤地震帶深部電性結構在橫向和縱向的分辨率都有明顯的提升.為了檢驗基于VP/VS波速比模型約束的大地電磁反演算法在張渤地震帶應用的可靠性和有效性，本文考慮設計形態(tài)與反演結果相似的5個獨立的高阻體(圖11中H1—H5)和3個上涌的低阻構造體(圖11中L1—L3)，通過理論模型合成數(shù)據(jù)試算，采用與實測數(shù)據(jù)反演一樣的VP/VS波速比模型(圖7)進行約束反演，對比無約束和有約束反演結果，從理論模型上檢驗該算法是否能提高張渤地震帶深部電阻率橫向和縱向分辨率.

為了能滿足上述5個獨立的高阻體和3個上涌的高導構造體的理論模型需求，本文將約束反演結果(圖11c)中電阻率值大于103.5Ωm的電阻率值設置為恒定值103Ωm，將小于101.8Ωm的電阻率值設置為恒值10 Ωm，其他區(qū)域的值設置為102Ωm.圖13a為本文設計的理論電阻率模型，從圖中可以看出理論模型在20km以內(nèi)設計有5個圈閉的高阻異常體(圖13a 中H1—H5)，中下地殼3個低阻隆起(圖13a 中L1—L3).

測點位置(圖13a黑色三角形)和每個測點的頻段選擇與實際觀測一樣，正反演網(wǎng)格也跟上述一樣.正演每個測點的理論視電阻率和相位曲線，同時加入與實測數(shù)據(jù)相同的誤差，形成觀測數(shù)據(jù).

圖13 理論電阻率模型(a)、無約束(b)和有約束(c)反演結果(a) 理論電阻率模型; (b) 無約束反演電阻率模型; (c) 有約束反演電阻率模型. H1—H5為高阻異常體編號; L1—L3為低阻異常體編號.Fig.13 Synthetic model (a), unconstrained (b) and constrained (c) inverted models(a) Synthetic model; (b) Unconstrained inverted model; (c) Constrained inverted model. H1—H5 represent the numbered high resistivity anomalies; L1—L3 represent the numbered low resistivity anomalies.

圖13b和c分別為大地電磁法二維無約束和有約束的反演結果，對比這兩種方法的反演結果，可以發(fā)現(xiàn)無約束反演的結果與上述實際數(shù)據(jù)無約束反演結果出現(xiàn)的問題相似，主要表現(xiàn)在：獨立高阻異常體H1—H4在無約束反演(圖13b)中表現(xiàn)為連接在一起、規(guī)模大的高阻異常；低阻體L1和L2幾乎沒有體現(xiàn)出來，L3的低阻異常體形態(tài)與真實異常體形態(tài)差別很大；同時異常體的底界面相比真實模型有明顯的往下拖尾的現(xiàn)象.從無約束反演結果可以看出因大地電磁法受限于方法本身體積效應、深度分辨率隨深度增加而降低等局限，無約束反演獲得的電阻率模型在橫向和垂向分辨率有限.

圖13c為基于VP/VS波速比模型約束反演的結果，從圖中可以明顯的看出淺部高阻異常體(H1—H5)的橫向分辨率和縱向分辨率都有明顯的提高，5個高阻異常體可以清晰的展示出獨立異常體的特征，異常體的拖尾現(xiàn)象得到明顯的改善；對于低阻體L1、L2和L3異常體的反演精度有明顯的提高.因此，通過這個理論模型合成數(shù)據(jù)試算，驗證了基于VP/VS波速比模型約束的二維大地電磁反演法在張渤地震帶應用可以提高深部電阻率的橫向和縱向分辨率，從而驗證了該算法在張渤地震帶應用的可靠性和有效性.

5 討論

根據(jù)反演得到的張渤地震帶地殼電性結構，結合地質(zhì)資料繪制電性構造解釋圖(斷裂信息參考徐錫偉等(2002)).圖14中斷裂F1—F9的深部形態(tài)分別用黑色虛線展示，斷裂地表位置與圖2中斷裂(F1—F9)位置對應.從整體上來看，電阻率結構表現(xiàn)出明顯的橫向不均勻性，與地表的各個構造相呼應.對比圖11a地形曲線和構造分區(qū)，可以發(fā)現(xiàn)測線上斷裂帶在電性結構上表現(xiàn)為明顯的電性分界面或電性梯度帶，在張家口斷裂帶、延慶—懷柔盆地、夏墊斷裂帶和唐山斷裂帶在中上地殼普遍存在高阻異常體，而在中下地殼表現(xiàn)為低阻異常特征.從M≥3.0地震分布特征來看，發(fā)現(xiàn)該區(qū)域的地震活動與斷裂帶的分布息息相關，具有明顯的分段性，以下將分段討論剖面電性結構特征.

唐山震群主要有三條斷裂帶控制(寧河—昌黎斷裂F1、唐山斷裂F2和豐臺—野雞坨斷裂F3).唐山7.8 級地震發(fā)生在主斷裂帶(F2)上，該區(qū)域為高阻塊體與低阻塊體的相交地帶、偏向高阻體一側.根據(jù)人工地震反射資料可知唐山斷裂為主斷裂，延伸至Moho面，地殼伴生的斷裂構造表現(xiàn)為走滑斷裂(花狀構造)特征.圖中可以發(fā)現(xiàn)高阻體的底界面大約在20 km左右，地震的發(fā)生主要集中在20 km以內(nèi)，位于高阻區(qū)域或高低阻交界區(qū).人工深反射剖面(劉保金等, 2011a)和天然地震觀測(劉啟元等, 2007; 楊歧焱等, 2018)發(fā)現(xiàn)唐山震區(qū)在深度<20 km內(nèi)呈現(xiàn)高速異常，在約深20 km處存在明顯的低速層(體)，該低速層(體)一直延伸到40 km以下的深度.劉國棟等(1983)發(fā)現(xiàn)唐山地震源區(qū)下約20 km深處存在高導異常.嘉世旭和張先康(2005)通過深地震測深資料也發(fā)現(xiàn)唐山強震區(qū)中下地殼存在低速異常.巖漿底侵和置換作用是下地殼和地幔物質(zhì)交換、相互作用的重要形式，推測唐山震區(qū)中下地殼的低速低阻異常體與流體和熱物質(zhì)作用有關，大規(guī)模的幔源物質(zhì)侵入，形成了中上地殼的孕震環(huán)境.

在三河—平谷斷裂帶區(qū)域，存在明顯的殼淺部到深部低阻特征，中下地殼存在上涌低阻體(L3)，淺部和深部的高導體不連續(xù)，地震分布主要位于該不連續(xù)層中，說明該區(qū)域深淺高導異常的成因可能不同.鄧前輝等(2001)通過大地電磁法也揭示夏墊斷裂帶中下地殼高導體和地幔上涌現(xiàn)象，淺部和深部的構造存在差異，推測震源體的形成和發(fā)震與深部高角度隱伏斷裂有關.對應地質(zhì)資料可知，三河—平谷斷裂帶為第四紀活動隱伏高傾角斷裂帶，地表覆蓋有上新統(tǒng)-第四系沉積層，淺部低阻異常與沉積層對應.從人工地震深反射剖面(劉保金等, 2011b; 趙金仁等, 2004)可知該斷裂帶是一條深淺共存的斷裂構造帶，在Moho存在復雜的楔形反射帶，說明該區(qū)域地殼發(fā)生過強烈的擠壓、變形，同時也反映出巖漿活動對地下地殼結構經(jīng)過強烈改造，形成了復雜的地殼深部結構.徐錫偉等(2002)也從人工地震反射剖面發(fā)現(xiàn)Moho面存在反射過渡帶，推測是早期地殼伸展過程中形成斷裂帶，填充了來自上地幔的基性巖墻.從圖14中深部上涌低阻異常(L3)可能反映的是上地幔熱物質(zhì)侵入，Zhu等(2012, 2015)指出西太平洋板塊的俯沖時，俯沖板塊的脫水交代上覆地幔楔使其發(fā)生熔融，俯沖帶后撤后導致巖石圈強烈伸展，導致地幔熱物質(zhì)上涌，形成了中上地殼的孕震環(huán)境.

黃莊—高麗營斷裂帶(F5)和南口斷裂(F6)為京西丘陵的山前斷裂帶，是大型拆離斷裂帶，在電阻率結構上表現(xiàn)為明顯的電阻率變化梯度帶.在山前斷裂帶到夏墊斷裂帶之間淺部有幾公里的反映沉積層特征的低阻層，在低阻層下方、延伸到下地殼為規(guī)模較大的高阻異常，未見低速體上涌現(xiàn)象.詹艷等(2011)從石家莊地區(qū)大地電磁剖面也發(fā)現(xiàn)太行山山前斷裂從淺地表到下地殼表現(xiàn)為高電阻體，無殼內(nèi)高導層發(fā)育.京西丘陵的山前斷裂帶是華北克拉通東部與中部的邊界，在重力梯度帶、地形、巖石圈厚度等資料顯示在該區(qū)域附近變化顯著(Zhu et al., 2012).規(guī)模大、延伸深的高阻異常(H4)可能是華北克拉通東部和中部地區(qū)構造運動形成的高密度的輝長巖、鐵鎂質(zhì)巖石或者下地殼麻粒巖(翟明國等, 2005; 翟明國, 2019).華北平原地區(qū)強震主要發(fā)生在NE和NW兩組構造帶，黃莊—高麗營斷裂帶(F5)和南口斷裂(F6)山前斷裂帶自新生代以來對北京凹陷起著重要的控制作用(焦青等, 2005;焦青和邱澤華, 2006)，是北京地區(qū)值得關注的強震發(fā)生重點監(jiān)視防御區(qū).

在南口斷裂(F6)以西區(qū)域為京西丘陵地區(qū)，主要有延慶—懷來盆地和張家口斷裂帶.延懷盆地是由山西裂谷走滑運動形成的NE向次級張性構造區(qū)，地表延慶—礬山次級盆地北緣斷裂(F7)和涿鹿—懷來盆地北緣斷裂(F8)控制盆地構造特征.在電阻率結構成像中(圖14)，在淺部表現(xiàn)明顯的低阻異常，在中下地殼有明顯的低阻體沿著盆地下方通道侵入的現(xiàn)象.從地震反射剖面(張先康等, 1996)可知在延慶—礬山盆地北緣斷裂正下方存在低速體，其成因可能與該區(qū)域地殼伸展運動、幔源物質(zhì)上涌有關，在地表形成盆嶺構造.

圖14 張渤地震帶電性結構解釋圖白色圈為2007年至今M≥3.0地震(由中國地震臺網(wǎng)中心提供). F1 寧河—昌黎斷裂; F2 唐山斷裂; F3 豐臺—野雞坨斷裂; F4 夏墊斷裂; F5 黃莊—高麗營斷裂; F6 南口斷裂; F7 延慶—礬山盆地北緣斷裂; F8 懷(來)—涿(鹿)盆地北緣斷裂; F9 張家口斷裂. F1—F8走向為NE或NEE；F9走向為NW.Fig.14 The geology interpretation maps of the electrical structure beneath Zhangbo seismic beltWhite circles: M≥3.0 earthquakes since 2007 (provided by China earthquake Networks Center). F1 Ninghe-Changli fault; F2 Tangshan fault; F3 Fengtai-Yejituo fault; F4 Xiadian fault; F5 Huangzhuang-Gaoliying fault; F6 Nankou fault; F7 North margin fault of Yanqing-Fanshan basin; F8 North margin fault of Huailai-Zhoulu basin; F9 Zhangjiakou fault. The strike of F1—F8 is NE or NEE; the strike of F9 is NW.

張家口斷裂帶位于張渤地震帶和山西地震帶的交匯區(qū)(趙博等, 2011)，區(qū)域內(nèi)斷裂相交和切割，形成復雜的孕震構造格局，主要由北西西向和近東西向兩組、 多段斷裂組成， 北西西向斷裂為斷裂帶主體(圖14中F9)，與研究區(qū)數(shù)條北北東至北東向的斷裂顯著不同(高戰(zhàn)武等, 2001)，是首都圈地震重點監(jiān)視防御區(qū).張學民等(2012)通過地震體波層析成像發(fā)現(xiàn)太行山北段多構造交界的張家口地區(qū)存在局部低速異常，推斷該區(qū)域受到西太平洋板塊俯沖的影響，表現(xiàn)為拉張應力場，深部地幔熱物質(zhì)活動及淺部應力累積引發(fā)活躍的構造活動.地質(zhì)資料顯示張家口地區(qū)出露中生代形成的源于古老下地殼或地幔的堿性花崗巖(楊進輝等, 2006; 韋忠良等, 2008; 汪洋和程素華, 2010)，從巖石學、地球化學角度論證了該區(qū)域經(jīng)歷了大規(guī)模的殼幔相互作用.從本文獲取的電性結構圖上可以發(fā)現(xiàn)在該斷裂帶上深部有低阻體上涌，在張家口中下地殼高導層的規(guī)模比懷來盆地大，這可能暗示張家口斷裂帶的深部物質(zhì)作用更為強烈.

6 結論

本文實現(xiàn)了基于VP/VS波速比模型約束的大地電磁二維反演算法，通過理論模型合成數(shù)據(jù)試算，檢驗了算法的可靠性.將約束反演算法應用于張渤地震帶31個大地電磁實測資料，研究發(fā)現(xiàn)唐山斷裂帶中上地殼表現(xiàn)為高阻異常，地震主要分布于高阻異常體周邊，在下地殼底部有上涌的低阻異常體，推測該區(qū)域深部低阻異常體與幔源物質(zhì)侵入有關.三河—平谷斷裂帶淺部表現(xiàn)為沉積層的低阻異常特征，深部高導異常與夏墊斷裂帶展布、殼幔物質(zhì)相互作用有關，地震事件主要發(fā)生在淺部和深部高導異常不連續(xù)區(qū)域.太行山山前斷裂電阻率結構上表現(xiàn)為明顯的電阻率變化梯度帶，在與華北平原過渡帶中下地殼存在規(guī)模大、延伸深的高阻異常，推測與該區(qū)域因構造運動形成高密度的輝長巖、鐵鎂質(zhì)巖石或者下地殼麻粒巖有關.延懷盆地在淺部表現(xiàn)明顯的低阻異常，在中下地殼有明顯的低阻體沿著盆地下方通道侵入的現(xiàn)象，其成因可能是地殼伸展運動，幔源物質(zhì)上涌，形成地表的盆嶺構造.張家口斷裂帶中下地殼高導異常層比懷來盆地深部高導異常規(guī)模大，這可能暗示張家口斷裂帶的深部物質(zhì)作用更為強烈.

致謝感謝中國地質(zhì)大學(北京)王瑜教授和中國地震局地球物理研究所尤惠川研究員對本文的地質(zhì)構造解釋提供的幫助.感謝審稿專家和編輯部對本文提出的重要修改意見.