Detecting shallow sedimentary structures in Sanhe and Tangshan seismic regions using H/V spectral ratio method
-
摘要:
利用一条穿过三河和唐山地震区的密集地震探测剖面数据，采用H/V谱比法获得了三河、唐山地震区及邻近区域的地壳浅部沉积结构特征、场地共振频率及易破坏程度等参数。结果表明：研究区松散沉积层厚度约为100—800 m，呈现出NW浅而SE深的变化特征；探测剖面西北部通县隆起的沉积层厚度约为350—450 m，沉积层界面起伏平缓；大厂凹陷沉积层厚度约为300—600 m，横向变化特征显著。测线东南部唐山地震区地壳浅部存在上下两组较明显的沉积界面：上层界面深度约为100 m，呈水平展布；下层界面深度约为300—800 m，且向东南方向逐渐加深。三河和唐山地震区的场地放大系数约为3—4，场地易破坏程度均大于20，显示出地表建筑物易破坏程度较高的特征。
Abstract:In the Sanhe and Tangshan seismic regions and their adjacent areas, the shallow crustal sedimentary structures detecting using H/V spectral ratio method has important scientific significance for regional seismic hazard risk assessment and engineering seismic fortification. In this paper, according to the data of a dense seismic sounding profile passing through the Sanhe and Tangshan seismic regions, the characteristics of the shallow crustal sedimentary structure, the site resonance frequency and the degree of vulnerability to damage in the Sanhe and Tangshan seismic regions and adjacent areas are obtained by using the H/V spectral ratio method. The results show that the thickness of loose sediment layer in the studied area is about 100–800 m, showing the characteristics of shallow in the northwest and deep in the southeast. The sedimentary thickness of the Tongxian uplift in the northwest of the exploration profile is about 350–450 m, the sedimentary interface fluctuates gently, and the sedimentary thickness of Dachang depression is about 300–600 m, with significant lateral variation characteristics. There are two sets of obvious sedimentary interfaces in the shallow part of the Tangshan earthquake region in the southeastern part of the survey line. The depth of the upper layer interface is about 100 m and is horizontally distributed. The depth of the lower layer interface is about 300–800 m, and gradually deepens toward the southeast. The site magnification factor of the Sanhe and Tangshan earthquake regions is about 3–4, and the site’s vulnerability to damage is greater than 20, showing that the surface buildings are more vulnerable to damage.
-
引言
在地震破裂过程研究中，破裂模型反演结果的可靠性依赖于断层参数的准确性。研究人员一般通过反演震源机制解获得地震断层参数。在震源机制反演研究方面已积累了大量方法和技术，例如P波初动、P/S振幅比（梁尚鸿等，1984）、面波波形反演（Patton，1980）、体波波形反演（Wallace，Helmberger，1982）、近震全波场反演（Dreger，Helmberger，1993）、CAP （cut and paste）方法（Zhao，Helmberger，1994）和W震相方法（Kanamori，1993）等。对于地震破裂模型研究，则一般采用远场地震波、近场强震地震波、全球定位系统（global positioning system，缩写为GPS）、合成孔径雷达干涉 （interferometric synthetic aperture radar，缩写为InSAR）等数据进行反演。尽管这些反演工作在大地震研究中已经比较成熟，但在中等强度地震的应用中仍然存在一些困难。由于震级相对较小，中等地震的地震波和形变资料信号偏弱，信噪比不高，可能会对破裂模型反演结果造成影响。另一方面，对我国而言，中等强度的地震足以造成破坏，导致人员伤亡和财产损失，同样需要对震源特征开展细致的研究工作。
2017年精河MS6.6地震即是一次典型的中等强度地震。据中国地震台网中心（China Earthquake Networks Center，缩写为CENC）测定，北京时间2017年8月9日上午7时27分52秒，在新疆精河县发生MS6.6地震，震中位置为（44.27°N，82.89°E），震源深度为11 km。截至8月20日23时59分，共记录到324次余震，其中最大震级MS4.7。地震导致36人受伤，未造成人员死亡（阿里木江·亚力昆等，2017；常想德等，2017）。此次地震发生于天山北脉博罗科努山北沿的库松木楔克山前断裂的东段附近。天山南北两侧的构造主要受印度板块与欧亚板块的挤压控制。根据GPS测定结果，天山南北向的缩短速率约为20 mm/a （Abdrakhmatov et al，1996），地壳持续的缩短与增厚使天山不断向南北两侧盆地推覆扩展。库松木楔克山前断裂发育于天山推覆体前缘的逆断层（陈建波等，2007）。根据野外考察结果，本次地震未观测到明显的地表破裂。对本次地震序列重定位的研究（白兰淑等，2017；姜祥华等，2017；刘建明等，2017；徐志国等，2019；翟亮等，2019）认为，库松木楔克山前断裂可能为本次地震的发震构造（图1）。何骁慧等（2020）结合区域地质资料、卫星影像等，也认为本次地震的发震断层为位于库松木楔克山前断裂以北的精河南断层。精河地震附近区域地震活动性较强，自1900年以来，在距本次地震震中200 km范围内，一共发生过6次M≥6.0地震。在精河县境内，最近发生的中等强度地震为2018年10月16日的MS5.4地震（44.19°N，82.53°E）。
![]()
图 1 2017年8月9日精河MS6.6地震及其附近地区的构造示意图Figure 1. Geological background of the Jinghe MS6.6 earthquake on August 9，2017对于精河地震，数据质量可能影响了其破裂过程的反演结果。由于地震震级不大，断层尺度较小，远震地震波可能难以分辨较小空间尺度内的破裂过程；同时，近场强震台站方位覆盖较为有限，且区域沉积层较厚，较强的近地表或场地效应难以在反演中较好地考虑；InSAR数据记录的地表位移较小，造成观测位移的相对误差较大等。已有的研究中，不同破裂模型也存在一定差异。张勇等（2017）使用远场地震波形资料反演得到了本次地震的破裂时空过程，单新建等（2017）、刘传金等（2018）、施贺青等（2019）和Gong等（2019）使用InSAR同震形变数据对本次地震发震断层的滑动分布进行了反演；Zhang等（2020）分别使用远场地震波形反投影和远震、近场强震、InSAR数据联合反演得到了本次地震的破裂特征。在InSAR反演结果中，刘传金等（2018）、Gong等（2019）的破裂区域相对更浅；张勇等（2017）的远震反演结果，其破裂区域距离震源较近；Zhang等（2020）的联合反演结果中破裂则更倾向于向深处扩展。为确证此次地震的更细节的破裂行为，本文发展了一套确定断层面参数的流程，采用多种类型数据联合反演本次地震的破裂模型，以期获得更为详细的、准确的震源破裂细节。
精河MS6.6地震后，一些机构发布了其测定的震源位置和矩张量解（表1）。不同结果中的震中位置较为接近，震源深度位于11—24 km不等。根据本次地震发震构造研究结果（徐志国等，2019；何骁慧等，2020），真实的发震断层面应为近南倾的节面I，然而不同结构给出的南倾节面参数在一定程度上存在差异，无从确定哪个结果更为准确、更适合于破裂过程联合反演。由于不准确的断层参数可能对反演结果造成较大影响，因此本文试图通过密集分布的InSAR同震形变数据来确定更为准确的断层几何模型。
表 1 研究机构发布的2017年精河MS6.6地震震源定位结果和矩张量解Table 1. Epicenter locations and focal mechanism solutions of the 2017 Jinghe MS6.6 earthquake released by different research institutes研究机构 北纬/° 东经/° 震源深度/km 节面Ⅰ 节面Ⅱ 走向/° 倾角/° 滑动角/° 走向/° 倾角/° 滑动角/° CENC 44.270 82.890 11 76 44 80 269 47 99 USGS 44.302 82.832 20 92 60 92 269 30 87 GFZ 44.330 82.840 24 85 47 81 277 44 99 注：GFZ为德国地学中心German Research Centre for Geosciences的缩写. 1. 数据与方法
我们从美国地震学研究联合会（Incorporated Research Institutions for Seismology，缩写为IRIS）数据中心下载了全球地震台网震中距介于40°—90°范围内的远震宽频带垂直向P波记录。为保证良好的方位角覆盖，本文以最小方位角间隔5°对台站进行进一步挑选，最终挑选出27个台站的地震记录（图1左上插图），之后对这些地震记录以0.01—0.4 Hz的频带进行滤波，并将波形重采样为4 sps。
从中国地震局工程力学研究所国家强震动台网中心获取了32个近场强震台站的三分量地震波数据。为了尽可能减弱浅层速度结构对地震波的放大效应，我们使用较低的0.04—0.1 Hz频带对其进行滤波，同样将波形重采样至4 sps。经过多次试反演，最终从32个台站中选出11个波形拟合较好的台站，用于最终的联合反演（图1左下插图）。
利用升降轨道哨兵-1 （Sentinel-1）雷达数据，两个同震干涉像对（Track85和Track63）被处理用于约束断层几何参数和滑动分布（表2）。采用Feng等（2016）发展的自动处理程序gInSAR处理和分析InSAR数据。数据处理中，采用30 m分辨率的SRTM （Shuttle Radar Topography Mission）地形数据用于地形相位的去除和数据配准。为减轻反演迭代中的计算压力，采用四叉树算法对原始干涉图进行降分辨率重采样，将Track85和Track63的采样点分别降为3 285个和4 271个（表2）。
表 2 本文所用雷达数据信息Table 2. Information of the radar data used in this study轨道信息 轨道方向 时间信息 空间基线/m 时间基线/d 采样点个数 震前图 震后图 T85 升轨 2017-08-08 2017-08-14 − 92.2 37 3 285 T63 降轨 2017-08-07 2017-08-13 −119.6 37 4 271 对于地震波记录，基于CRUST1.0模型（Laske et al，2013），采用Wang （1999）的Qseis程序计算格林函数。对于InSAR记录，采用基于弹性位错理论的Okada （1985）方法计算地表形变。
2. 断层模型
有限断层破裂过程反演中，断层面为通过震源（点源）的平面断层。需将其沿走向和倾向划分为若干个正方形的子断层，震源位于其中一个子断层的中心。考虑到地震发生在隐伏断层上，断层几何参数不能很好地确定。为获得合理的断层参数，基于InSAR同震形变数据设计了一套确定断层面走向、倾角和震源深度的格点搜索流程。基于InSAR形变数据，采用不同的断层几何参数和位置参数进行静态滑动分布反演，根据残差最小的反演结果所对应的走向、倾角和震源位置，构建更合理的断层几何模型，并得到震源在断层面上的相对位置。
在以上搜索流程中，我们将对断层面位置参数的搜索转化为对震源深度的搜索。首先，设置一个较大且能较好地覆盖潜在滑动区域的断层面；其次，将断层面的上边界固定于地表，断层面沿倾向方向的移动等同于改变震源所在子断层在倾向方向上的序数，这样震源深度即可由震源到断层面上边界的距离乘以断层面倾角的正弦值得到。通过这种做法，只需对震源所在子断层在倾向方向上的序数进行搜索，即实现了对震源深度的搜索，同时也完成了对断层面沿倾向方向位置的搜索。
我们将三维搜索空间中的任意一个格点（走向、倾角、序数）记为FP （fault-plane parameters）。理想情况下，存在一个全局最优点，具有最小残差，三维搜索空间中的其它点的残差均随着最小残差点的距离单调递增，即最小残差点是一个稳定的点。
InSAR是一种相对测量技术，即任意干涉像对中理论上只能得到相对形变。通常如精河MS6.6地震这样的小形变事件，可以假定远场形变为零，进而得到近场绝对形变量。但由于大气和电离层的影响，得到的形变场中仍然可能含有一定的零漂和线漂成分。考虑一般情况，本文尝试从反演流程中自动估计数据偏移，实现在反演中进行去漂移处理。本文在Zhang等（2012）的InSAR同震滑动反演方法的基础上，采用线性拟合的方式去除数据中的漂移。InSAR观测数据可以拆分为去漂移后的数据和漂移项，即
$$ {{\boldsymbol{Q}}_{{\rm{ob}}}}{\text{＝}}{{\boldsymbol{Q}}_0} {\text{＋}} {{\boldsymbol{Q}}_{{\rm{drift}}}}{\text{，}} $$ (1) 式中，Qob为观测数据，Q0为无漂移数据，Qdrift为漂移大小。Q0可以表示为
$$ {{{\boldsymbol{Q}}_{{0}}}} {\text{＝}}{\boldsymbol{B}} {\boldsymbol{f}} {\text{，}} $$ (2) 式中，B为由断层面上各子断层单位滑动形成的地表形变矩阵，f为包含子断层滑动大小的未知参数。假设漂移是线性的，则Qdrift可以表示为
$$ {{\boldsymbol{Q}}_{{\rm{drift}}}}{\text{＝}}{{\boldsymbol{s}}_0} {\text{＋}} {\boldsymbol{X}}{{\boldsymbol{s}}_x} {\text{＋}} {\boldsymbol{Y}}{{\boldsymbol{s}}_y}{\text{，}} $$ (3) 式中：X和Y分别为平面坐标系中InSAR数据覆盖区域内各点的坐标值；s0，sx和sy分别为拟合参数。则考虑漂移的InSAR同震滑动反演方程即可列为
$$ {{{\boldsymbol{Q}}_{{\rm{ob}}}}}{\text{＝}}\left[\!\!\!\! {\begin{array}{*{20}{c}} {\boldsymbol{B}}&{\boldsymbol{I}}&{\begin{array}{*{20}{c}} \!\!\!\!{\boldsymbol{X}}&{\boldsymbol{Y}} \end{array}} \end{array}} \!\!\!\!\!\!\!\!\right]\left[\!\!\!\!\!\!\!\! {\begin{array}{*{20}{c}} {\boldsymbol{f}}\\ {{{\boldsymbol{s}}_0}}\\ {\begin{array}{*{20}{c}} {{{\boldsymbol{s}}_{{x}}}}\\ {{{\boldsymbol{s}}_{{y}}}} \end{array}} \end{array}} \!\!\!\!\!\!\!\!\!\!\right]{\text{，}} $$ 式中，I为元素全为1的列向量。基于该方程，我们可以反演出InSAR数据的漂移项，最终得到去漂移后的InSAR数据。
在去漂移的InSAR反演中，需要给定断层面参数（initial fault-plane parameters，缩写为IFP）。通过给定不同的IFP并对InSAR数据进行去漂移处理后看出，不同的IFP得到的无漂移InSAR数据Q0存在明显的区别。可能的原因是，精河MS6.6地震破裂产生的地表形变量较小，接近InSAR数据的固有误差水平（约0.1 m），即InSAR数据的相对误差较大，使得IFP较小的变动也会导致去漂移结果的明显变化。为解决此问题，我们尝试了多组IFP，得到多组去漂移的InSAR数据，然后对每一组处理后的InSAR数据进行三维格点搜索，得到多个FP值，记为SFP （searched fault-plane parameters），再根据得到的SFP点的分布确定合适的、适用于最终联合反演的SFP。图2给出了基于InSAR数据滑动分布反演的三维格点搜索流程。
![]()
图 2 基于InSAR滑动分布反演的格点搜索流程示意图Figure 2. Flow chart of the grid search based on the InSAR slip inversions反演结果显示，精河地震在断层面上约21 km×14 km范围内产生滑动，为避免边缘效应，我们将断层面尺度适当拓展为30 km×30 km，并沿走向方向和倾向方向分别划分为15个子断层，即每个子断层大小为2 km×2 km。断层面滑动角在90°±45°之间变化。经过若干次反演尝试，我们将震源固定在沿走向方向的第12个子断层上，并采用中国地震台网中心的定位结果（44.27° N，82.89° E）作为震中位置。将准备的256个IFP作为输入：走向取值分别为75°，80°，85°，90°，95°，100°，105°和110°；倾角取值分别为37°，40°，42°，45°，47°，50°，52°和55°；震源在断层倾向方向上的序数分别为8，9，10，11。采用这256组初始断层面参数对原始InSAR数据分别进行去漂移反演，得到256组去漂移的InSAR数据，随后对256组去漂移的InSAR数据进行基于InSAR数据滑动分布反演的三维格点搜索：各组数据IFP走向搜索范围为±16°，间隔1°；倾角搜索范围为±16°，间隔1°；序数搜索范围为±4，间隔为1。搜索完成后，得到256个SFP点的分布。从所有输出结果为该SFP点的IFP点中寻找残差最小的IFP点，并采用该IFP点去漂移得到的InSAR数据用于联合反演。
图3给出了搜索得到的256个SFP点的个数和密度分布，其中每一点的密度计算方法为该点SFP数量加上边与之相接的四个点SFP数量之和的四分之一，再加上顶点与之相接的四个点SFP数量之和的八分之一。从三维搜索空间的三个侧面均能看到，SFP点存在一定的聚集。通过计算得到的256个SFP点在走向、倾角和序数的平均值为[95，47，10]，并将其作为搜索得到的断层面参数。为考察该点是否稳定，我们比较了该点附近256个SFP点的平均残差。结果显示，从三个侧面来看，[95，47，10]在搜索空间中均为其附近的最小残差点，结果具有稳定性。然后，我们找到输出结果为SFP＝[95，47，10]的最小残差的IFP点，即IFP＝[100，40，8]，并将该IFP点的去漂移InSAR数据用于接下来的联合反演。
![]()
图 3 搜索得到的256个SFP点的分布图（a）—（c）分别为倾角与走向、震源所在子断层沿倾向方向序数（简称序数）与走向、序数与倾角剖面所显示的SPF三维空间分布，圆圈表示SFP点，其尺寸对应残差大小；图（d）—（f）分别为与图（a）—（c）对应的SFP分布密度；图（g）—（i）分别为已选取的SFP＝[95，47，10]为中心，搜索空间中的走向截面、倾角截面和序数截面的平均残差Figure 3. Distribution of the searched 256 SFPsFigs.（a）− （c） show the SFP distribution in the cross sections of strike-dip，ordinal number of subfault along dip the epicenter located in （“ordinal”for short）-strike，and ordinal-dip，respectively，where the circles indicate the SFPs. The size of the circle denotesthe residuals，and circle color represents the number of SFPs. Figs.（d）−（f） show the density distribution of SFP correspon-ding to Figs.（a）− （c），respectively. Figs.（g）−（i） show the average residuals in the search space （selected SPF＝[95，47，10] centered） in the cross section of strike，dip and ordinal，respectively3. 破裂过程联合反演
基于搜索流程给出的结果，我们进一步明确联合反演的断层模型参数：断层面走向为95°，倾角为47°，为一个30 km×30 km的矩形，并将其划分为225个2 km×2 km的子断层。震源位于沿走向方向的第12个子断层也是沿倾向方向的第10个子断层的中心点。反演中滑动角允许在90°±45°之间变化。最大破裂速度限定在2.8 km/s以内，允许每个子断层上的最长破裂持续时间为4 s。
我们采用远震数据、近场强震数据和InSAR数据，基于Zhang等（2012）的方法进行震源破裂过程的联合反演。反演中另行添加了空间光滑约束、时间光滑约束（Horikawa，2001；Yagi et al，2004）和最小标量地震矩约束（Hartzell，Iida，1990），采用非负共轭梯度法（Ward，Barrientos，1986）求解反演方程组。经尝试，给定强震数据3倍于远震数据的权重，并赋予InSAR数据和地震波数据同等权重。同时，为获得更好的波形拟合结果，允许地震波形在时间上进行平移：当首次联合反演完成后，对于远震波形，在时间域上对每个台站的观测波形进行不超过5 s的平移，使得观测波形与合成波形的相关系数达到最大；对于强震波形，将同一台站的三分量波形同时进行不超过10 s的平移，使得三分量的观测波形与合成波形的平均相关系数最大。平移完成后，再次联合反演，确定最终的反演结果。
图4展示了联合反演得到的滑动分布和震源时间函数。反演结果显示，精河MS6.6地震为一次单侧破裂事件，滑动集中分布于震源以西5—15 km，从地表沿断层面向下15—25 km，考虑倾角因素，对应为11—18 km的深度范围。最大滑动量约为0.8 m，所处位置在震源以西8 km，深度与震源相近。滑动整体没有出露地表。断层面上的平均滑动角为106°。震源时间函数显示本次地震只有一次主事件，破裂整体持续约9 s：从破裂开始地震矩快速释放，第2.5 s时达到顶峰，随后平缓下行至第9 s。从破裂过程的快照，能够清楚地看到破裂向西传播的过程。反演得到的标量地震矩为3.6×1018 N·m，对应矩震级为MW6.3，与中国地震台网中心、USGS和德国地球科学研究中心（German Research Center for Geosciences，缩写为GFZ）的结果一致。图5显示了地震波形和InSAR数据的拟合情况。
![]()
图 4 基于远震资料、近场强震资料和InSAR资料的地震破裂过程联合反演结果（a） 滑动分布和震源时间函数；（b） 子断层上的子震源时间函数；（c） 每2 s积累的滑动量分布快照Figure 4. Joint inversion results of rupture process based on teleseismic，near-field strong-motion and InSAR data（a） Slip distribution and source time function；（b） Subfault source time functions；（c） Snapshots ofthe slip accumulated at each 2-second interval![]()
图 5 联合反演的观测与合成资料比较（a） 观测（黑线）与合成（红线）远震地震波；（b） 观测（黑线）与合成（红线）近场强震地震波；（c） 观测与合成InSAR数据Figure 5. Comparison of the observed and synthetic data of the joint inversion（a） Observed （in black） and synthetic （in red） teleseismic data； （b） Observed （in black） and synthetic （in red） strong-motion data；（c） Observed and synthetic InSAR data图6给出了地震破裂过程时空图像。主要考虑矩心的破裂速度，由于反演结果的不确定性，为了去除破裂滑动中的非主要信息，在计算每1 s内的矩心位置时，我们将1 s内滑动量小于整个破裂过程最大累积滑动量15%的子断层剔除。同时，考虑到反演结果的不确定性所导致的矩心瞬时速度的不稳定，采用了平均破裂速度（自破裂发生起矩心的累积位移除以已经发生的破裂时间），并计算以每1 s内地震矩大小为权重的加权平均破裂速度。平均破裂速度计算结果显示，破裂速度从第2 s的1.0 km/s缓慢增加至第6 s的2.1 km/s，后缓慢降至第8 s的1.9 km/s；从加权平均破裂速度来看，破裂速度从第2 s的1.0 km/s增加至第5 s的2.6 km/s，后降至第8 s的2.4 km/s。综合来看，本次地震的破裂速度介于2.1—2.6 km/s。
![]()
图 6 矩心破裂传播速度（a） 矩心的位置、时间和平均移动速度，图中时间为破裂开始后经过的时间 ；（b） 矩心的位置、时间和地震矩加权平均移动速度；（c） 破裂传播速度-时间曲线Figure 6. Rupture velocity of the moment centroid（a） The location，time and average velocity of the centroid，the time represents the duration after the rupture；（b） The location，time and moment-weighted average velocity of the centroid；（c） Velocity-time curves4. 讨论
对于精河MS6.6地震这样的中等强度地震，破裂模型反演结果可能会受到数据的影响。因此，采用大地测量数据确定断层面参数，并综合远震地震数据、近场强震数据和大地测量数据进行联合反演，有助于得到更好的破裂模型。
本文也尝试采用三维格点搜索流程对真实的发震断层进行判别（图7）。在对北倾节面的FP进行搜索后，同样得到了256个SFP点的分布。北倾节面对应的SFP相对残差范围为0.081—0.136，平均值为0.096；而南倾节面SFP的相对残差范围为0.080—0.148，平均值为0.098。在北倾节面对应的SFP中选取[261，38，8]作为联合反演使用的参数，该SFP点对应的最小残差为0.094，同时南倾节面中选取的SFP对应的最小残差为0.090。从残差分布范围来看，北倾节面与南倾节面重合度较高；北倾节面平均残差略低于南倾节面，但二者差别较小。综合来看，我们采用的搜索流程难以对发震断层进行有效地判别，其原因是此次地震破裂位置较深，造成的有效地表形变较小，导致InSAR数据的信噪比较低。基于距离破裂区域较远且具有一定误差水平的数据，在区分断层的倾向上存在困难。
![]()
图 7 南倾节面与北倾节面对应的SFP残差分布Figure 7. Residuals of the SFPs for south-dipping and north-dipping nodal planes需要说明的是，一些因素可能会影响联合反演得到的破裂模型。首先，区域沉积层对近场强震台站记录的地震波形有较大影响，导致强震波形拟合结果相较于远震拟合更差。虽然本文采用的CRUST1.0模型包含区域沉积层速度结构信息，但由于该模型的最小分辨尺度相对本次地震破裂尺度大得多，其包含的沉积层速度信息可能与实际速度结构存在一定差别，因而在计算理论格林函数时，我们对区域沉积层速度结构的考虑有可能不够全面。若能获得震中附近区域更为精确的速度结构信息，则有可能提高强震资料的拟合。其次，强震台站大部分位于震中南侧，方位覆盖较差，这样的台站布局可能对反演结果产生一定影响。但从反演结果来看，虽然方位角分布不均匀，但强震数据反演得到的破裂模型与其它数据的反演结果具有一致性；同时，强震数据有助于更好地约束地震破裂位置和位错大小，在时空上的分辨能力比远震数据更强。因此，即使强震台站分布不均匀，强震数据在破裂过程反演中仍然具有一定的优势。最后，震源定位的结果直接影响破裂模型中断层面上滑动区域的位置和破裂速度的计算结果。因此，若能采用其它方法对震源进行更精确的重定位，则有可能提高反演得到破裂模型的准确程度。
为考察在联合反演中不同类型的数据所起的作用，我们分别采用远震数据、强震数据和InSAR数据进行了单独反演。反演结果如图8所示。其中，强震数据反演、InSAR反演和联合反演显示出较为一致的破裂方向性，即破裂从震源向西传播了10—15 km，而远震反演结果中的滑动集中在震源附近，未显示明显的破裂方向性，这与张勇等（2017）的结果相似。这表明，对于中等强度地震产生的破裂，远震P波（本文中为震中距40°以上）可能难以有效地分辨破裂方向性，而距离较近（本文中为震中距3°以内）的强震台站和InSAR数据能够较好地分辨。由于在反演中，InSAR数据对破裂区域的绝对位置的约束力较好，因此我们倾向于认为本次地震破裂向西传播。虽然在本震例中远震数据对于破裂方向性难以分辨，但由于远震地震波穿过地幔，有助于约束震级，另外远震地震波对于断层的不规则性和弯曲不敏感，能够更好地约束断层宏观破裂行为，因此在联合反演中仍起到一定作用。联合反演中，经过多次尝试，我们将远震数据与强震数据的权重比设置为1 ∶ 3。
![]()
图 8 不同类型数据反演得到的滑动分布和震源时间函数（a） 远震、近场强震和InSAR数据联合反演；（b） 远震数据反演；（c） InSAR数据反演；（d） 近场强震数据反演Figure 8. Slip distributions and source time functions of inversions with different data（a） Joint inversion of teleseismic，strong-motion and InSAR data；（b） Inversion of teleseismic data；（c） Inversion of InSAR data；（d） Inversion of strong-motion data我们将联合反演结果与张勇等（2017）的远震快速反演结果、Gong等（2019）的InSAR反演结果、Zhang等（2020）的联合反演结果进行比较。从反演得到的矩震级看，本文结果、张勇等（2017）、Gong等（2019）、Zhang等（2020）依次为MW6.3，MW6.3，MW6.2，MW6.3，结果基本一致；最大滑动量分别为0.8 m，1.1 m，0.3 m和0.6 m，其中远震反演结果最大，InSAR反演结果最小，联合反演结果居中，与我们进行的不同类型数据反演结果一致；从破裂方向性上看，除张勇等（2017）的远震反演结果外，其余反演结果均显示出较为明显的西向破裂，显示远震数据对本次地震破裂方向性的分辨能力较差。此外，Zhang等（2020）的联合反演结果显示破裂向深部扩展，而本文联合反演和Gong等（2019）显示破裂向浅部扩展，这可能是Zhang等（2020）基于北倾节面反演所致，但各反演结果中破裂未出露地表的主要特征是一致的。
在联合反演中，我们在初次反演后，将观测地震波形在时间上进行平移，使得合成波形与相关波形之间具有最大的相关系数，并在此基础上再次进行反演得到最终结果。这一操作的目的是只拟合波形信息，而不考虑到时误差。这有助于减小理论到时的不准确性对反演结果的干扰，也是震源波形反演中常见的操作（Yagi et al，2004）。对于远震数据，由于地震波传播距离较远，到时通常存在最大1—3 s的不确定性，因此在远震数据反演中允许波形时间平移是较为合理的。对于强震数据，由于精河地震发生区域沉积层较厚，浅层速度结构对强震波形有较大的影响，不但影响到时的准确性，还导致强震数据质量偏低，这也是我们在32个强震台站中最终只挑选出11个台站进行反演的原因之一。本文中，我们允许强震波形的时间平移，即选择放弃使用到时信息，以提高波形的拟合程度。同时注意到，在联合反演中InSAR数据对破裂绝对位置的约束较好，部分补偿了强震数据到时信息缺失导致的对破裂绝对位置分辨能力下降的问题。因此，允许强震波形时间平移，能尽可能地降低到时信息不准确造成的负面影响。
精河MS6.6地震主震破裂分布与余震分布可能存在一定的互补关系。图1显示了主震滑动分布和余震分布在地表上的投影。我们共选取2017年8月14日16:00之前中国地震台网中心记录到的209次余震事件，其位置已由白兰淑等（2017）采用双差定位法（Waldhauser，Ellsworth，2000）进行重定位确定。从地表上看，大部分余震集中在破裂区域的西边，小部分余震集中在破裂区域的西南角。余震分布和主震破裂分布的互补性说明主震对余震具有很强的触发作用。主震破裂后，在破裂区域周围引起较强的应力扰动和库仑应力增加，可能是震后余震密集发生的主要原因（Harris，1998；Kilb et al，2002）。
5. 结论
本文中，我们利用InSAR同震形变数据，设计了一套基于滑动分布反演的格点搜索流程，得到了2017年8月9日新疆精河MS6.6地震断层面的走向、倾角和震源深度；在此基础上，利用远震资料、近场强震资料和InSAR同震形变资料，联合反演得到了此次地震的震源破裂时空过程。根据搜索和反演结果，本次地震断层面走向为95°，倾角为47°，震源深度为14 km。整个地震主要以单侧破裂为主，破裂过程持续约9 s，破裂朝西以2.1—2.6 km/s的速度传播，最终在震源以西5—15 km、沿倾向15—25 km范围形成了一个主要的滑动区域，最大滑动量约0.8 m，平均滑动角为106°。整个破裂过程释放的标量地震矩为3.6×1018 N·m，对应矩震级为MW6.3。由于此次地震的破裂深度主要集中在10 km以下，未来可能需要关注震源区浅部0—10 km处发生潜在地震的可能。
本文研究所用的远震数据下载自美国地震学研究联合会IRIS，强震资料由中国地震局工程力学研究所国家强震动台网中心提供，InSAR资料来自于欧洲空间局哨兵-1雷达数据。作者在此对以上机构一并表示感谢。
-
![]()
图 1 三河—唐山短周期密集地震测线分布图
F1：黄庄—高丽营断裂；F2：顺义—良乡断裂；F3：通县—南苑断裂；F4：夏垫断裂；F5：香河断裂；F6：丰台—野鸡坨断裂；F7：唐山断裂；F8：宁河—昌黎断裂（引自刘保金等，2009，2011），下同
Figure 1. Distribution map of short period dense seismic lines in Sanhe-Tangshan
F1：Huangzhuang-Gaoliying fault；F2：Shunyi-Liangxiang fault；F3：Tongxian-Nanyuan fault；F4：Xiadian fault；F5：Xianghe fault；F6：Fengtai-Yejituo fault；F7：Tangshan fault； F8：Ninghe-Changli fault （after Liu et al，2009，2011），the same below
![]()
![]()
![]()
图 4 三河地震区地壳浅部H/V谱比法二维伪深度剖面图与地震反射结果图对比
剖面位置为图1剖面线的西段，白色虚线为本文获得的三河地震区地壳浅部H/V谱比法二维伪深度剖面图中沉积层厚度走向
Figure 4. Comparison between 2D pseudo-depth profile by H/V spectral ratio method and seismic reflection result map of shallow crust in Sanhe seismic region
（The profile location is the west section of the profile line shown in Fig. 1，the white dashed line shows the thickness trend of sedimentary layer in the shallow crust 2D pseudo-depth profile by H/V spectral ratio method obtained in the Sanhe seismic region）
![]()
图 5 唐山地震区地壳浅部H/V谱比法二维伪深度剖面图（a）与地震反射结果图（b）（刘保金等，2011）对比
剖面位置为图1剖面线的东段，图（a）中紫色虚线为Bao等（2018）的剖面沉积层厚度走向，白色虚线为本文的沉积层厚度走向
Figure 5. Comparison between 2D pseudo depth profile by H/V spectral ratio method （a） and seismic reflection result map （b）（Liu et al，2011） of shallow crust in Tangshan seismic area
The profile location is the east section of the profile line shown in Fig. 1 ，in Fig. （a），the purple dotted line is the trend of sediment thickness from Bao et al （2018），and the white dotted line is the trend of sediment thickness in this paper
![]()
图 6 唐山断裂、夏垫断裂附近两个观测点的H/V谱比曲线
Figure 6. H/V spectral ratio curves of two observation points
![]()
图 7 研究区场地放大系数（a）和场地易破坏程度（b）分布图
红色圆点代表每个台站的数据；蓝色三角形代表10 km内各台站数据的平均值
Figure 7. Distribution map of site magnification factor （a） and site vulnerability （b） in the studied area
The red dots represent the data of each station；the blue triangles represent the average value of the data from stations within 10 km
-
虢顺民,李志义,程绍平,陈献程,陈孝德,杨主恩,李如成. 1977. 唐山地震区域构造背景和发震模式的讨论[J]. 地质科学,12(4):305–321. Guo S M,Li Z Y,Cheng S P,Chen X C,Chen X D,Yang Z E,Li R C. 1977. Discussion: On the regional structural background and the seismogenic model of the Tangshan earthquake[J]. Chinese Journal of Geology,12(4):305–321 (in Chinese).
高文学, 马瑾. 1993. 首都圈地震地质环境与地震灾害（附图册）[M]. 北京: 地震出版社: 8–9. Gao W X, Ma J. 1993. Seismo-Geological Background and Earthquake Hazard in Beijing Area （Atlas）[M]. Beijing: Seismological Press: 8–9 （in Chinese）.
国家地震局《一九七六年唐山地震》 编辑组. 1982. 一九七六年唐山地震[M]. 北京: 地震出版社: 4–23. Editorial Group of “The 1976 Tangshan Earthquake” , State Seismological Bureau. 1982. The 1976 Tangshan Earthquake[M]. Beijing: Seismological Press: 4–23 （in Chinese）.
刘保金,胡平,孟勇奇,酆少英,石金虎,姬计法. 2009. 北京地区地壳精细结构的深地震反射剖面探测研究[J]. 地球物理学报,52(9):2264–2272. doi: 10.3969/j.issn.0001-5733.2009.09.010 Liu B J,Hu P,Meng Y Q,Feng S Y,Shi J H,Ji J F. 2009. Research on fine crustal structure using deep seismic reflection profile in Beijing region[J]. Chinese Journal of Geophysics,52(9):2264–2272 (in Chinese).
刘保金,曲国胜,孙铭心,刘亢,赵成彬,徐锡伟,酆少英,寇昆朋. 2011. 唐山地震区地壳结构和构造:深地震反射剖面结果[J]. 地震地质,33(4):901–912. doi: 10.3969/j.issn.0253-4967.2011.04.014 Liu B J,Qu G S,Sun M X,Liu K,Zhao C B,Xu X W,Feng S Y,Kou K P. 2011. Crustal structures and tectonics of Tangshan earthquake area:Results from deep seismic reflection profiling[J]. Seismology and Geology,33(4):901–912 (in Chinese).
刘渊源,崇加军,倪四道. 2011. 基于井下摆天然地震数据测量首都圈近地表波速结构[J]. 地震学报,33(3):342–350. Liu Y Y,Chong J J,Ni S D. 2011. Near surface wave velocity structure in Chinese capital region based on borehole seismic records[J]. Acta Seismologica Sinica,33(3):342–350 (in Chinese).
彭菲,王伟君,寇华东. 2020. 三河—平谷地区地脉动H/V谱比法探测:场地响应、浅层沉积结构及其反映的断层活动[J]. 地球物理学报,63(10):3775–3790. doi: 10.6038/cjg2020O0025 Peng F,Wang W J,Kou H D. 2020. Microtremer H/V spectral ratio investigation in the Sanhe-Pinggu area:Site responses,shallow sedimentary structure,and fault activity revealed[J]. Chinese Journal of Geophysics,63(10):3775–3790 (in Chinese).
彭一民,李鼎容,谢振钊,王安德,刘清泗. 1981. 北京平原区同生断裂的某些特征及其研究意义[J]. 地震地质,3(2):57–64. Peng Y M,Li D R,Xie Z Z,Wang A D,Liu Q S. 1981. Some features of contemporaneous faults in Beijing plain and their significance[J]. Seismology and Geology,3(2):57–64 (in Chinese).
王伟君,陈棋福,齐诚,谭毅培,张项,周青云. 2011. 利用噪声HVSR方法探测近地表结构的可能性和局限性:以保定地区为例[J]. 地球物理学报,54(7):1783–1797. Wang W J,Chen Q F,Qi C,Tan Y P,Zhang X,Zhou Q Y. 2011. The feasibilities and limitations to explore the near-surface structure with microtremor HVSR method:A case in Baoding area of Hebei Province,China[J]. Chinese Journal of Geophysics,54(7):1783–1797 (in Chinese).
张宝龙,李志伟,包丰,邓阳,游庆瑜,张森琦. 2016. 基于微动方法研究五大连池火山区尾山火山锥浅层剪切波速度结构[J]. 地球物理学报,59(10):3662–3673. Zhang B L,Li Z W,Bao F,Deng Y,You Q Y,Zhang S Q. 2016. Shallow shear-wave velocity structures under the Weishan volcanic cone in Wudalianchi volcano field by microtremor survey[J]. Chinese Journal of Geophysics,59(10):3662–3673 (in Chinese).
Agostini L,Boaga J,Galgaro A,Ninfo A. 2015. HVSR technique in near-surface thermal-basin characterization:The example of the Caldiero district (North-East Italy)[J]. Environ Earth Sci,74(2):1199–1210. doi: 10.1007/s12665-015-4109-0
Bao F,Li Z W,Yuen D A,Zhao J Z,Ren J,Tian B F,Meng Q J. 2018. Shallow structure of the Tangshan fault zone unveiled by dense seismic array and horizontal-to-vertical spectral ratio method[J]. Phys Earth Planet Inter,281:46–54. doi: 10.1016/j.pepi.2018.05.004
Bard P Y,Campillo M,Chávez-Garcia F J,Sánchez-Sesma F. 1988. The Mexico earthquake of September 19,1985:A theoretical investigation of large- and small-scale amplification effects in the Mexico City valley[J]. Earthq Spectra,4(3):609–633. doi: 10.1193/1.1585493
Guo H,Jiang W L,Xie X S. 2011. Late-Quaternary strong earthquakes on the seismogenic fault of the 1976 MS7.8 Tangshan earthquake,Hebei,as revealed by drilling and trenching[J]. Science China Earth Science,54(11):1696–1715. doi: 10.1007/s11430-011-4218-x
Ibs-von Seht M,Wohlenberg J. 1999. Microtremor measurements used to map thickness of soft sediments[J]. Bull Seismol Soc Am,89(1):250–259. doi: 10.1785/BSSA0890010250
Jemberie A L,Langston C A. 2005. Site amplification,scattering,and intrinsic attenuation in the Mississippi embayment from coda waves[J]. Bull Seismol Soc Am,95(5):1716–1730. doi: 10.1785/0120040203
Konno K,Ohmachi T. 1998. Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components of microtremor[J]. Bull Seismol Soc Am,88(1):228–241. doi: 10.1785/BSSA0880010228
Nakamura Y. 1989. A Method for Dynamic Characteristics Estimation of Subsurface Using Microtremor on the Ground Surface[R]. Tokyo: Railway Technical Research Institute: 25–30.
Nakamura Y. 1997. Seismic vulnerability indices for ground and structures using microtremor[C]//World Congress on Railway Research. Florence, Italy: World Congress on Railway Research: 1–7.
Nakamura Y. 2009. Basic structure of QTS （HVSR） and examples of applications[C]//Increasing Seismic Safety by Combining Engineering Technologies and Seismological Data. Dordrecht: Springer: 33–51.
Ni S D,Li Z W,Somerville P. 2014. Estimating subsurface shear velocity with radial to vertical ratio of local P waves[J]. Seismol Res Lett,85(1):82–90. doi: 10.1785/0220130128
Parolai S,Bormann P,Milkereit C. 2002. New relationships between vS,thickness of sediments,and resonance frequency calculated by the H/V ratio of seismic noise for the Cologne area (Germany)[J]. Bull Seismol Soc Am,92(6):2521–2527. doi: 10.1785/0120010248
Seale S H,Archuleta R J. 1989. Site amplification and attenuation of strong ground motion[J]. Bull Seismol Soc Am,79(6):1673–1696.
Withers M,Aster R C,Young C,Beiriger J,Harris M,Moore S,Trujillo J. 1998. A comparison of select trigger algorithms for automated global seismic phase and event detection[J]. Bull Seismol Soc Am,88(1):95–106. doi: 10.1785/BSSA0880010095
-
期刊类型引用(1)
1. 郑兴群，陶正如，白凯. 面向地震动估计需求的区域传播介质参数. 地震地质. 2024(05): 1091-1105 . 
百度学术
其他类型引用(4)
下载:
计量
- 文章访问数: 202
- HTML全文浏览量: 78
- PDF下载量: 85
- 被引次数: 5
