Temporal and Spatial Analysis of Land Subsidence in Beijing Plain Based on TPCA

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个监测数据。

Figure 1.  Data Structure of TPCA

X是M个样本、N个变量构成的原始矩阵，其表达式为：

第i个变量可表示为X_i=（x_1i，x_2i…x_Mi）^T，通过TPCA对X线性变换后，得到新变量矩阵Z：

其中，

为了用Z₁代替原来的N个变量X₁，X₂…X_N，要求Z₁尽可能多地反映原来N个变量的信息，信息量大小用方差来描述。方差Var（Z₁）越大，表示Z₁包含的信息越多。

当变换矩阵满足UU$\mathrm{\text{'}}$=I，存在满足条件的U₁使得Var（Z₁）达到最大，Z₁就称为第一主成分（或者主分量）；使得方差Var（Z₂）达到第二最大，则Z₂成为第二主成分，以此类推求得剩余成分。具体过程如下：（1）对矩阵X均值化：求矩阵X的均值E‍（X）‍=μ，X减去每一时间的平均值；（2）对均值化的矩阵进行标准化，得到标准化矩阵S；（3）基于标准化矩阵求得相关系数矩阵R为：

求式（2）中Z₁=XU₁，也就是求特征向量U_1‍=（u₁₁，u₂₁…u_N1）$\mathrm{\text{'}}$，使U₁U₁′=I时，Var（Z₁）达到最大；等价于求相关矩阵R的特征值和特征向量。设λ=λ₁是R的最大特征值，对应的特征向量即为所求U₁。X的第i主成分得分Z_i可通过R的第i特征向量U_i得到，即：

式中，X是M×N原始矩阵；U是N×N特征向量矩阵，也称为载荷矩阵；Z_i是新变量得分矩阵，称为主成分得分或者分量，代表原始数据X的线性组合，图 2为TPCA分析过程中的矩阵。

Figure 2.  Matrix from TPCA

TPCA解算可以得到特征值、特征向量、主成分得分3个变量，其意义分别为：（1）特征值是矩阵的一种数学特性，用于决定提取分量的个数，测量主成分分量所解释的方差大小。特征值越大，主成分分量的方差越大，分量解释的信息占比（即方差贡献率=主成分特征值/特征值总和）越多。（2）特征向量是特征值对应的向量。时间主成分分析的特征向量可以表征主成分的时间特征，从而揭示可能存在的形变规律。（3）依据方差大小选择特征向量即主分量方向，将原始数据投影到主分量方向上，得到主成分得分。表征每个PS点的主要信息，这里主要是指PS点的空间分布特征。经TPCA分解出的特征向量和主成分得分能够反映数据变化的时间空间结构，按照特征值大小排列，量化可解释原信息的大小。从定义上来看，PCA目的是线性组合找到新变量代替原变量即时间变量，因而TPCA的新变量侧重描述空间特征。

本文以北京平原区为研究对象，技术路线如下：（1）借助Sarproz平台，对2003—2010年的Envisat ASAR数据及辅助数据数字高程模型（digital elevation model，DEM）进行PS-InSAR处理，经过预处理、生成干涉图、解算形变量，获取北京平原区地面沉降信息，得到约10万个PS点及其时序沉降结果。（2）将沉降结果与水准测量结果进行对比验证。（3）精度验证合格后，利用TPCA方法分析并提取北京平原区地面沉降的时空演化特征，如图3所示。

Figure 3.  Technical Route to Analyze Land Subsidence Characteristics Using TPCA

从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年的沉降结果精度相对可靠。

Figure 4.  Accuracy Verification of Land Subsidence and Lev‍eling Data

对10万个永久散射体的时序累计沉降数据进行时间主成分分析，得到2003—2010年北京平原区地面沉降的主成分特征值（见表 1）和主成分特征向量（见图 5）。除了第1、第2、第3成分，其他各成分的方差贡献率不到3%，因此本文取前3个成分作为主成分进行研究。

Figure 5.  Time Series Loading Diagram of PC1-‍PC3

特征值越大，对应主成分分量解释的信息占比越多。由表 1可知，PC1特征值的方差贡献率达到86.18%，提取了原数据的绝大部分信息和变化特征，说明PC1代表北京平原区的主要地面沉降变化特征。

由于主成分Z_j=u_1jX₁+u_2jX₂+…+u_NjX_N，cov‍（X_i，Z_j）=cov（u_i1Z₁+u_i2Z₂+…+u_iNZ_N，Z_j）=u_ijλ_j，由此可得X_i（所有原变量）与Z_j（某一主成分得分）的相关系数为：

式中，σ_i为X_i方差；$\sqrt[]{{\lambda }_j}$为Z_j方差。

可见，第j个主成分得分Z_j与变量X_i的相关程度与对应主成分线性组合系数（特征向量）的大小有关，特征向量代表了主成分得分与原始变量的相关程度。

特征向量表示线性组合后的主成分得分与原始变量的相关程度，也表示主成分特征的时间变化趋势。图 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年一直发展，不均匀沉降显著。

Figure 6.  Relationship Between PC1 Score and Average Annual Subsidence Rate

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与可压缩层厚度相关性较强。

Figure 7.  Relationship Between PC2 Component and Compressible Thickness(> 130 m)

PC1得分绝对值远大于PC2，将两者叠加会使PC2的空间特征被PC1的空间主趋势掩盖。通过比较两分量的相位关系，研究两者的变化关系，发现PC1得分为负、且PC2为正的PS点覆盖了-30 mm/a沉降漏斗范围的严重沉降区域，如图 8（a）的Ⅱ象限。其余两分量的关系分别为：图 8（a）的I象限代表PC1、PC2均为正值，覆盖了非沉降区域；III象限代表PC1、PC2均为负值，覆盖一部分沉降区域和非沉降区域；IV象限代表PC1为正值，PC2为负值，覆盖一部分非沉降区域。

Figure 8.  Relationship Between Two Principal Components

由于严重沉降区域对经济社会生活有重大影响，因此着重分析严重沉降区域PS点的规律。通过空间自相关分析发现严重沉降区域的分类情况。对图 8（a）II象限的监测点的PC1得分进行空间自相关分析，得到该PS点的聚集分类图，如图‍图 8（b）所示。研究区的北部沉降区域（海淀苏家坨、昌平沙河-八仙庄、顺义）呈现高高集聚（见图 8（b）橙色部分），南部沉降区域（朝阳来广营、金盏-东坝、黑庄户）呈现低低集聚（见图 8（b）绿色部分）。

研究区北部和南部沉降区域体现出不同的季节性沉降。通过统计一年内前后半年的累计沉降量来分析季节性沉降，将作差之后的沉降量（即春夏季节沉降量-秋冬沉降量）称为半年差分沉降量。半年差分沉降量为负值，则春夏季节沉降量多于秋冬季节沉降量；为正值，则秋冬季节沉降量多于春夏季节沉降量。从半年差分沉降量（见图 9）的统计中发现，北部沉降区域的春夏季节沉降量较多，南部沉降区域的秋冬季节沉降量较多。严重沉降区域的高高集聚、低低集聚分类图（见图 8（b））与严重沉降区域北部、南部季节性沉降的空间分布相似（见图 9）。

Figure 9.  Seasonal Variation of Cumulative Land Subsidence in 2004 and 2009

进一步分析研究区北部沉降区和南部沉降区的差异。沿北部沉降区域的红虚线AA$\mathrm{\text{'}}$、南部沉降区域的黑实线BB$\mathrm{\text{'}}$（见图 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）），其季节性沉降不明显。

Figure 10.  Seasonal Land Subsidence

Figure 11.  Subsidence Rate Profile Over the Years

本文利用TPCA对地面沉降数据进行分解，PS点作为观测点，时间作为变量。TPCA实质上是从时间维度上进行降维，将具有相同空间变化特征的观测点分离，线性组合的新变量代替原始的时间变量，得到空间特征。得到以下结论：（1）‍TPCA应用于具有空间位置和时序属性的数据集，可以识别出该空间下的主要特征。特征值决定了各分量解释的信息量大小。本文主成分特征向量为一个时间序列，代表PC空间模式与沉降之间的相关关系。（2）主成分得分是TPCA分解得到的空间模式，表示的是不同的空间特征；本文的PC1得分与年平均沉降速率空间分布一致，揭示地面沉降的空间变化主趋势；PC2得分与可压缩层厚度相关性较强。（3）从主成分得分关系中得出严重沉降区域的分类情况，主要体现南北沉降区域的不同聚类现象，该现象与季节性沉降相关。同时研究区的北部、南部沉降区域存在较明显的季节性差异，与前面所描述的集聚分类空间分布相似。

TPCA属于数据驱动方法，不需要先验模型做出假设，主要利用多元统计原理，仅根据数据内在结构提取具有统计意义的区域信息。所以相对假设驱动方法具有更好的自适应性和信号发掘能力，具有比较强的客观性，客观揭示数据内在信息。缺陷在于提取的分量不一定与现实吻合，缺乏一定的现实意义。由于TPCA是线性组合，通过找到方差最大的方向进行投影，得到的变量仅仅是不相关，并不相互独立。PCA只用到了原始数据的二阶统计信息，而忽略了其高阶统计信息。因此，需要借助旋转主成分进行优化，找到更多主成分蕴含的物理意义。