地面沉降是北京平原区主要的地质灾害,其沉降范围和幅度呈现逐年扩大的趋势^[1]。研究表明,不合理开采地下水是诱发北京平原区发生沉降的主要因素^[2,3,4]。除此之外,密集的高大建筑群、高速立体交通网所带来的地表动静载荷的增加,加速了平原区不均匀地面沉降的发生^[5]。已有研究表明,自北京平原出现地面沉降灾害以来就呈现出不均匀的时空特征^[6,7,8,9,10]。宫辉力等^[1]结合北京平原地质资料,初步分析了黄庄高丽营、南口孙河断层两侧不均匀沉降与地层、构造的关系; 杜钊锋等^[6]采用坡度分析方法获取了北京北部平原区不均匀沉降的分布特征,发现地下水位变化导致土壤受压不匀形成不均匀沉降; 陈蓓蓓等^[7]发现不同的浅地表空间利用对不均匀沉降的影响存在差异; 周超凡等^[11]采用莫兰指数方法对北京地面沉降的空间分布格局特征进行了分析,发现北京平原区地面沉降存在过渡带,是不均匀沉降的主要发生地。已有研究主要集中于不均匀沉降的演化特征和形成原因,但对于沉降不均匀程度的差异研究较少。

为了研究北京平原区不均匀沉降态势,采用雷达干涉测量(synthetic aperture Radar interferometry,InSAR)技术获取北京平原区的地面沉降信息。InSAR可以有效地获得地面沉降的时空信息,并且已在全球许多地面沉降灾害城市成功应用^[11]。Ferrett等^[13,14,15]提出的永久散射体干涉测量(permanent scatterer InSAR,PSInSAR)方法,弥补了传统InSAR技术的时空失相关、大气延迟等缺陷,能得到更加精确的沉降信息,且在长时间缓慢累积地表形变监测中有独特的优势。本文采用48景Radarsat-2数据,基于PSInSAR技术获取了北京平原区2011—2015年间的地表形变信息。采用莫兰指数定量表征北京平原典型沉降区的不均匀沉降程度,结合交叉小波探究地下水及浅地表空间利用差异对沉降不均匀程度的影响。

北京市位于华北平原的西北部边缘,整个地势西北高、东南低。北京平原为典型的冲洪积缓倾平原,由永定河、潮白河、温榆河和泃河等河流携带的松散物质在山前基底上沉积而成。山前倾斜地带多为黏砂碎石、砂卵石组,自山前向平原沉积物厚度由薄变厚,颗粒由粗变细,由单一砂砾石层逐渐变为多层砂与黏土互层。本文以北京平原区为研究区,其数字高程模型(digital elevation model,DEM)如图1所示。

为了完整地获取2011—2015年间北京平原的地面沉降信息,选用在同等空间分辨率模式下,影像覆盖面积较大的48景Radarsat-2降轨数据,见图1。Radarsat-2是30 m空间分辨率的宽幅数据,采用降轨右视成像,SS&Dual极化方式,时间跨度从2010年11月22日—2015年11月20日期间。

PSInSAR方法通过识别区域内的永久散射体(permanent scatterers,PS),对这些PS点进行时间序列分析,消除大气的影响,准确获取区域地表形变信息。本次研究中,为进一步减小大气延迟相位所引起的误差,将气象数据添加到大气校正分析中,可更精确地分离出大气延迟相位,获得高时空分辨率的地面沉降结果,提高形变测量的精度与效率。PSInSAR方法原理公式为

式中: φ_flat为平地相位; φ_topo为DEM误差带来的残余地形相位; φ_atm为大气延迟相位; φ_noise为热噪声、配准误差等带来的噪声; φ_def为地面沉降相位。去除其他4种相位即可获得地面沉降信息。

北京平原区地面沉降时空信息获取的具体步骤如图2所示。

地理学第一定律指出,任何事物都是与其他事物相关的,相近的事物关联更紧密^[16]。空间自相关分析即要对地理样本之间的相关程度进行表征。若所测值随样本间距离减小相似性越来越大,则样本间关系为空间正相关; 若所测值随样本间距离减小差异越来越大,则称为空间负相关; 若所测值随样本间距离的改变无变化,则为空间随机分布。空间自相关分析通过空间自相关指数的计算,使得对地理数据的空间自相关分析由定性变为定量描述。主要的空间自相关指数之一即为全局莫兰指数(Moran’s I)。其表达式为

式中: x_i为变量{x_i}在位置或区域i的值; x̅为变量的均值; S²为变量的方差; S₀为所有变量空间权重之和; n为变量观测值的总个数,也是观测值所对应区域或位置的总个数; W_ij为空间权重矩阵W中的元素,指区域或位置i与j之间的空间权重,设定若空间单元i与j相邻时,W_ij=1; 否则W_ij=0。

Moran’s I一般在[-1,1]之间,[-1,0)代表负相关的程度,(0,1]代表正相关的程度,值越大空间自相关程度越大,若Moran’s I指数为0,则地理变量的空间分布为随机状态。在验证统计学分析结果是否有意义时,前提是提出一个0假设。0假设代表统计分析结果为随机分布,莫兰指数通过P值和Z得分决定是否拒绝0假设,P值代表数据集的可靠性,Z得分和莫兰指数的I值表明数据集有明显的聚集或离散规律。表1为未经验证的P值、Z得分和置信度之间的关系。

本文通过空间自相关分析(Global Moran’s I)实现莫兰指数的计算,采用距离关系中的设定距离生成空间权重矩阵,地理变量为2011—2015年间每个典型区域各年累计沉降量和时序沉降量。

对覆盖北京平原区的48景Radarsat-2影像(2011—2015年)进行差分干涉处理和PS点提取,最终获得高空间密度的PS点共210 107个,密度为33个/km²,平原区年均沉降速率范围为-129.97~12.73 mm/a。

2011—2015年间北京平原区累计沉降量的变化如图4所示,其中最大累计沉降量达到667.85 mm,最小累计沉降量为70.88 mm。

选取2011—2013年间30个水准点的观测结果对同期InSAR监测结果进行检验得到,均方根误差最大为22.74 mm/a,最小为2.13 mm/a,整体均方根误差在10.48 mm/a左右。对InSAR结果与水准测量结果做回归分析,R²达到0.986 1,表明InSAR反演结果与实测地面水准点测量结果相关性很强,如图5所示。

将2011—2015年间累计沉降量大于300 mm的区域作为沉降漏斗区,以沉降漏斗中心作为圆心,半径为6 km,划定5个面积大小相同区域(118 km²)为典型区,分别位于北京平原西北部、中部和东南部,具体分布位置如图6所示。

为探讨沉降的不均匀特征,分别计算了5个典型区2011—2015年间累计沉降量和时序沉降量的莫兰指数发现: 5个典型区的Z得分均大于临界值+1.96,P值均小于显著水平0.05,表明莫兰指数结果置信度大于95%,拒绝0假设。Moran’s I结果均为正值,且达到了0.94以上,表明5个典型区2011—2015年间每一年的累计沉降量和时序沉降量存在正的空间自相关,且各典型区累计沉降量和时序沉降量在空间分布上存在显著集聚的现象。由图7(a)可知,2011—2015年间5个典型区的累计沉降量的莫兰指数分别超过0.972,0.969,0.981,0.948和0.985。5 a间各区域累计沉降量的莫兰指数大小呈现相同规律,由大到小均为: I₅>I₃>I₁>I₂>I₄。累计沉降量的莫兰指数越大,表明沉降的空间相似程度越大,地面沉降的不均匀程度越小。

即区域4累计沉降不均匀程度最大,区域5累计沉降不均匀程度最小。由图7(b)可知,2011—2015年间,5个典型区时序沉降量的莫兰指数分别超过0.972,0.969,0.981,0.948和0.985。5 a间各区域时序沉降量的莫兰指数大小亦呈现相同规律,由大到小均为: I₅>I₃>I₁>I₂>I₄。即区域4各年的时序沉降量的不均匀程度均最大,区域5各年时序沉降量的不均匀程度均最小。此规律与累计沉降量的不均匀程度规律相同。累计沉降量与时序沉降量的不均匀程度差异规律相互验证,表明5个区域沉降的总体不均匀程度差异规律是可靠的,即: I₅>I₃>I₁>I₂>I₄。

由图8可知,典型区1在2015年,时序和累计沉降量的莫兰指数均达到最大,分别为0.979和0.981,表明地面沉降在该年呈现的空间关系最强,即地面沉降空间相似程度最大,表明地面沉降不均匀程度最小; 时序和累计沉降量的莫兰指数最小的年份均为2011年,均为0.972,表明典型区1在该年地面沉降空间相似程度最小,即地面沉降不均匀程度最大。典型区2时序和累计沉降量的莫兰指数分别在2014年和2015年达到最大,分别达到0.979和0.978,地面沉降不均匀程度最小; 时序和累计沉降量的莫兰指数最小的年份均为2011年,均为0.969,表明典型区2在该年地面沉降不均匀程度最大。典型区3时序和累计沉降量的莫兰指数均在2014年达到最大,分别为0.985和0.986,地面沉降不均匀程度最小; 时序和累计沉降量的莫兰指数最小的年份均为2011年,均为0.981,表明典型区3在该年地面沉降不均匀程度最大。典型区4时序和累计沉降量的莫兰指数分别在2014年和2015年达到最大,分别为0.956和0.986,地面沉降不均匀程度最小; 时序和累计沉降量的莫兰指数最小的年份均为2011年,均为0.948,表明典型区4在该年地面沉降不均匀程度最大。典型区5时序和累计沉降量的莫兰指数均在2013年达到最大,分别为0.989和0.988,地面沉降不均匀程度最小; 时序和累计沉降量的莫兰指数最小的年份均为2011年,均为0.985,表明典型区5在该年地面沉降不均匀程度最大。

进一步分析不均匀沉降特征变化趋势,2011—2015年间,5个典型区时序沉降和累计沉降的莫兰指数变化虽有波动,但2015年相较于2011年,全局莫兰指数整体均呈现增大趋势,表明5个典型区地面沉降不均匀程度均呈现减弱趋势。

2011—2015年10月间各典型区每月的时序沉降量变化如图9所示。认为曲线波动剧烈的时间段内,时序沉降量的变化剧烈,相应时间段内沉降的不均匀性明显。时间越短,波动越大,不均匀程度越大; 时间越长,波动越小,不均匀程度越小。区域1内,2011年8月—2012年2月、2013年12月—2014年8月、2015年1月—7月3个时间段内沉降变化量波动较大,说明区域1内这3个时间段地面沉降的不均匀性相对较大。同理,区域2在2011年7月—2012年1月、2013年12月—2014年7月、2014年12月—2015年6月3个时间段地面沉降的不均匀性相对较大。区域3在5年内沉降变化相对比较稳定,只在2015年1月—9月沉降的不均匀性相对明显。区域4在2014年1月—7月、2015年1月—10月2个时间段内沉降的不均匀性相对明显。区域5在2014年1月—7月、2015年1月—10月2个时间段内沉降的不均匀性相对较大。

将典型区与浅地表空间利用情况进行叠加分析可得(图10和表2),5个区域浅地表空间利用相对复杂程度由复杂到简单依次是: 区域3、区域2、区域1、区域4和区域5。

5个典型区地质背景情况如表3所示,由表3可知,区域5以单层结构黏性土为主,其他4个区域均为或绝大部分为多层结构土体。区域1可压缩层厚度小于50 m外,其他区域可压缩层厚度均大于50 m,典型区地质条件背景较为一致。区域1和区域4内分别有2条和3条地裂缝,其他3个区域范围内无地裂缝。

按照不均匀程度由大到小的顺序依次分析每个典型区域的不均匀沉降特征。区域4的浅地表空间利用情况较简单,但不均匀性是5个区域中最明显的,说明区域4沉降空间分布的不均匀程度除受空间利用复杂程度影响外,可能受地裂缝影响较大。区域2空间利用复杂程度仅次于区域3,其地质条件与区域3和5相同,但其沉降的不均匀程度大于区域3和5,仅次于区域4,说明复杂的浅地表空间利用加剧了区域2沉降的不均匀程度。区域1空间利用复杂程度为中等,其不均匀程度也为5个区域中的中等,与浅地表空间利用复杂程度呈正相关。区域3空间利用情况最复杂,导致地面沉降整体均呈现严重状态,虽然地面沉降不均匀程度较小,其复杂的空间利用仍然可能加剧地面沉降的发展,仍需重点关注。区域5浅地表空间利用情况为5个区域中最简单的。土层结构虽然与其他4个区域不同,但土层结构简单,地质条件单一,累计沉降量和年均沉降速率均为5个区域中最小,沉降发育相对缓慢且沉降程度相对轻微,相应地其沉降的不均匀程度不明显,不均匀程度为5个区域中最低。

综上,区域1,2和5沉降不均匀程度大小与浅地表空间利用复杂程度呈正相关。其中特别是浅地表空间利用最复杂的区域3不均匀性较低,但整体地面沉降呈现严重态势,需重点关注。区域4受地质条件影响较大。此外,空间相对复杂区域与相对简单区域相比,相对复杂区域的不均匀态势更显著。

在北京平原地区,造成地面沉降的最主要原因是不合理开采地下水。本文通过比较时序沉降变化量与地下水位变化量的波动趋势研究地面沉降时序不均匀性与地下水流场变化的响应关系,通过交叉小波方法验证结果,同时量化2个时间序列间的时频域关系。

区域1内,2011年5月—2012年1月间地下水位变化量的波动趋势与2011年8月—2012年2月间的沉降量变化趋势相同(图11),表明地下水位变化影响了时序沉降量的变化,2个不均匀沉降较明显的时段均在这个时间区间内。2015年1月—12月期间,地下水位变化滞后于沉降量变化1个月左右。第三个不均匀沉降较明显的时段在此区间内,表明此时间段内影响不均匀沉降程度的主要因素不是地下水,存在其他因素加速了沉降的变化。通过交叉小波进行验证,地面沉降变化与地下水位变化2个时间序列波动较大的共同时域范围为2012年1月—2013年8月,在此期间地面沉降滞后地下水变化1.1~2.8个月。

区域2内,地下水位变化在2011年3月—2015年3月间的波动趋势与沉降变化量在2011年7月—2015年6月间的波动趋势大体相同(图12),3个时序不均匀程度较大的时间段均位于此区间内。因此,2011—2015年间区域2时序沉降的不均匀性主要受地下水位变化影响。交叉小波分析结果表明,地面沉降变化与地下水位变化2个时间序列波动较大的共同时域范围为2012年3月—2013年11月,在此期间地面沉降滞后地下水变化3~4.7个月。

区域3内,地下水位变化量在2014年5月—10月间的波动趋势与沉降变化量在2015年1月—9月的波动趋势大体相同(图13),区域3不均匀沉降较明显的时段位于此区间内。交叉小波分析结果表明,地面沉降变化与地下水位变化2个时间序列波动较大的共同时域范围为2012年2月—2015年2月,在此期间地面沉降滞后地下水变化3~4.9个月。

区域4内,地下水位变化在2011年2月—2013年4月和2014年7月—2015年3月的波动趋势分别与沉降变化量在2011年9月—2013年4月和2014年7月—2015年10月的波动趋势大体相同(图14)。其中一个不均匀沉降较明显的时间段位于此区间内。2013年4月—9月和2014年1月—7月间,地下水位与沉降量变化一致,存在共振现象。其中一个不均匀沉降较明显的时间段位于此区间内,地下水位变化直接影响了沉降量的变化。交叉小波分析结果表明,2个时间序列波动较大的共同时域范围为2012年3月—2014年12月,在此期间地面沉降滞后地下水变化2.3~3.5个月。

区域5内,地下水位变化量在2011年4月—12月间和2013年6月—2015年4月的波动趋势分别与沉降变化量在2011年7月—12月和2013年6月—2015年间的波动趋势大体相同(图15)。2个不均匀性较明显的时间段均位于此区间内。交叉小波分析结果表明,2个时间序列波动较大的共同时域范围为2012年3月—2015年2月,在此期间地面沉降滞后地下水变化2.8~4.6个月。

本文在总结国内外学者采用地面沉降监测技术,进行不均匀沉降特征及影响因素研究成果基础上,利用PSInSAR方法获取了2011—2015年间北京平原区地面沉降信息。根据北京平原地面沉降发育情况,选取5个典型地面沉降中心区,结合空间自相关分析方法,选取全局莫兰指数,定量揭示了各典型区差异沉降演化特征; 结合地下水监测数据及浅地表空间利用情况,分析影响因素与沉降不均匀程度的关系。

1)通过空间自相关分析,获得了5个区域2011—2015年间累计沉降量和时序沉降量的莫兰指数,结果表明5个典型区地面沉降的不均匀态势均呈现减弱趋势,其中区域4的不均匀性最大,区域5的不均匀性最小。结合5个区域浅地表空间利用情况进行分析,发现区域4浅地表空间利用虽较简单,但区域内的多条地裂缝加剧了其不均匀态势的发展; 区域3浅地表空间利用虽最复杂,但沉降的不均匀程度较小,整体地面沉降严重,从长时间序列来看复杂的空间利用仍然会加剧该区域地面沉降发展态势; 其他区域不均匀程度与浅地表空间利用复杂程度呈正相关。

2)通过分析5个典型区域地面沉降时序变化量波动趋势,得到了每个区域时序不均匀沉降特征,5个典型区域分别有不均匀态势比较明显的时段; 结合地下水位时序变化,研究了地下水流场与地表形变之间的响应关系; 通过小波分析方法,将2个时间序列共同波动较大的时域范围具体量化,发现地下水流场动态变化是每个区域地面沉降时序不均匀态势的主要影响因素,地下水水位变化与地面沉降变化2条时间序列间呈现共振变化关系或滞后变化关系,即地下水位波动越大,持续时间越短,沉降的时序不均匀程度越大,地下水位波动越小,持续时间越长,沉降的时序不均匀程度越小。

本研究中存在的不足之处主要表现在不均匀沉降影响因素分析方面。由于北京平原区地质背景的复杂性,对地层岩性、厚度及结构特征精细分析不够。后续研究应搜集平原区更多物探、钻孔数据,细化水文地质对沉降不均匀性的影响分析。