山東地區地殼三維縱橫波速度和泊松比結構及其對地震活動性的影響＊

范建柯丁志峰郭慧麗蘇道磊張斌，

1） 中國北京 100081 中國地震局地球物理研究所

2） 中國山東青島 266071 中國科學院海洋研究所

3） 中國山東青島 266071 中國科學院海洋地質與環境重點實驗室

4） 中國濟南 250001 濟南市地震監測中心

引言

山東地區位于中國大陸東部，自中生代以來經歷了復雜的地質構造運動，如中生代早期華北地塊與揚子地塊的碰撞、中生代中晚期的華北克拉通破壞以及新生代以來的拉張變形等.陸地區域內發育沂沭斷裂帶（郯廬斷裂帶山東段）和聊考斷裂帶兩大斷裂系統，北部海域發育張家口——蓬萊斷裂帶，導致地震活動頻繁，并伴有大地震的發生，如1668年郯城M8.5地震，造成了重大的人員傷亡和經濟損失.

2003年6月5日和2020年2月18日，山東省最大的兩個城市青島和濟南分別發生了一次M4.1地震（圖1），雖然未導致較大的人員傷亡和經濟損失，但由于其發生于人口均近千萬的兩大城市，且有感范圍較大，余震次數較多，而上述兩個地區通常被認為屬于“弱震”區，因此，這兩次地震仍然引發了非常大的關注.前人綜合余震重定位結果、震源機制解等數據，推測這兩次地震可能屬于相對完整巖體條件下的一次新破裂活動（潘元生等，2004）或區域構造應力作用下附近斷裂或次級派生斷裂活動的結果（張斌等，2020）.另外，自2000年以來，山東半島地區相繼發生了嶗山震群、乳山震群和長島震群，但目前上述發震區尚未有高精度的層析成像數據，速度結構與發震機理仍不清楚.因此，有必要查明上述地區的淺部及深部速度結構形態，探討速度結構與地震發震機理之間的關系，為后期的防震減災提供理論支持.

圖1 山東及周邊地區構造分區圖（修改自蘇道磊等，2016）Fig.1 Tectonic settings in and around Shandong area（modified from Su et al，2016）

1 數據與方法

本文所用數據來源于山東及鄰區96個臺站（圖2a）記錄的區域內地震事件的P波和S波到時數據，主要分為兩部分：第一部分為1975年1月至2014年1月發生的天然地震事件（蘇道磊等，2016）；第二部分為2016年1月至2019年12月發生的天然地震事件.原始數據共包含7 271 個地震事件.

圖2 臺站（a）及地震（b）分布圖Fig.2 Distribution of seismic stations （a） and local earthquakes （b） used in this study

為保證反演結果的準確性，需要對地震事件進行嚴格挑選，設定標準如下：① 每個地震事件至少被4個臺站接收；② 震相走時殘差絕對值小于3.0 s；③ 重定位前后，發震時刻偏差小于2.0 s，水平向偏差小于6 km，震源深度偏差小于8 km.最終篩選出4 652個地震事件，包括3萬6 482個P波震相到時和3萬2 600個S波震相到時用于成像反演（圖2b）.

本文采用TOMOG3D三維層析成像反演方法（Zhaoet al，1992）.該方法采用三維網格節點表示空間的速度分布，空間內每一點的速度值由周圍八個節點的速度值進行線性插值獲得，允許三維空間內存在間斷面，通過偽彎曲算法快速準確地計算射線路徑和走時，并能夠同時處理近遠震及后續震相到時數據.

合適的一維初始速度模型對最終反演結果的準確性至關重要.蘇道磊等（2016）分別測試了三個初始速度模型：① 魯西地區地殼速度模型，同時考慮莫霍面起伏（嘉世旭，張先康，2005）；② 山東地震臺網定位用地殼速度模型；③ 華北地區地殼速度模型（陳立華，宋仲和，1990）.測試結果表明，魯西地區地殼速度模型的走時殘差均方根最小，也更符合真實地層情況，而華北地區地殼速度模型走時殘差遠大于其它兩個速度模型.因此，本文僅對比前兩個初始速度模型的走時殘差，并同時考慮莫霍面起伏的影響.其中，莫霍面埋深數據來自于CRUST1.0模型（Laskeet al，2013）和鄭宏等（2021）利用接收函數反演得到的山東地區莫霍面埋深等（圖3）.初始S波速度模型由P波速度除以1.732得到.

圖3 基于不同模型反演所得的莫霍面埋深（a，b）及本研究所用莫霍面埋深（c）Fig.3 Moho depths inverted from different models （a，b） and used in this study （c）

經過計算，上述兩個初始速度模型的殘差均方根分別為0.841 s和0.860 s，走時殘差分布如圖4所示.對比發現，魯西地區速度模型要優于山東地震臺網定位用地殼速度模型，與蘇道磊等（2016）的結論相似.因此，后續的層析成像反演采用魯西地區速度模型，并考慮莫霍面起伏的影響.

圖4 使用魯西地區地殼速度模型（a）和山東地震臺網定位用地殼速度模型（b）在地震重定位前后所得的走時殘差分布圖Fig.4 Travel time residual distribution based on crustal velocity model in Luxi area （a） and Shandong seismic network （b） before and after relocation

本研究中，初始模型三維網格節點橫向間隔為0.2°×0.2°，縱向上分別在1，10，20和30 km深度處設置節點層.反演前，基于初始速度模型和原始到時數據對所有地震事件進行了重新定位.重定位前后的走時殘差分布如圖4所示，重定位后的總體走時殘差均方根由0.841 s降低到 0.649 s，其中 P 波走時殘差均方根由 0.792 s降低到 0.620 s，S 波走時殘差均方根由0.893 s降低到0.680 s，說明地震定位精度有了較大幅度的提高.走時殘差絕大多數位于±2 s以內，因此選擇走時殘差在±2 s以內的震相參與最終的成像反演.反演采用帶阻尼和平滑因子的最小二乘（least squares QR-factorization，縮寫為 LSQR）方法（Paige，Saunders，1982）得到最終的三維速度結構.通過大量的測試，P波和S波速度反演采用的最優阻尼和平滑因子均為 5.0 和 50.0 （圖5）.

圖5 P 波（a，b）和 S 波（c，d）成像中不同阻尼和平滑因子對應的三維速度模型標準差與走時殘差均方根關系曲線Fig.5 Trade-off curves between the standard deviation of the 3-D velocity model and the root-mean-square of travel time residual according to the damping （a and c） and smoothing （b and d） parameters for P- （a，b） and S-wave （c，d） tomographies

2 分辨率測試

在分析成像結果之前，需要對成像結果的可靠性進行評估，一般采用檢測板（checkerboard）方法（Zhaoet al，1992）.在該方法中，首先將三維空間內相鄰網格節點分別設置±3%的速度擾動，然后利用相同的地震和臺站分布計算理論走時，并在計算理論走時過程中加入標準差為0.1 s的隨機誤差來檢測計算穩定性.隨后對得到的理論走時在初始一維速度模型的基礎上進行反演，通過對比三維網格節點處擾動值的恢復情況對成像結果的可靠性進行評估.如果反演后擾動值與反演前擾動值的分布相似，則說明分辨率較好.

圖6顯示了橫向間隔分別為 0.33°×0.33°，0.4°×0.4°和 0.5°×0.5°的分辨率測試結果.結果顯示，對于P波和S波成像結果，研究區內大部分地區的分辨率能達到0.4°×0.4°，部分地區如魯西南地區和山東半島地區可達到0.33°×0.33°的分辨率.其中：網格間距為0.33°×0.33°時，P波和S波檢測板擾動值恢復準確性分別達到74%和81%，振幅恢復達到70%的比例約為47%和59%，振幅恢復達到100%的比例約為32%和41%；網格間距為0.4°×0.4°時，P波和S波檢測板擾動值恢復準確性分別達到76%和82%，振幅恢復達到70%的比例約為52%和62%，振幅恢復達到100%的比例約為36%和41%；網格間距為0.5°×0.5°時，P波和S波檢測板擾動值恢復準確性分別達到77%和82%，振幅恢復達到70%的比例約為55%和64%，振幅恢復達到100%的比例約為37%和42%.

圖6 檢測板分辨率測試結果Fig.6 Results of checkerboard resolution tests

3 結果與討論

經過反演，P波和S波走時殘差均方根由反演前的0.620 s和0.680 s分別降低到0.378 s和0.417 s，減少近40%.利用P波、S波反演結果計算得到了泊松比成像結果，如圖7所示.蘇道磊等（2016）的P波層析成像結果揭示出：在1——10 km深度切片上，沂沭斷裂帶沿構造走向表現出強烈的橫向不均一性，高低速異常交替出現，膠東半島北部海域、膠萊盆地和濟陽坳陷主要表現為低速異常；20——30 km深度，魯西地區存在較大范圍的低速異常（蘇道磊等，2016）.上述結果與本文的P波成像結果非常一致（圖7a）.另外，本文的S波成像結果與Li等（2018）基本一致，如在10 km深度上，沂沭斷裂帶北部以低速異常為主（圖7b），20——30 km山東半島表現為大范圍的低速異常（圖7a）.但由于本文應用了更多的P波和S波到時數據，因此本文的成像結果具有更高的分辨率，如在1 km和10 km深度，沂沭斷裂帶南部西側的兩個低速異常區相互分離，斷裂帶東側表現為高速異常（圖7a），這在蘇道磊等（2016）的研究中揭示得并不明顯.同時，本文的檢測板測試結果也要優于蘇道磊等（2016），魯西南和山東半島地區分辨率能夠達到0.33°×0.33°.山東半島地區中下地殼存在大范圍低速異常，說明該地區可能存在較強烈的地幔上涌（李志偉等，2006）.本文同時利用Liu和Zhao （2018）提出的方法計算了1——30 km的P波、S波和泊松比的平均值.平均泊松比異常與鄭宏等（2021）利用接收函數反演得到的山東地區的泊松比分布基本一致，進一步說明了本文成像結果的可靠性.

圖7 山東地區不同深度 P 波速度擾動（a）、S 波速度擾動（b）及泊松比相對擾動（c）結果底部3幅圖分別為1——30 km的P波、S波速度擾動及泊松比相對擾動的平均值Fig.7 Map views of P-wave （a），S-wave （b） and Poisson’s ratio （c） perturbations in Shandong area The averages of P-wave，S-wave and Poisson’s ratio anomalies from 1 to 30 km are shown at the bottom

泰山作為山東地區的最高點，自新生代以來經歷了多期快速抬升（李理，鐘大賚，2006）.接收函數研究結果顯示，泰山地區具有較薄的地殼厚度（32 km左右）和較大的泊松比（0.27），且沒有明顯的方位各向異性（鄭宏等，2021）.本文的成像結果顯示，泰山地區下方存在明顯的低速異常（圖8左下四幅圖），這說明該地區仍存在較強的地幔上涌，導致泰山新生代以來出現顯著抬升，并且現在仍處于抬升階段（鄭宏等，2021）.另外，在淺部1——10 km深度（圖8左上四幅圖），泰山北部為高速異常，南部為低速異常，這與泰山的巖性和地質構造是一致的：以泰山山前斷裂為界，斷裂北側為泰山主體，主要出露前寒武紀泰山變質雜巖，以隆升為主；山前斷裂南側以沉降為主，主要為泰安——萊蕪盆地巨厚的第三系碎屑巖及第四系沉積物（李理，鐘大賚，2006）.

圖8 長清地震及泰山周邊不同深度處 P 波速度擾動（a）、S 波速度擾動（b）及泊松比相對擾動（c）結果Fig.8 Map views of P-wave （a），S-wave （b） and Poisson’s ratio （c） perturbations around Changqing earthquake and Mount Tai

2020年2月18日，濟南長清發生了一次MS4.1地震，震源深度約為2.7 km，并引發近40次的余震活動（張斌等，2020）.本文的成像結果顯示，該地震震中位于P波、S波高低速異常和高低泊松比異常過渡帶，可能與該地震有關的長清斷裂也位于P波、S波和泊松比高低異常過渡帶（圖8上半部和圖9）.震中東部的高速異常與該地區的背景噪聲成像結果一致，可能與濟南侵入巖體有關（雷霆，2020）.震源機制解顯示，該地震具有正斷兼走滑性質，其余震多為走滑性質（張斌等，2020）.另外，余震重定位結果顯示，主震兩側的破裂呈明顯不對稱分布，以西北側破裂為主（張斌等，2020）.因此，濟南長清MS4.1地震可能是區域構造應力下長清斷裂發生左旋走滑運動的結果（張斌等，2020），這與本文的層析成像結果相一致.另外，主震的西北向不對稱破裂的形成可能是由于震中東側高速異常體的存在阻礙了地震的東向破裂.

圖9 過長清地震的剖面 AA′（a）的 P 波速度擾動（b）、S 波速度擾動（c）和泊松比擾動（d）Fig.9 Vertical cross section of P-wave （b），S-wave （c） and Poisson’s ratio （d） perturbations along a profile AA′ （a） across the Changqing earthquake shown in the inset map

2003年6月5日青島嶗山地區發生了ML4.1地震（圖10a左側），2004年11月1日附近地區又發生一次ML3.6地震，這兩次地震都形成了震群序列.根據定位結果，這兩次震群序列大致呈北西走向，與附近的主要斷裂近似垂直、與次級斷裂展布方向基本一致但并不重合，可能屬于相對完整巖體條件下的一次新的破裂（潘元生等，2003，2005）.根據本文的成像結果，這兩次地震及震群序列發生在P波高低速異常過渡帶、S波低速異常區和高低泊松比過渡區且偏向于高泊松比區域（圖10a——c 1km），同時深部存在明顯的P波和S波低速異常（圖10a，b 10——30 km），這與1995年日本神戶地震（ML7.2）震源區的速度和泊松比結構非常相似（Zhaoet al，1996）.另外，震中附近水資源豐富，同時為花崗巖侵入形成地區，且構造裂隙發育（趙廣濤等，1996）.因此本文認為，可能是深部的流體填充相對完整巖體內的裂隙并引發破裂，從而導致地震的發生.需要注意的是，雖然這兩次地震及震群序列與本地區北東走向的主斷層沒有直接關系，但本文的成像結果顯示，該地區北東走向的主斷層基本位于低速區內以及高低泊松比異常的過渡帶（圖10），因此仍需警惕該地區發生震級更大的中強地震的可能性.

乳山地區位于大別——蘇魯超高壓變質帶的東部，可能存在大量的中生代巖漿侵入體（郭敬輝等，2005）.自2013年10月開始，記錄到的乳山震群地震數已超過1萬余次，最大震級為M5.0，目前震群活動可能已結束.震群周邊的斷裂大多以NE-NNE向為主，而近期的重定位結果和震源機制解顯示，乳山地震序列主要呈NW向展布，絕大多數地震發生于中上地殼，發震斷層可能為傾角近直立的左旋走滑斷層，與區域內的主要斷裂并不一致（張斌等，2017）.本文的成像結果顯示，乳山震群周邊1 km深度以高速異常和高泊松比異常為主，10 km深度主要為高低速異常過渡區，20——30 km深度以低速異常和低泊松比異常為主（圖10）.另外，震群附近地熱資源較為豐富（田粟，2012）.因此，本文認為，乳山地區深部存在熱地幔物質上涌，內部含有的流體注入淺部相對較完整的侵入體內或侵入體之間，在區域應力場的作用下導致侵入體的破裂或侵入體之間隱伏斷裂的活動，從而引發乳山震群.這與前人對乳山震群進行的震源譜參數反演、震中空間分布、震中隨時間的演化規律以及震源區應力狀態的研究結果是一致的（鄭建常等，2016；王鵬，2019）.

圖10 山東半島地區不同深度 P 波速度擾動（a）、S 波速度擾動（b）及泊松比相對擾動（c）結果Fig.10 Map views of P-wave （a），S-wave （b） and Poisson’s ratio （c） perturbations in and around Shandong Peninsula

長島震群的發生可能主要受控于NWW向的張家口——蓬萊斷裂，這是一條深大斷裂，可能已切穿莫霍面甚至巖石圈，成為地幔熱物質或基性物質上涌的通道（張嶺等，2007），斷裂主要以正斷兼走滑運動為主（索艷慧等，2013）.歷史上，附近海域曾發生過M6.0和M7.0左右的大地震（王志才等，2006）.本文的成像結果顯示（圖10），張家口——蓬萊斷裂帶的地殼速度結構特征與沂沭斷裂帶具有非常大的相似性（蘇道磊等，2016），地殼速度和泊松比結構在1——10 km深度處具有非常強烈的橫向不均一性，斷裂一側高低速異常交替分布，斷裂帶整體位于高低速異常與高低泊松比異常的過渡帶，20 km以下以低速和低泊松比異常為主，但部分地區仍有顯著的高泊松比異常，可能反映了地幔熱物質的上涌.速度和泊松比異常在10 km左右發生明顯變化，這與前人得到的b值在8.5 km左右發生轉折是一致的（申金超等，2019）.本文認為，深部地幔熱物質沿斷裂帶上涌所產生的構造應力，和/或地幔熱物質內部含有的流體沿斷裂上涌或侵入裂隙導致了長島震群、甚至周邊地區的強震活動.

4 結論

利用山東及鄰區的地震臺站記錄的P波和S波到時數據反演獲得了研究區內高精度的縱橫波速度結構和泊松比異常分布形態.研究結果揭示了研究區內的地殼結構具有強烈的橫向不均一性.2020年濟南長清MS4.1地震可能是區域構造應力下長清斷裂發生左旋走滑運動的結果，震中東側的高速異常體可能阻礙了地震的東向破裂.2003年青島嶗山ML4.1地震可能是由深部的流體填充相對完整巖體內的裂隙并引發破裂所致.嶗山震群、乳山震群和長島震群可能與深部流體有非常強的相關性.