张家口-宣化盆地位于张渤带和山西断陷带交汇部位，盆地内地震活动受张家口断裂和宣化盆地南缘断裂控制，区内曾发生一系列中强以上地震，但目前该盆地及周边区域地震活动特征表现为中等地震活动超常平静。国内学者利用地质学方法、水准和GPS观测数据等对张家口断裂和宣化盆地南缘断裂的滑动速率和活动性开展相关研究^[1-5]，但未给出2条断裂沿断层面的倾滑速率和闭锁深度，无法对张家口-宣化盆地中-长期地震危险性分析提供参考。

InSAR技术作为一种新兴的大地测量手段，具有全天候、全天时、覆盖广、高空间分辨率和高测量精度等优势，可获取连续的空间地表形变场^[6]。本文使用2018~2020年升降轨模式下3条轨道的Sentinel-1数据，利用PS-InSAR技术结合先验条件的最小二乘迭代逼近法，获取张家口-宣化盆地地表三维形变速率场，并以此为基础反演张家口断裂和宣化盆地南缘断裂的倾滑速率和闭锁深度，为区域构造模式研究和中-长期地震危险性评价提供依据。

1 数据与处理方法 1.1 InSAR数据处理

本文选用欧空局3条轨道的Sentinel-1A和Sentinel-1B雷达影像数据，成像波段为C波段，TOP模式为干涉宽幅(IW)模式，时间跨度为2018~2020年，影像参数见表 1。47轨道、120轨道和142轨道分别选取2019-03-10、2018-11-27和2019-06-03的SLC影像作为主影像，空间基线见图 1。采用ISCE软件中topsStack数据处理模块进行影像配准和差分干涉处理，其中影像配准采用综合考虑主辅影像脉冲响应函数和重叠区域偏差的增强谱分集方法(ESD)，距离向和方位向多视比为4 ∶1，利用StaMPS软件实现干涉图的时序分析^[7]。InSAR外界误差源主要包括大气误差、轨道误差和地形误差等，本文采用PS-InSAR方法提取雷达波后向散射较强且时序上相对稳定的地物目标，数据信噪比较高，大气误差较小；C波段雷达卫星发射和接收电磁波主要受对流层影响而产生相位延迟，在PS点相位解缠之前使用GACOS大气延迟相位模型削弱对流层引起的干涉误差；分别使用精密轨道星历数据AUX_POEORB和30 m分辨率的STRM高程模型数据改进轨道精度和模拟地形相位，以提高永久散射体的形变速率精度^[7]。最终得到3条轨道的LOS向最优形变速率。

1.2 三维形变场构建

由于Sentinel-1卫星升降轨使用同一卫星平台，成像时各物理参数几乎未发生变化，因此影像质量一致性较好。升轨和降轨SAR卫星成像的最大差异为雷达电磁波束入射角和卫星飞行方向方位角不同，根据该特征，构建地表形变监测结果的三维分解模型^[8]。InSAR观测结果为地表真实三维形变在LOS向的投影，同时LOS向形变又可被分解到NS向、EW向和UD向，其几何关系式为：

$ \begin{gathered} d_{\mathrm{LOS}}=D_U \cos \theta-D_N \sin \theta \cos (\alpha-3 {\rm{ \mathsf{ π} }} / 2)- \\ D_E \sin \theta \sin (\alpha-3 {\rm{ \mathsf{ π} }} / 2) \end{gathered} $

(1)

式中，d_LOS为LOS向观测值；D_N、D_E和D_U分别为NS向、EW向和UD向形变值(约定U向、N向和E向为正)；θ为雷达侧视角；α为卫星轨道方位角。

理论上，融合3个及更多平台或轨道的InSAR形变观测值就可对地表三维形变场进行构建。本文使用结合先验条件的最小二乘迭代逼近法，该方法利用已知条件或假设可信已知条件作为多余观测，迭代使用最小二乘方法求解矛盾方程，以获取变量近似解^[8]。模型方程及约束条件为：

$ \boldsymbol{F}_{n \times 3} \times \boldsymbol{D}_{3 \times 1}=\boldsymbol{d}_{n \times 1} $

(2)

$ \boldsymbol{B}_{n \times 3} \times \boldsymbol{D}_{3 \times 1}+\boldsymbol{W}_{n \times 1}=\bf{0} $

(3)

式中，d_n×1为n组视线向观测值组成的矩阵；F_n×3为投影系数矩阵，由SAR成像几何条件确定；W_n×1为约束条件；B_n×3为条件系数矩阵。

现有SAR卫星多采用近极地轨道飞行和侧视成像模式，LOS向形变观测值在三维分解过程中对NS向形变信息非常不敏感，对UD向形变敏感度最高。可将NS向位移值向极小逼近作为初始约束条件，结合式(2)进行最小二乘解算，根据LOS向形变观测值中三维分量贡献程度和解算精度的高低，依次选择可信度最高的分量作为约束条件，使用最小二乘解算得到研究区真实可靠的三维形变场。

2 InSAR三维形变场分析

本文利用PS-InSAR技术和结合先验条件的最小二乘迭代逼近法获取张家口-宣化盆地及周边区域三维形变速率场。

2.1 三维形变结果与检验

使用文献[9-10]中2009~2015年水平向和垂直向GPS速率与基于InSAR观测数据并结合先验条件的最小二乘迭代逼近法获取的三维形变速率进行对比，结果如表 2所示。由表可知，二者差值大多在2 mm/a以内，这可能是由局部地表变形在不同监测时段的运动速率差所致。统计结果表明，NS向、EW向和UD向的均方根误差分别为0.76 mm/a、1.64 mm/a和1.51 mm/a，说明本文三维解算结果总体上具有可靠性。

2.2 三维形变场空间变化特征

图 2(a)为张家口-宣化盆地NS向形变速率场，由图可知，形变速率整体为负，但绝对值较小，盆地所在区域表现为整体向南的微弱运动趋势，且盆地内部由西向东N向运动速率逐渐变大。图 2(b)为张家口-宣化盆地EW向形变速率场，其中张家口断裂北侧和宣化盆地南缘断裂南侧区域整体为正，仅在张家口市区北部、张家口断裂与宣化盆地南缘断裂之间区域为负，葛峪堡村、南滩村附近区域为正且绝对值较大，经野外实地调查发现均为人为干扰。图 2(c)为张家口-宣化盆地UD向形变速率场，其中张家口断裂中段北侧区域整体呈微弱抬升趋势，断裂东段北侧区域整体为沉降运动，沉降量最大区域为崇礼区东坪黄金矿区附近；张家口断裂和宣化盆地南缘断裂之间区域整体为正，且盆地东部形变速率值大于西部，但在陈家堡村和葛峪堡村分别受地下水开采和土方开采等人为干扰导致局部发生显著沉降，南滩村附近山体上大面积安装太阳能板使地表形变表现为显著抬升；宣化盆地南缘断裂南部区域整体表现为沉降趋势。

张家口断裂中段为该断裂主体部分且构造行迹清晰。为进一步分析张家口-宣化盆地、张家口断裂中段及宣化盆地南缘断裂的形变速率变化特征，在跨断裂两侧12 km近场区域选取3条垂直于断裂走向的速率剖面，剖面宽1 km，均跨越活动断裂且尽可能远离人为干扰较大区域。如图 2所示，AA1剖面位于张家口断裂中段(西太平山)，长21 km；BB1剖面线位于张家口断裂中段(人头山村)，长24 km；CC1剖线位于宣化盆地南缘断裂(龙门坡村东)，长24 km。

根据InSAR形变速率剖面，提取跨断层两侧的形变速率均值，估算得到张家口断裂中段和宣化盆地南缘断裂LOS向形变速率(图 3~5)。由图可知，张家口断裂中段两侧LOS向跨断裂形变速率差异为0.19~0.66 mm/a，宣化盆地南缘断裂两侧速率差异为0.01~0.66 mm/a，在断层近场区域均表现出较小的运动差异。根据反演得到的张家口-宣化盆地三维形变速率剖面，分别提取跨断层两侧NS向、EW向和UD向平均形变速率，由于断层附近形变点位数量较大，部分点位因误差过大并不能真实反映断层实际运动状态，因此仅保留3倍标准差以内的点位。可以看出，3条剖面在断层近场处均有微小差异，其中张家口断裂中段和宣化盆地南缘断裂跨断层水平向速率剖面显示，断层两侧地表存在微小形变差异。由NS向跨断层形变速率最佳拟合曲线(图 6(a)~8(a))可知，盆地南部区域和张家口断裂中段北部区域表现为微弱拉张运动，这与张家口-宣化盆地NNW-SSE向拉张性质水平构造应力场特征基本一致^[11]；由EW向跨断层形变速率最佳拟合曲线(图 6(b)~8(b))可知，张家口断裂中段和宣化盆地南缘断裂两侧地表相对运动存在微小差异；将断层两侧NS向和EW向形变速率投影到平行断层走向方向，最佳拟合曲线(图 6(c)~8(c))结果显示，张家口断裂中段呈右旋走滑运动特征，速率为0.29~0.57 mm/a，宣化盆地南缘断裂呈左旋走滑运动特征，速率为1.01 mm/a。图 9为3条剖面跨断层UD向形变速率最佳拟合曲线，可以看出，张家口-宣化盆地内部和张家口断裂中段北侧区域构造活动整体表现为抬升，宣化盆地南缘断裂南侧区域则表现为沉降；图 9(a)显示张家口断裂中段在该处两盘均呈上升趋势，且南盘相对北盘下降，表现为正断特征；图 9(b)显示张家口断裂中段在该处南盘相对北盘呈微弱上升趋势，表现为逆断特征；图 9(c)显示宣化盆地南缘断裂北盘上升，南盘下降，表现为正断特征。

3 断裂滑动速率和闭锁深度反演

张家口断裂和宣化盆地南缘断裂构造运动以垂向为主^[11]，本文基于跨断层剖面提取的跨断层两侧垂向平均形变速率，采用由刃型位错提出的倾滑断层震间形变数学表达式进行剖面反演^[12]。反演模型数学表达式为：

$ u_x=\frac{s}{{\rm{ \mathsf{ π} }}}\left[\cos \delta \tan ^{-1}\left(\frac{x}{d}\right)+\frac{\sin \delta-\frac{x}{d} \cos \delta}{1+\left(\frac{x}{d}\right)^2}\right] $

(4)

$ u_z=-\frac{s}{{\rm{ \mathsf{ π} }}}\left[\sin \delta \tan ^{-1}\left(\frac{x}{d}\right)+\frac{\cos \delta+\frac{x}{d} \sin \delta}{1+\left(\frac{x}{d}\right)^2}\right] $

(5)

式中，u_x为水平面垂直断层方向位移，u_z为垂直地面方向位移，s为沿断层面倾滑位移，d为断层闭锁深度，δ为断层倾角，x为观测点到断层的垂直距离。

由于式(5)高度非线性，无法确定参数的近似值，直接使用最小二乘法进行反演容易使方程产生奇异并无解。为解决该问题，本文使用拟牛顿法和通用全局优化法来确定倾滑速率、闭锁深度和断层倾角3个参数的最优估计结果。在反演过程中，根据区域历史地震的震源深度，将断层闭锁深度范围设置为1~30 km；根据已有的断层结构研究结果^[1]，将断层倾角范围设置为55°~65°；研究区内断裂平均垂直滑动速率小于1 mm/a，因此将断层倾滑速率搜索范围设置为0.1~1 mm/a^[2]。利用上述方法反演得到的张家口断裂中段和宣化盆地南缘断裂的断层参数如表 3所示，断层闭锁深度分布如图 10所示。

4 讨论与分析 4.1 与地质、水准和GPS结果进行对比

基于由刃型位错提出的倾滑断层震间形变数学表达式，利用2018~2020年获取的InSAR垂直形变速率场数据进行拟合，反演得到张家口断裂中段和宣化盆地南缘断裂的倾滑速率分别为0.42~0.79 mm/a和0.51 mm/a，与前人采用地质学方法和水准观测得到的断层垂直滑动速率结果^[1-3]基本一致。

使用剖面投影法分别获得断层两侧水平向运动差异特征，其中张家口断裂中段和宣化盆地南缘断裂的水平向滑动速率分别为0.29~0.57 mm/a和1.01 mm/a，与陈长云等^[4]利用1999~2007年GPS观测数据计算得到的张家口断裂左旋走滑速率结果基本一致，而本文计算的张家口断裂中段水平向走滑运动特征与已有研究结果相反^{[2, 4]}。张家口-宣化盆地所处的华北地区在印度洋板块和太平洋板块的共同作用下，现代构造应力场表现为以NEE-EW向挤压为主，兼NNW向拉张，使得NW向和近EW向断裂表现为左旋错动，但区域构造应力场在演化过程中会出现应力调整作用^[13]。

洞体应变观测是监测地壳运动和变形的重要手段之一，主要用来测量地壳表面两点之间的应变量(观测曲线上升为拉张，下降为挤压)，进而揭示弹性应变能在局部区域逐步积累和应力、应变逐步演变特征。本文利用位于张家口断裂中段北侧附近的张家口地震台洞体应变仪2017~2020年观测数据，该仪器NS分量观测数据于2019-07-13~2020-04-01出现与往年同期不同的区域应力异常现象，且持续时间较长，局部区域构造应力场在NS向呈挤压变化趋势，并持续至2020-05-30(图 11)。张家口台洞体应变仪NS分量观测数据反映出张家口-宣化盆地北部局部区域构造应力场在2019~2020年发生变化，导致NWW向张家口断裂中段在2018~2020年沿断层走向方向表现出微弱的右旋走滑特征。

4.2 断裂滑动速率、闭锁深度与地震危险性

张家口断裂中段和宣化盆地南缘断裂的闭锁深度分别为9.05~12.98 km和8.38 km，与张希等^[14]使用负位错模型反演得到的鄂尔多斯块体及其周缘断层的闭锁深度结果基本一致，说明2条断裂均具有发生中强以上地震的应变能积累背景条件。

5 结语

本文利用PS-InSAR技术分别选用3条轨道(升轨142、降轨47和降轨120)的Sentinel-1数据获取张家口-宣化盆地及周边区域地表形变速率场，并利用结合先验条件的最小二乘迭代逼近法计算盆地及周边区域三维地表形变场。分别使用剖面投影法和由刃型位错提出的倾滑断层震间形变数学表达式计算和分析张家口断裂中段及宣化盆地南缘断裂现今水平向滑动速率、倾滑速率及断层闭锁深度，得出以下结论：

1) 张家口断裂中段和宣化盆地南缘断裂的水平向滑动速率分别为0.29~0.57 mm/a和1.01 mm/a，倾滑速率分别为0.42~0.79 mm/a和0.51 mm/a；张家口断裂中段表现为右旋走滑兼正(逆)断倾滑活动特征，宣化盆地南缘断裂表现为正断倾滑兼左旋走滑活动特征。

2) 张家口断裂中段和宣化盆地南缘断裂的断层闭锁深度分别为9.05~12.98 km和8.38 km，二者相差较小，均具有发生中强以上地震的应变能积累背景条件。