-
北京地区最早发生的地面沉降出现在20世纪30年代,位于西单-东单区域。近年来,北京的快速发展需求使得地下水长期超采,导致地面沉降的幅度和范围逐年扩大。2003—2010年的最大年沉降速率达到110 mm/a,最大累计沉降量达到了723 mm,年均沉降速率达到30 mm/a的覆盖区域面积为480 000 000 m2。北京平原区形成海淀苏家坨、昌平沙河-八仙庄、顺义、朝阳来广营、东郊八里庄、大兴榆垡6大沉降区[1]。由于地面沉降严重威胁城市安全,因此需要分析该地区的地面沉降时空演化特征,预测演化趋势。
有关地面沉降时空演化特征方面的研究方法分为时序分析与空间分析两类。时序演化特征分析通常采用典型点时序图法,从原始数据中观察地面沉降在时间序列上的变化特征或者通过统计某区域的年度沉降量来分析时序特征。刘凯斯[2]采用时序排列熵法分析北京地铁1号和6号线的地面沉降时间序列演化特征。段光耀等[3]用Mann-Kendall检验对北京平原区时空变化进行分析,研究历年发生突变现象的机理。空间演化特征分析一般采用剖面分析、梯度分析来分析空间上的不均匀地面沉降特征[4-6]。Zhou等[7]用空间分析方法等扇分析探究2012—2018年北京平原区地面沉降的扩张趋势,发现地面沉降由向东扩张变为向东和向北扩张。Zuo等[8]运用标准差椭圆方法发现北京市地面沉降漏斗的移动,揭示不均匀地面沉降。上述方法分别从时间或空间的角度研究沉降演化特征,在时间和空间上相分离,不能从时空角度发现数据中隐藏的信息和可能存在的规律。本文采用高维数据分析中的主成分分析(principal component analysis,PCA)方法研究地面沉降的时空演化特征,充分利用合成孔径雷达干涉测量(interferometric synthetic aperture radar,InSAR)所得的沉降信息具有长时序、覆盖范围广的优势。
PCA常用于数据降维,本文的应用目的是通过降维挖掘地面沉降的主要时空特征。PCA应用于地学领域中,能够有效地从时空数据中提取出某信号的时间序列与空间分布,主成分分析模式分为6种,其中T模式(temporal mode)[9]已应用在GPS站点数据、电磁测距和潮汐计数据上,用来分离出瞬态形变事件[10-11]。气象领域中,可对得出的多种环流模式进行解释[12],分析气象雷暴日的规律[13]。
T模式时间序列分析能识别多个时间序列之间的相似空间模式[14]。Rudolph等[15]通过时间主成分分析(temporal principal component analysis,TPCA)从InSAR时序数据中提取主要的时间行为模式;Lipovsky[16]运用TPCA提取长时序形变的季节信号;Chaussard等[17]对小范围且量级小的InSAR监测结果进行TPCA分析,得到第一成分为沉降主趋势,第二主成分表征为季节性形变,与承压水空间覆盖范围较一致,第三主成分表现的空间特征与断裂带位置相关;吴玉苗[18]用类似于PCA的经验正交函数得到隧道两个方向的变形时空特征;邹正波等[19]基于重力场数据识别日本地震,并研究2002—2015年的重力场时空变化特征;Jiang等[20]在对沧州中部承压含水层系统的含水层参数和地下水储量变化定量研究过程中,通过多通道奇异谱分析对地表形变和地下水数据中的季节信号进行分离,推算出弹性骨架存储率,此外,分别重构出总地下水储量、可恢复地下水储量和不可逆地下水储量。
综上所述,TPCA可应用于地学领域中,在未知先验知识条件下,提取时空数据中的时间序列和空间分布特征。本文使用TPCA方法来分析2003—2010年的北京平原区地面沉降,定量提取时空特征并进行合理解释。
HTML
-
TPCA是一种处理多因素复杂问题的多元统计方法,其主要用途就是进行数据降维和数据解释。作为统计分析“由表及里”的数学手段,常用于指标评价研究中,近些年也应用于地理数据的模式识别与数据挖掘中。
利用合成孔径雷达(synthetic aperture radar,SAR)数据获取的地面沉降数据具有大量的空间点位、长时间的监测序列等特点。本文将TPCA方法应用到高维度的地面沉降数据中,可理解为对众多观测点的多维时间数据集以时间为变量进行PCA分析,得到数据本身所蕴含的时空特征和规律。
-
通过TPCA构建M×N的矩阵X,如图 1(a)所示,M指观测点个数,代表永久散射体干涉测量技术(permanent scatterer interferometric synthetic aperture radar,PS-InSAR)识别的永久散射体(permanent scatterer,PS)点;N指日期,代表SAR影像的获取时间。本文的数据即图 1(b)中每一个地理位置的PS点,共近10万个点;每个PS点有51个监测数据。
X是M个样本、N个变量构成的原始矩阵,其表达式为:
第i个变量可表示为Xi=(x1i,x2i…xMi)T,通过TPCA对X线性变换后,得到新变量矩阵Z:
其中,
为了用Z1代替原来的N个变量X1,X2…XN,要求Z1尽可能多地反映原来N个变量的信息,信息量大小用方差来描述。方差Var(Z1)越大,表示Z1包含的信息越多。
当变换矩阵满足UU
=I,存在满足条件的U1使得Var(Z1)达到最大,Z1就称为第一主成分(或者主分量);使得方差Var(Z2)达到第二最大,则Z2成为第二主成分,以此类推求得剩余成分。具体过程如下:(1)对矩阵X均值化:求矩阵X的均值E(X)=μ,X减去每一时间的平均值;(2)对均值化的矩阵进行标准化,得到标准化矩阵S;(3)基于标准化矩阵求得相关系数矩阵R为: 求式(2)中Z1=XU1,也就是求特征向量U1=(u11,u21…uN1)
,使U1U1′=I时,Var(Z1)达到最大;等价于求相关矩阵R的特征值和特征向量。设λ=λ1是R的最大特征值,对应的特征向量即为所求U1。X的第i主成分得分Zi可通过R的第i特征向量Ui得到,即: 式中,X是M×N原始矩阵;U是N×N特征向量矩阵,也称为载荷矩阵;Zi是新变量得分矩阵,称为主成分得分或者分量,代表原始数据X的线性组合,图 2为TPCA分析过程中的矩阵。
-
TPCA解算可以得到特征值、特征向量、主成分得分3个变量,其意义分别为:(1)特征值是矩阵的一种数学特性,用于决定提取分量的个数,测量主成分分量所解释的方差大小。特征值越大,主成分分量的方差越大,分量解释的信息占比(即方差贡献率=主成分特征值/特征值总和)越多。(2)特征向量是特征值对应的向量。时间主成分分析的特征向量可以表征主成分的时间特征,从而揭示可能存在的形变规律。(3)依据方差大小选择特征向量即主分量方向,将原始数据投影到主分量方向上,得到主成分得分。表征每个PS点的主要信息,这里主要是指PS点的空间分布特征。经TPCA分解出的特征向量和主成分得分能够反映数据变化的时间空间结构,按照特征值大小排列,量化可解释原信息的大小。从定义上来看,PCA目的是线性组合找到新变量代替原变量即时间变量,因而TPCA的新变量侧重描述空间特征。
1.1. TPCA原理
1.2. TPCA的结果
-
本文以北京平原区为研究对象,技术路线如下:(1)借助Sarproz平台,对2003—2010年的Envisat ASAR数据及辅助数据数字高程模型(digital elevation model,DEM)进行PS-InSAR处理,经过预处理、生成干涉图、解算形变量,获取北京平原区地面沉降信息,得到约10万个PS点及其时序沉降结果。(2)将沉降结果与水准测量结果进行对比验证。(3)精度验证合格后,利用TPCA方法分析并提取北京平原区地面沉降的时空演化特征,如图3所示。
-
从Envisat ASAR数据集的地面沉降监测结果中提取水准点附近的监测点,将水准点附近的InSAR地面沉降数据与2005—2010年的水准点监测数据进行对比验证(见图 4),发现两者年均沉降速率的相关系数达到0.97。以水准监测值为基准,PS-InSAR监测值相对于水准监测值的平均绝对值误差为3.5 mm/a,在轻微沉降区(年均沉降量小于10 mm/a),二者相差在2.5 mm/a以内;在一般沉降区和严重沉降区(年均沉降量大于10 mm/a),二者相差10 mm/a以内。由此可知,2003—2010年的沉降结果精度相对可靠。
2.1. 技术路线
2.2. 地面沉降信息精度验证
-
对10万个永久散射体的时序累计沉降数据进行时间主成分分析,得到2003—2010年北京平原区地面沉降的主成分特征值(见表 1)和主成分特征向量(见图 5)。除了第1、第2、第3成分,其他各成分的方差贡献率不到3%,因此本文取前3个成分作为主成分进行研究。
成分 特征值 解释占比/% 累计占比/% PC1 43.95 86.18 86.18 PC2 4.42 8.66 94.84 PC3 1.21 2.37 97.21 Table 1. Principal Component Eigenvalues of Land Subsidence
特征值越大,对应主成分分量解释的信息占比越多。由表 1可知,PC1特征值的方差贡献率达到86.18%,提取了原数据的绝大部分信息和变化特征,说明PC1代表北京平原区的主要地面沉降变化特征。
由于主成分Zj=u1jX1+u2jX2+…+uNjXN,cov(Xi,Zj)=cov(ui1Z1+ui2Z2+…+uiNZN,Zj)=uijλj,由此可得Xi(所有原变量)与Zj(某一主成分得分)的相关系数为:
式中,σi为Xi方差;
为Zj方差。 可见,第j个主成分得分Zj与变量Xi的相关程度与对应主成分线性组合系数(特征向量)的大小有关,特征向量代表了主成分得分与原始变量的相关程度。
特征向量表示线性组合后的主成分得分与原始变量的相关程度,也表示主成分特征的时间变化趋势。图 5中PC1的特征向量在0.15附近保持稳定,说明地面沉降的发展趋势在该段时间内保持一致。PC2和PC3的特征向量整体变化幅度相对PC1要大一些,PC2和PC3主成分通过计算可以揭示地面沉降的部分季节性变化特征。
2005年之前,前3个特征向量的系数有明显的变化幅度,这是因为2003—2005年的采集时间稀疏,处理数据是累计沉降数据,前段时间的差异变化相对大,后期细微变化表现得不明显,导致前段时间的系数变化较大。
-
特征向量代表最大的变化方向。将标准化数据投影到特征向量上,得到的投影值是新变量,即主成分得分。不同主成分得分代表不同的空间分布特征。PC1得分揭示地面沉降的空间变化主趋势。由图 6所知,PC1得分(图 6(b))与年平均沉降速率(图 6(a))空间分布一致,均体现出北京平原区的主要沉降漏斗和不均匀沉降,且相关系数达0.98(图 6(c)),PC1方向上,地面沉降速率分区明显(图 6(d));PC1特征向量的系数差异不大,变化趋势稳定;PC1的方差贡献率最大,代表最大的变化特征。PC1表示地面沉降空间特征在2003—2010年一直发展,不均匀沉降显著。
-
PC2与可压缩层厚度相关性较强。观察PC2得分的空间分布(图 7(b)),黑色虚线范围为-30 mm/a沉降漏斗区(下文同),发现平原区北部沉降区和东部漏斗沉降区存在不同的空间特征。PC2得分的正值分布范围与130 m以深的可压缩层空间分布相似(图 7(a));且两者相关系数达0.61,相关性较大(图 7(c));PC1和PC2得分的散点图中,PC2方向上,较厚(> 130 m)与较薄(60~130 m)的PS点有较为明显的分区现象(图 7(d)),说明可压缩层厚度在PC2方向上差异明显。因此PC2与可压缩层厚度相关性较强。
-
PC1得分绝对值远大于PC2,将两者叠加会使PC2的空间特征被PC1的空间主趋势掩盖。通过比较两分量的相位关系,研究两者的变化关系,发现PC1得分为负、且PC2为正的PS点覆盖了-30 mm/a沉降漏斗范围的严重沉降区域,如图 8(a)的Ⅱ象限。其余两分量的关系分别为:图 8(a)的I象限代表PC1、PC2均为正值,覆盖了非沉降区域;III象限代表PC1、PC2均为负值,覆盖一部分沉降区域和非沉降区域;IV象限代表PC1为正值,PC2为负值,覆盖一部分非沉降区域。
由于严重沉降区域对经济社会生活有重大影响,因此着重分析严重沉降区域PS点的规律。通过空间自相关分析发现严重沉降区域的分类情况。对图 8(a)II象限的监测点的PC1得分进行空间自相关分析,得到该PS点的聚集分类图,如图图 8(b)所示。研究区的北部沉降区域(海淀苏家坨、昌平沙河-八仙庄、顺义)呈现高高集聚(见图 8(b)橙色部分),南部沉降区域(朝阳来广营、金盏-东坝、黑庄户)呈现低低集聚(见图 8(b)绿色部分)。
-
研究区北部和南部沉降区域体现出不同的季节性沉降。通过统计一年内前后半年的累计沉降量来分析季节性沉降,将作差之后的沉降量(即春夏季节沉降量-秋冬沉降量)称为半年差分沉降量。半年差分沉降量为负值,则春夏季节沉降量多于秋冬季节沉降量;为正值,则秋冬季节沉降量多于春夏季节沉降量。从半年差分沉降量(见图 9)的统计中发现,北部沉降区域的春夏季节沉降量较多,南部沉降区域的秋冬季节沉降量较多。严重沉降区域的高高集聚、低低集聚分类图(见图 8(b))与严重沉降区域北部、南部季节性沉降的空间分布相似(见图 9)。
进一步分析研究区北部沉降区和南部沉降区的差异。沿北部沉降区域的红虚线AA
、南部沉降区域的黑实线BB (见图 9(a))做剖面分析,得到沿线2004—2009年的半年地面沉降量剖面图(见图10)。由图 9和图 10发现:(1)北部沉降区域整体上在春夏季节的沉降更多。2005年数据缺失严重,未体现春夏季节的沉降多;2006年、2008年、2009年均有部分季节性形变不明显的地方,其沉降速率相较前一年变化大(图 11(a)),导致季节性不明显,沉降速率变化大可能与该位置地下水的水位异常变化有关系。(2)南部沉降区域整体上秋冬季节沉降更多。2005年同样未体现秋冬季节的沉降多;2004年、2006年、2008年、2009年,范围约100~8 000 m部分对应来广营沉降漏斗,该沉降漏斗区别于其他两个沉降漏斗,可能与其地下水开采有较大差异有关;2007年相对前一年的沉降速率变化最大(见图 11(b)),其季节性沉降不明显。
3.1. TPCA特征值解释原始信息大小
3.2. TPCA主成分得分揭示地面沉降空间变化特征
3.2.1. PC1得分分析
3.2.2. 与可压缩层厚度空间分布相似的PC2
3.3. 主成分关系及季节性形变
3.3.1. 基于主成分关系对沉降区域PS点进行分类
3.3.2. 沉降区域的季节性特征
-
本文利用TPCA对地面沉降数据进行分解,PS点作为观测点,时间作为变量。TPCA实质上是从时间维度上进行降维,将具有相同空间变化特征的观测点分离,线性组合的新变量代替原始的时间变量,得到空间特征。得到以下结论:(1)TPCA应用于具有空间位置和时序属性的数据集,可以识别出该空间下的主要特征。特征值决定了各分量解释的信息量大小。本文主成分特征向量为一个时间序列,代表PC空间模式与沉降之间的相关关系。(2)主成分得分是TPCA分解得到的空间模式,表示的是不同的空间特征;本文的PC1得分与年平均沉降速率空间分布一致,揭示地面沉降的空间变化主趋势;PC2得分与可压缩层厚度相关性较强。(3)从主成分得分关系中得出严重沉降区域的分类情况,主要体现南北沉降区域的不同聚类现象,该现象与季节性沉降相关。同时研究区的北部、南部沉降区域存在较明显的季节性差异,与前面所描述的集聚分类空间分布相似。
TPCA属于数据驱动方法,不需要先验模型做出假设,主要利用多元统计原理,仅根据数据内在结构提取具有统计意义的区域信息。所以相对假设驱动方法具有更好的自适应性和信号发掘能力,具有比较强的客观性,客观揭示数据内在信息。缺陷在于提取的分量不一定与现实吻合,缺乏一定的现实意义。由于TPCA是线性组合,通过找到方差最大的方向进行投影,得到的变量仅仅是不相关,并不相互独立。PCA只用到了原始数据的二阶统计信息,而忽略了其高阶统计信息。因此,需要借助旋转主成分进行优化,找到更多主成分蕴含的物理意义。