京津冀地區地殼應力場特征

王曉山 馮向東 趙英萍

摘要：基于標量斷層類型值，對京津冀地區及鄰區2 187個中小地震震源機制解進行分類，統計結果顯示研究區震源機制類型以走滑斷層和正斷層為主，P軸優勢方位為NEE—EW和SWW—EW向;采用MSATSI軟件包反演該區1°×1°網格的精細地殼應力場，結果表明：最大主壓應力軸最優解的優勢方向為NEE—EW向，與P軸優勢方位一致;所有網格的相對應力大小R值均小于0.5，表明京津冀地區應力狀態偏拉張性質，而且最小主壓應力軸的不確定度變化范圍相對穩定，表明現今京津冀地區地殼應力場處于一個相對統一的NNW—SSE向的拉張作用控制下。39°N以北地區最大主壓應力軸方位最優解顯示一定角度的偏轉，同時最大、中等、最小主壓應力軸最優解推斷的應力狀態由西向東存在一個正斷層—走滑斷層—正斷層的轉換過程;而39°N以南地區的現今構造應力場保持穩定，最優主壓應力軸呈NEE—SWW向，大部分網格應力狀態顯示走滑型。構造應力場的反演結果與活動構造、GPS主應變方向和剪切波分裂的快波偏振方向等相關研究結果基本一致。

關鍵詞：京津冀地區;震源機制;應力張量阻尼反演;地殼應力場

中圖分類號：P315.727文獻標識碼：A文章編號：1000-0666（2020）04-0610-10

0引言

地殼應力場的時空分布特征，特別是地殼深部應力狀態，是解決地震發生機理、區域地殼穩定性、震源物理過程等地球動力學有關科學問題的基礎。地質調查與大地測量結果表明，華北平原現今不存在大規模的拉張變形，區內主要斷裂帶的滑移速率低（Shen et al，1996;郭良遷等，2008），是一個地殼變形相對較弱的地區。在這種緩慢的變形區域中，沉積蓋層或淺部巖石中的直接應力測量通常不能代表深部的應力場，而震源機制解是可用于估計緩慢變形區域中震源深度處應力場的少數信息之一。

京津冀地區是我國東部強震活動頻度和強度最強的地區，同時也是新構造活動最強烈的地區之一，區內分布有張家口—渤海地震構造帶（張渤地震帶）和華北平原地震構造帶，中強地震絕大多數都發生在這2條構造帶或其交會部位。該區是國內使用地震記錄開展地殼應力場較早的地區，李欽祖等（1973）利用河北紅山臺和沙城臺記錄的小地震P波初動資料確定了臺站周邊區域的應力場;許忠淮等（1979，1983）對上述方法進行改進，形成了“綜合震源機制解法”，求解了京津唐張和華北地區的構造應力場方向。首都圈數字地震臺網于2002 年開始正式運行（劉瑞豐等，2008），為地震精確定位和震源機制解反演提供了有利條件，研究人員利用P波初動極性、振幅比和波形反演求解了大量的中小震震源機制解，且其質量也更加可靠（胡新亮等，2004）。基于此，不同研究人員使用不同方法開展了包括首都圈在內的京津冀地區的地殼應力場研究（蘭從欣等，2005;李瑞莎等，2008;張紅艷等，2009;武敏捷等，2011，2012;胡幸平，崔效鋒，2013;黃驥超，萬永革，2015;劉靜等，2016;劉麗等，2017;樊文杰等，2019）。其中大部分研究根據構造特征和地震分布將研究區域分成幾個應力小區，每個應力小區的范圍都和許忠淮等（1983）對京津冀地區的分區基本一致，所不同的是張紅艷等（2009）在此基礎上進行了細分，但是所得的構造應力場方向基本相同;劉靜等（2016）和劉麗等（2017）對河北地區整體構造應力場的研究結果也驗證了該地區處于華北地區統一地殼應力場的控制下（李欽祖，1980;魏光興等，1982;許忠淮等，1983）。

在京津冀地區震源機制解數據日益豐富的當下，武敏捷等（2011）采用細分網格加滑動步長的方式初步討論了華北北部地區的構造應力場，由于每個網格都是獨立反演，分區方式不同和分區大小的變化都可能對反演的應力張量結果產生影響（Maury et al，2013）。相鄰2個網格主應力方向的不同或許反映了應力場真實的變化，但也可能是震源機制數據的誤差或反演的約束較差造成的假象。地殼中的真實應力場是連續分布的，為了解決應力反演模型依賴于分區的問題，更準確地反映研究區域的應力場空間變化特征，Hardebeck和Michael（2006）構建了一組可調整的阻尼參數，引入平滑約束來抑制相鄰網格應力模型的差異，提出了區域尺度的應力張量阻尼反演方法。本文基于京津冀地區及鄰區的中小地震震源機制解，將該方法應用于對該區精細地殼應力場特征的研究。

1數據和方法

1.1震源機制解數據

本文所用震源機制解共2 187個，震級區間2.0≤ML≤5.5，覆蓋范圍包括京津冀地區以及山西、遼寧、內蒙部分地區。數據來源分為2個部分：第一部分為2002—2008年京津冀地區ML≥2.0地震的震源機制解，使用梁尚鴻等（1984）提出的利用區域地震臺網垂直向P，S最大振幅比資料測定小震震源參數的方法，該方法反演結果的可靠性在前人的對比論證中得到了證實（胡新亮等，2004;付虹等，2009），且已經在首都圈的應力場研究中得到應用（李瑞莎等，2008;張紅艷等，2009;武敏捷等，2011，2012）;第二部分為2009—2019年10月京津冀地區ML≥2.0地震的震源機制解，使用格點嘗試法（許忠淮等，1983）利用P波初動極性計算震源機制解，P波初動極性資料主要來自“地震編目系統”下載的震相報告和從事件波形中人工讀取的P波初動。為了保證震源機制解的可靠性，參與計算的P波初動極性數據一般要求10個以上，且在震中周圍均勻分布，矛盾比小于0.20。胡新亮等（2004）對首都圈數字地震臺網內的地震分別使用格點嘗試法和垂直向最大振幅比方法測定震源機制，二者結果基本一致，表明2種方法計算結果可以聯合起來進行京津冀地區地殼應力場特征分析。

1.2震源機制解分類方法

對震源機制解分類有多種方法，如使用P，T軸參數進行三角形分類（Frohlich，1992，2001）、根據震源機制解3個應力軸傾角的世界應力圖劃分標準（Zoback，1992）、以3個應力軸的傾角最大值判斷震源機制破裂類型（刁桂苓等，2011）等，以上都是以應力軸參數對震源機制解分類。本文采用Shearer 等（2006）提出的依據節面滑動角ri對震源機制解進行分類的方法。該方法的優勢在于提供單個標量斷層類型值用于表征震源錯動類型，標量斷層類型值fptype為：fptype=r90（1）式中：當r190ri，ri≤90（i=1，2）（2）標量斷層類型值從-1（正斷層）→0（走滑斷層）→1（逆斷層）變化。Lin 和Okubo（2016）取一個介于0～1的Xf值來區分斷層類型：-1≤fptype<-Xf，正斷層-Xf≤fptype≤Xf，走滑斷層Xf

1.3應力場反演方法

目前，使用震源機制解進行應力場反演應用最廣泛的2種方法為：震源機制應力反演（FMSI）的網格搜索方法（Gephart，Forsyth，1984;Gephart，1990）和自助線性應力反演（LSIB）方法（Michael，1984，1987）。Hardebeck和Hauksson（2001）使用含噪聲的合成數據對這2種方法進行測試，2種方法都能精確地確定應力方向;FMSI方法通常對于高質量數據更準確，而LSIB方法對于離散度大的數據更準確;LSIB方法產生的置信區間通常較為合理，而FMSI方法的置信區間通常過大。合成數據測試結果表明，只要滿足應力場均勻和斷層的滑動方向與剪應力方向一致的模型假設，由LSIB方法產生的置信區間也應適用于真實數據。Hardebeck和Michael（2006）引入阻尼最小二乘反演方法對LSIB方法加以改進，提出了區域尺度的應力張量阻尼反演方法SATSI（Spatial And Temporal Stress Inversion）算法;Martínez-Garzón等（2014）將SATSI算法移植到Matlab環境中，形成一款結合可靠經典方法的、新簡化用戶處理及可視化工具的應力反演Matlab軟件包。由于本文所用震源機制解數據量較大，且2種不同的方法計算得到的數據離散度較大，故選用MSATSI軟件包反演京津冀地區地殼應力場，以期可以準確地約束應力張量的方向（Hardebeck，Michael，2006）。

2震源機制解統計特征

本文參考Lin和Okubo（2016）的做法，對京津冀地區中小地震震源機制解的破裂類型比例進行定量估計，取Xf=0.33，將[-1，1]分成3個相等的部分，震源機制解也相應分為走滑斷層、正斷層和逆斷層（圖1）。對上述地震進行統計后發現，京津冀地區有1 315個走滑型震源機制解，占整體數量的60.1%，628個正斷型地震，占28.7%，244個逆斷型地震，占11.2%，該區震源錯動類型主要以走滑型和正斷型地震為主。不同類型的震源機制解空間分布（圖1）顯示，走滑型地震在整個區域彌散分布，而正斷型和逆斷型地震基本上沿著張渤地震帶和華北平原地震帶分布，沒有形成明顯的單一類型的叢集區。

將震源機制解按節面走向、傾角、滑動角和P軸方位角繪制成直方圖進行統計（圖2）。由于絕大多數地震震級太小，無法區分真正的斷層面和輔助斷層面，故將2 個節面等同看待，參數合并統計。節面走向大體呈各個方位均勻分布，在NNE—SSW和NWW—SEE向稍顯優勢（圖2a）;節面傾角以接近直立的高傾角最多，表明以走向滑動為主，幾乎沒有近水平的節面（圖2b）;節面滑動角以走向滑動為主，左旋、右旋走滑數量相仿，正傾滑明顯多于逆傾滑（圖2c）。P軸方位角在NEE—EW和SWW—EW向存在顯著優勢取向，可以看作是京津冀地區最大水平壓應力方向（圖2d）。上述震源機制解的統計特征與前人研究結果（李欽祖，1980;許忠淮等，1983;武敏捷等，2012;劉靜等，2016;劉麗等，2017）基本一致。

3地殼應力場反演

采用1°×1°的步長，將京津冀地區（36.0°～43.0°N，113°～120°E）劃分網格，使用MSATSI軟件包反演該地區精細的地殼應力場。該方法通過在反演過程中加上一個合理的阻尼，使反演的相鄰區域的應力張量變化最小，消除人為劃分反演區域導致的應力偏轉假象，更好地凸顯整體應力場特征（Hardebeck，Michael，2006）。在反演之前要先確定阻尼系數e，它控制著理論值與觀測數據之間擬合殘差和應力反演模型長度的相對權重，對反演結果起著至關重要的作用。圖3為模型長度—數據擬合殘差的折中曲線，其拐點在1.2，因此最佳阻尼系數取e=1.2。反演中每個網格采用的震源機制解至少為6個，隨機選取其中一個節面進行構造應力場反演。反演應力場參數的置信水平的可設定范圍為68%～95%，本文設定為95%。重采樣的迭代次數對于解的不確定度評價很重要，重采樣次數過少導致評估結果沒有統計意義，過多會導致計算量大，反演中重采樣次數的默認范圍是1 000～5 000次（Martínez-Garzón et al，2014），本文對所有數據進行2 000次抽樣估計。經過計算得到最優狀態下的3個應力主軸的方向和相對應力大小（R值）（Gephart，Forsyth，1984）為：R=（σ2-σ1）/（σ3-σ1）（4）式中：σ1，σ2，σ3分別表示最大、中間、最小主壓應力。

從表1和圖4可以看出，最優主壓應力軸σ1的優勢方向為NEE—EW向，且有大約半數的網格傾角接近直立，表明應力狀態為正斷層類型，而網格編號05，09，23的最優主壓應力軸σ1方向為NNW—NWW向，傾角都接近直立，表明最大主應力軸垂直，這與前人用中強震震源機制解和綜合震源機制解研究華北地區地殼應力場的結果有所差異（李欽祖，1980，1982;魏光興等，1982;許忠淮等，1983），主要表現為最優主應力軸的傾角直立或水平，這可能與本文所用的震源機制解多為小震、易受局部構造運動控制、且具有隨機性特征有關，反映了整體背景應力場下的非均勻性和復雜性;最優主張應力軸σ3的優勢方向為SSE—NWW向，且傾角均接近水平。

從圖4中最大主壓應力σ1軸的指示方向可以看出，在39°N以北的首都圈地區最大主壓應力軸自西向東存在由NE向近EW向的偏轉，這與前人對首都圈構造應力場的研究結果基本一致（張紅艷等，2009;武敏捷等，2012;黃驥超，萬永革，2015;樊文杰等，2019），而且本文的分區更加精細，應力場的變化連續可靠。由最大、中間、最小主壓應力軸最優解推斷的張渤地震帶陸地段的應力狀態表現為，自西向東存在一個正斷層—走滑斷層—正斷層的轉換過程，這與張紅艷等（2009）、武敏捷等（2012）、樊文杰等（2019）以走滑型為主的應力場特征存在一定差異;同時，在黃驥超和萬永革（2015）對首都圈的構造應力場反演結果中，北京—張家口—大同地區和唐山及鄰區都表現出少量正斷層性質。39°N以南地區現今構造應力場保持穩定，最優主壓應力軸σ1為NEE—SWW向，大部分網格應力狀態顯示走滑型。

相對應力大小R值是中間主應力σ2軸的值，是更接近于最大主壓應力軸σ1還是最小主壓應力軸σ3的大小的度量。當R值趨于中間值0.5時，表示3個應力主軸方位可完全分辨;而當R值趨于取值端點0或者1時，表示壓縮或拉張應力狀態完全不能分辨，只能確定它在垂直于拉張或壓縮應力主軸的平面內。一般認為當R值小于0.5時，獲得的應力狀態偏拉張性質，反之則為偏壓縮性質（Guiraud et al，1989）。從表1和圖5可以看出，所有網格的R值均小于0.5，表明京津冀地區的應力狀態呈現偏拉張性質。同時從圖4可知，最大主應力軸σ1的不確定范圍最大，最大主壓應力軸σ1和中間應力軸σ2表現的壓應力狀態是一致的，兩軸無法區分（Guiraud et al，1989;萬永革等，2011）;而最小主壓應力軸σ3的不確定度變化范圍相對穩定，我們可以推斷京津冀地區現今地殼應力場處于一個相對統一的NNW—SSE向的拉張作用控制下。4結論和討論

基于京津冀地區及鄰區2 187個中小地震震源機制解，通過計算標量斷層類型值將其劃分為3種震源錯動類型，采用MSATSI軟件包反演該地區1°×1°網格的精細地殼應力場，研究結果表明：研究區震源機制類型以走滑斷層和正斷層為主，占比約90%;不同類型的震源機制解沒有形成明顯的單一類型的叢集區，反映出小地震的發生除了受背景應力場的控制，還受局部構造條件控制且具有一定的隨機性。這為我們利用小地震的震源機制解來推斷構造應力主軸的方向提供了更大的優越性，震源機制解的一對節面解中既含有構造應力場的信息，也可能含有局部介質不均勻結構的影響（許忠淮等，1983）。

地殼應力場反演結果顯示，所有網格的相對應力大小R值均小于0.5，表明京津冀地區應力狀態偏拉張性質，而且最小主壓應力軸σ3的不確定度變化范圍相對穩定，表明現今京津冀地區地殼應力場處于一個相對統一的NNW—SSE向的拉張作用控制下。這一結果與前人的研究結果比較一致，高名修（1979）認為華北地區處于引張應力場環境，有限元模擬顯示主張應力大約為主壓應力2～8倍（劉峽等，2010）。利用反演結果中最大、中間、最小主壓應力軸最優解推斷的應力狀態顯示，39°N以北地區的最大主壓應力軸方向存在一定角度的偏轉，同時應力狀態由西向東存在一個正斷層—走滑斷層—正斷層的轉換過程;而39°N以南地區的現今構造應力場保持穩定，大部分網格應力狀態顯示走滑型。

39°N以北地區與前人研究中的張渤帶陸地段和首都圈范圍基本一致，大致分為3段：西段、中段和東段。西段為張家口—大同地區，應力狀態表現為正斷層類型，該區為山西斷陷帶右旋走滑在NE端形成的一個張性構造區（徐錫偉等，2002），活動構造、遙感和垂直形變研究結果都表明山西斷陷帶中北段至張家口地區總體以張性正斷層活動為主（孫啟凱等，2018;吳玉濤等，2018;申星等，2019;彭遠黔等，2019），該區的中小地震震源機制解以正斷型和走滑型為主（孫小入等，2019）。中段為北京—天津地區，應力狀態為走滑類型和正斷層，該區歷史上曾發生過1679年三河—平谷地震，使用現今小震資料推斷該地震的錯動類型為走滑兼少量正斷（王曉山，2017），2018年發生在該區的河北永清地震的震源機制為走滑兼正斷類型（王曉山等，2018;李赫等，2020），GPS觀測資料亦顯示該區存在高剪切變形區（陳長云，2016）。東段為唐山地區，應力狀態為走滑類型和正斷層，楊雅瓊等（2016）對唐山地震序列的現今小震震源機制分段擬合最佳應力張量，得到的結果為正斷層和走滑類型兼有。綜合前人多種資料的論證表明，本文的研究結果與前人結果具有較好的一致性，但仍存在少許差異，可能由所使用的反演方法、網格劃分和數據量等導致，有待進一步研究分析。

GPS主應變率結果顯示：京津冀大部分地區地殼變形特征主要受NEE的主壓應變控制（陳長云，2016），京津冀地區的剪切波分裂的快波偏振方向為NEE—EW向（吳晶等，2008a，b;趙博等，2011;吳鵬等，2017）。該結果與該區使用中小地震震源機制解反演的構造應力場最優解基本一致，同時也和前人使用中強地震震源機制解和綜合震源機制解得到的華北地區一致性良好的應力場吻合，表明中小地震震源機制的結果也可以較好地用于區域構造應力場的研究。京津冀地區是中國大陸較早實現數字化地震觀測網絡的地區，目前積累了豐富的中小地震波形資料，開展中小地震震源機制解及構造應力場的動態演化研究可以為該區地震動力學研究提供具有物理意義的信息，更好地服務于該地區的防震減災綜合研究。

參考文獻：

陳長云.2016.張家口—渤海斷裂帶分段運動變形特征分析[J].地震，36（1）：1-11.

刁桂苓，徐錫偉，陳于高，等.2011.汶川MW7.9和集集MW7.6地震前應力場轉換現象及其可能的前兆意義[J].地球物理學報，54（1）：128-136，doi：10.3969/j.issn.0001-5733.2011.01.014.

樊文杰，崔效鋒，胡幸平，等.2019.首都圈地區震源機制解及現今構造應力場時空變化特征研究[J].地震學報，41（1）：33-45，doi：10.11939/jass.20180064.

付虹，王紹晉，龍曉帆，等.2009.利用數字化速度波形和仿真位移波形資料求解震源機制解的比較[J].地震研究，32（1）：18-24.

高名修.1979.華北斷塊構造區的現代引張應力場[J].地震地質，1（2）：1-12.

郭良遷，馬青，杜雪松，等.2008.華北地區斷層形變與地震的關系[J].大地測量與地球動力學，28（3）：14-20.

胡新亮，刁桂苓，馬瑾，等.2004.利用數字地震記錄的P，S振幅比資料測定小震震源機制解的可靠性分析[J].地震地質，26（2）：347-354.

Gephart J W，Forsyth D W.1984.An improved method for determining the regional stress tensor using earthquake focal mechanism data：application to the San Fernando Earthquake sequence[J].J Geophys Res，89（B11）：9305-9320.

Gephart J W.1990.FMSI：A FORTRAN program for inverting fault/slikenside and earthquake focal mechanism data to obtain the regional stress tensor[J].Computers and Geosciences，16（7）：953-989.

Guiraud M，Laborde O，Philip H.1989.Characterization of various types of deformation and their corresponding deviatoric stress tensors using microfault analysis[J].Tectonophysics，170（3-4）：289-316.

Hardebeck J L，Hauksson E.2001.Stress orientations obtained from earthquake focal mechanisms：what are appropriate uncertainty estimates？[J].Bull Seis Soc Amer，91（2）：250-262.

Hardebeck L L，Michael A J.2006.Damped regional-scale stress inversions：methodology and examples for southern california and the coalinga aftershock sequence[J].J Geophys Res，111（B11）：B11310，doi：10.1029/2005JB004144.

Lin G，Okubo P G.2016.A large refined catalog of earthquake relocations and focal mechanisms for the island of Hawaii and Its seismotectonic implications[J].J Geophys Res，121（7）：5031-5048，doi：10.1002/2016JB013042.

Martínez-Garzón P，Kwiatek G，Ickrath M，et al.2014.MSATSI：A MATLAB package for stress inversion combining solid classic methodology，a new simplified user-handling，and a visualization tool[J].Seismol Res Lett，85（4）：896-904.

Maury J，Cornet F H，Dorbath L.2013.A review of methods for determining stress fields from earthquakes focal mechanisms;application to the sierentz 1980 seismic crisis（upper Rhine graben）[J].Bull Soc geol，184（4-5）：319-334.

Michael A J.1984.Determination of stress from slip data：faults and Folds[J].J Geophys Res，89（B13）：11517-11526.

Michael A J.1987.Use of focal mechanisms to determine stress：a control study[J].J Geophys Res，92（B1）：357-368.

Shearer P M，Prieto G A，Hauksson E.2006.Comprehensive analysis of earthquake source spectra in Southern California[J].J Geophys Res，111（B6）：303，doi：10.1029/2005JB003979.

Shen Z K，Jackson D D，Ge B X.1996.Crustal deformation across and beyond the Los Angeles Basin from geodetic measurements[J].J Geophys Res，101（B12）：27957-27980.

Zoback M L.1992.First- and secondorder patterns of stress in the lithosphere：the world stress map project[J].J Geophys Res，97（B8）：11703-11728，doi：10.1029/92JB00132.

Characteristics of Crustal Stress Field in Beijing-Tianjin-Hebei Region