2020, Vol. 63
2. 武汉大学卫星导航定位技术研究中心, 武汉 400079
2. GNSS Research Center, Wuhan University, Wuhan 400079, China
地壳受到冰川均衡调整、水文负载与卸载、构造运动等作用而发生垂向形变，这些地壳垂向形变在连续空间和时间域内呈现了显著特征.过去二十几年，国内外学者采用GNSS、机载数据、InSAR等先进的测量技术进行了大量监测，实现了单站“点模式”和“区域影像模式”的形变监测，获取了时空域内地壳垂向形变的特征(Blewitt et al., 2010)；揭示地壳形变的时空变化特征将有助于理解地球物理过程，为研究地球内部相互作用机制提供技术支持.
在过去的研究中，GNSS技术是典型的单站“点模式”的形变监测，它主要通过利用时间序列建立速度模型实现.GNSS时间序列的形变监测兼具地表位移信息的精确性及时间域内连续性等优点，已成为研究地壳垂向形变随时间演变过程及其物理机制重要的基础数据源.然而，单站“点模式”受制于系统误差、构造因素、非构造因素与人为因素，不利于提取研究区域内地壳垂向运动的共性空间特征，导致其仍是大地测量领域内极具挑战的难题.
InSAR技术是实现“区域影像模式”形变监测的主要手段，具有全天候成像、空间分辨率高、覆盖范围广等优势，可实现区域形变整体特征提取(谭凯等，2016；单新建等，2015；许才军等，2015；蒋廷臣等，2013；Liu et al., 2011；Jiang et al., 2015).但是，InSAR容易受时变去相的影响，难以提供足够的时间分辨率.即便卫星硬件的不断发展使得重返周期日益缩短，但微小形变与噪声之间的混叠，导致噪声更加复杂(朱建军等，2017).因此，获取大范围、高时空分辨率的地壳垂向形变图仍是一项具有挑战的难题.
近20年来，在观测站数量日益增多、数据不断积累、覆盖范围扩大的背景下，充分利用研究区域内GNSS测站的共性特点，综合地质统计方法和图像处理技术，对GNSS时间序列速度场进行空间加密，保留变形区域的边界特征，更加细致全面地描绘地壳垂向形变在空间上的分布，由此构建反映研究区域内连续空间地壳垂向形变的图像，即GNSS影像(Hammond et al., 2016)，使得兼顾时间和空间分辨率的地壳垂向形变研究成为可能.
1 GNSS成像原理GNSS成像方法(Hammond et al., 2016)基于区域内多个测站时间序列的速度及其不确定性，通过空间结构函数描述区域内测站速度间的相关性，进而构建顾及空间结构的滤波器，剔除粗差、增强区域共性信号，并借鉴图像处理技术对速度场进行空间加密，从而实现地壳垂向形变的时空特征提取.GNSS成像方法既可保留变形区域的边界特征，又能描绘垂向形变在连续空间内随时间的演变特征.
1.1 速度及其不确定性估计受构造运动、非构造运动及人为因素影响，GNSS坐标时间序列呈现周期性谐波信号、粗差、阶跃以及数据的缺值等现象.常用的速度及不确定性估计方法是基于最小二乘准则的参数估计，该估计方法易受上述因素影响，对序列的数据预处理要求甚高，且难以对地壳垂向形变的速度进行准确估计.针对上述难题，本文采用基于中位数并顾及年际差异的非参数方法(Median Interannual Difference Adjusted for Skewness，MIDAS方法)(Blewitt et al., 2016)实现速度及不确定性的估计.
MIDAS方法是改进的泰尔森回归(Theil-Sen Regression)方法(Theil，1950；Sen，1968)，它利用绝对中位差与总体标准差消除粗差(Wilcox，2013)，并根据周年窗口选取数据对，消除周期性谐波信号的影响.其具体实现步骤：
(1) 数据对选取与速度计算.对GNSS坐标时间序列正序与逆序选取出所有时间间隔为1年的数据对，并计算其速度vk：
|
(1) |
其中，xi与xj是GNSS测站在ti与tj历元的垂向坐标.计算得到的所有速度表示为：θ:{v1, v2, …, vn}.
(2) 计算绝对中位差MAD，并估计总体标准差σ：
|
(2) |
|
(3) |
其中，vi表示由各数据对计算的速度值，
(3) 本文设置2σ作为阈值条件剔除粗差，即剔除不满足以下条件的值：
|
(4) |
此时，再次计算获得的中位数
(4) 计算有效历元个数(公式(5))并估计速度估值的标准差
|
(5) |
|
(6) |
|
(7) |
其中，Npairs表示步骤(1)中选取的数据对个数，N表示有效历元个数, σ表示剔除粗差后的总体标准差.
1.2 空间结构函数的构建空间结构函数(Spatial Structure Function，SSF)(Hammond et al., 2016)根据GNSS测站间距，将研究区域内测站对分组，并根据站间速度差的统计结果，构建站间距离与速度相关性之间的映射关系，反映区域内测站间速度的相关性.具体过程如下：
(1) 研究区域内任意两测站构成测站对，计算每组测站对速度差的绝对值δvij(i=1, 2, …, n; j=1, 2, …, n)，及其相应的测站对间距离Δij，并按站间距离从小到大排序.本文使用测站对所在大圆弧的圆心角作为测站间距离，速度差绝对值按照公式(8)计算：
|
(8) |
式中，vi与vj分别代表第i个测站与第j个的速度；
(2) 依据站间距离，划分出若干个数据池(通常以0.25log10Δ作为数据池的分隔距离，以保证各数据池中至少含有20个测站对为准则)，并计算各数据池内测站速度差绝对值δv的绝对中位差(MAD)：
|
(9) |
式中，i=1, 2, …, n; k=1, 2, …, m，MADk表示m个数据池中第k个数据池的绝对中位差，δvi表示该数据池内对应的n个速度差绝对值，δ
(3) 构建并标准化空间结构函数：
|
(10) |
|
(11) |
式中，ssfi指的是第i个数据池中点的空间结构函数值，MADi指的是在第i个数据池中计算得到的绝对中位差.
PBO(Plate Boundary Observation)站网的空间结构函数随站间距离变化以分段函数的形式呈现(如图 1所示).随着站间距离增大，空间结构函数的值单调递减，并在站间距离达到21.4 km左右时急剧减小为0，较好地描述了地壳运动在空间分布上的相关性.因此，空间结构函数可强化研究范围内对多个测站的共性描述，也有助于粗差剔除.
|
图 1 PBO GNSS网空间结构函数 Fig. 1 Spatial Structure Function of PBO GNSS network |
通过上述分析可知，具有相似空间结构的测站速度及其不确定性存在共性，本文通过顾及空间结构的中位数滤波(Median Spatial Filter，MSF)(Hammond et al., 2016)对各测站的速度值进行筛选，以增强这些区域GNSS测站垂向速度的共性.主要步骤包括：
(1) 建立测站间空间关系.本文主要通过构建球面Delaunay三角网(Renka，1997)来建立测站间的空间关系；
(2) 选定某一测站作为“滤波测站”，根据公式(12)，同时顾及速度不确定性与空间相关性，计算在三角网中与“滤波测站”相连各测站的权值并标准化(公式(13))：
|
(12) |
|
(13) |
其中，ssfi表示输入当前测站与评估测站间距离得到的空间结构函数值，σi指各测站估计得到的速度的不确定性，ωi与ω′i分别代表标准化前后的测站权值；
(3) 根据各测站权值与速度估值，计算其加权中位数(Cormen et al., 2001)，即为“滤波测站”滤波后的速度值.
1.4 格网点生成与插值为了将离散的GNSS观测站网进行加密，本文根据研究区域的范围，首先对研究区域格网化；随后，将各格网点作为1.3节中的“滤波测站”，根据其周围若干测站的速度估值及不确定性，利用MSF插值获得各格网点速度；最终生成GNSS影像.滤波前后的测站速度估值及其不确定性都能够获取GNSS影像结果；然而，滤波前的测站网速度中含有的粗差会干扰地壳垂向形变共性特征的提取，本文3.2.1节与3.2.2节对比分析了滤波前后GNSS影像结果的差异，并说明了MSF对GNSS影像生成的重要性.
1.5 棋盘格检测原理影像可靠性通常利用棋盘格检测进行验证.棋盘格检测对比分析棋盘格状虚拟速度场与GNSS影像反映的速度场，从定性角度验证GNSS成像方法的合理性与正确性，探寻GNSS影像正确反映地壳形变特征的能力，是GNSS成像方法中必不可少的一个环节.本文不仅使用棋盘格检测(Checkerboard Test)进行定性分析(Hammond et al., 2016)，还进一步提出了基于混淆矩阵总体精度(孙家抦，2009)的定量检验方法.其实现过程如下：
(1) 按照一定的间隔将研究区域划分为格网并交替填入给定绝对值的速度，本文以+3 mm·a-1与-3 mm·a-1作为输入速度，获得一个具有突变边缘的棋盘格状虚拟速度场；
(2) 对获得的虚拟速度场进行高斯平滑，使其成为棋盘格边缘速度为0 mm·a-1，内部为±3 mm·a-1，并在各棋盘格边缘与内部之间连续平缓变化的速度场；
(3) 在GNSS测站位置处取棋盘格内对应的速度值，获得采样结果图；
(4) 针对实测数据，将MIDAS方法估计的速度不确定性作为采样结果的不确定性；
(5) 以本文提出的总体精度作为量化指标，对比分析根据采样结果生成GNSS影像与虚拟速度场的差异，分析方法在下文中进行介绍.
目前，GNSS影像的可靠性判断仍以棋盘格检测的定性方法为主.为更好地表征结果的可靠性，本文提出了基于混淆矩阵的定量分析方法.首先计算虚拟速度场与GNSS影像各格网点的速度差值，设定与不确定性统计量一致的阈值，并比较速度差值与阈值，判定GNSS影像中各格网点是否正确地反映了虚拟速度场中的速度.若小于阈值，则认定该格网点正确反映了虚拟速度场的速度，否则认定该点中GNSS影像未能正确反映虚拟速度场的速度，统计正确的成像格网数量，并基于混淆矩阵的定量分析方法计算其总体精度pi：
|
(14) |
其中，nc指的是成像前后速度差值在阈值范围内的格网点数量，ntotal指的是总格网点数，0≤pi≤1.因此，总体精度的值能够反映GNSS影像还原虚拟速度场的总体情况，能够为分析GNSS影像的可靠性提供定量依据.
2 数据来源与预处理美国西部地区与中国大陆地壳活动复杂，是研究地壳垂向形变的天然实验场，两个区域地壳垂向形变在空间上的分布特征具有研究价值.因此本文分别选取美国西部地区(109°W—125°W，32°N—49°N)以及中国大陆(73°E—134°E，18°N—54°N)作为研究区域，使用板块边界观测网络(Plate Boundary Observation，PBO)与中国大陆构造环境监测网络(Crustal Movement Observation Network of China，CMONOC)两个测站网的GNSS坐标时间序列产品生成GNSS影像，以揭示这些典型区域的形变特征.
2.1 美国西部地区PBO GNSS站网本文从UNAVCO发布的GNSS数据产品(ftp://data-out.unavco.org/pub/products/position)中获取了PBO位于美国西部地区1209个测站的坐标时间序列.统计发现，有2个测站的时间序列跨度小于1年，50个测站的序列跨度介于1~5年之间，185个测站的序列跨度介于5~10年之间，673个测站的序列跨度介于10~15年之间，另有299个大于15年；此外，这些测站大多分布在地壳活跃的地区，提供了足够的数据，可用以生成GNSS影像.
在数据预处理工作中，我们将序列长度小于2.5年的测站以及可用数据率AR(Available Rate)小于0.9的测站进行剔除，剩余1132个测站，其测站分布如图 2所示.
|
图 2 PBO测站分布图 Fig. 2 Distribution of selected PBO stations |
其中，可用数据率AR的计算公式如下：
|
(15) |
其中，n有效为时间序列中的实际数据历元数，n总为时间序列中应有的数据历元数.
2.2 中国大陆CMONOC GNSS站网本文从中国地震局发布的GNSS数据产品(ftp://ftp.cgps.ac.cn/products/position)中获得了CMONOC中290个连续观测参考站的GNSS时间序列.使用与2.1节相同的标准进行数据预处理后，剩余241个测站，所有测站的时间序列跨度均大于5年，并且有28个测站的时间序列长度大于15年，能够满足速度与不确定性估计对时间序列跨度的要求，测站分布如图 3所示.
|
图 3 CMONOC测站分布图 Fig. 3 Distribution of CMONOC stations |
MSF作为空间滤波，其主要作用是剔除整网中与整体趋势不符的峰值点或突变点.为了验证MSF的可行性，3.1节对比了测站网中部分测站MSF滤波后的速度与原始速度.由于MSF可能删除部分信息，为了进一步判断GNSS影像的可靠性，3.2节利用棋盘格检测，检验GNSS影像对研究区域最小形变的空间响应.最后，生成GNSS影像，总结研究区域内地壳垂向形变的时空特征，并与已有的大地测量学、地球动力学等结果进行交叉验证.
3.1 MSF验证实验结果为了验证MSF的效果，本文采取以下步骤：(1)分别将上述两个测站网中的部分测站移除，称作“移除站网点”；(2)在未进行滤波的情况下，采用剩余的测站，通过MSF对“移除站网点”的速度进行插值恢复；(3)对比分析“移除站网点”测站速度原始与插值结果.
由于PBO的测站数量多，分布密集，因此本文将研究范围内所有测站根据经度依次排列，并将排在偶数个的测站移除，共移除566个测站，“移除站网点”分布如图 4a所示.根据MIDAS方法估计得到的测站速度及其不确定性，对“移除站网点”的速度进行插值，统计原始与插值速度场差值，判定差值是否在原始速度的3倍不确定性范围内，结果如图 4b所示(按照纬度排列).纬度36°N附近部分测站的速度不确定性不仅大于其他测站，其速度的还原效果也不理想.然而，仍有57.1%的插值速度在速度估值的1倍不确定性范围内，82.9%的插值结果在2倍不确定性范围以内，91.6%在3倍不确定性范围内.插值结果的均方根误差为3.89 mm·a-1；由于MSF具有滤波特性，为了获得更加普遍的结论，删去速度绝对值明显高于周边测站的9个测站，均方根误差迅速减小为1.27 mm·a-1.因此，本文认为在密集的测站网中，MSF较为准确地恢复了绝大部分“移除站网点”的速度，能够用于进行GNSS影像的生成.
|
图 4 PBO“移除站网点”分布(a)与插值结果(b) Fig. 4 Distribution of extracted sites (a) and interpolations (b) in PBO |
CMONOC的测站数量较少，范围广，且分布不均，任意移除点位可能会导致失去某个区域的地壳形变速度的信息.为了避免这种情况的发生，本文使用CMONOC测站进行MSF验证实验时，选取了测站较密集区域的点位进行移除与插值，其范围为经度95°E—122°E，纬度20°N—41°N.同样，将测站较密集区域所有测站根据经度依次排列，并将排在偶数个的测站移除，共移除79个测站，“移除站网点”分布如图 5a所示.统计原始与插值速度场差值，判定差值是否在原始速度的3倍不确定性范围内，结果如图 5b所示(按照经度排列).对插值结果进行统计，46.8%的插值结果在速度估值的1倍不确定性，72.1%的插值结果在2倍不确定性以内，88.6%在3倍不确定性范围内.然而，其均方根误差达到6.56 mm·a-1，观察各测站速度后，发现HECX与TJWQ两个测站的速度绝对值明显高于周边测站，因此认为该两个测站可能存在速度的峰值或粗差，MSF无法有效地插值.为了获得更加普遍的结论，删去这两个测站重新计算均方根误差；此时，均方根误差迅速下降为1.57 mm·a-1.因此本文认为，在较为稀疏的测站网中，MSF也能反映地壳形变的总体趋势，可以用于GNSS影像的生成；然而由于MSF属于空间滤波器，无法自动判断区域内的正常速度峰值或突变，在实际运用中尚需人工干预，具体描述见3.2.2节.
|
图 5 CMONOC“移除站网点”分布(a)与插值结果(b) Fig. 5 Distribution of extracted sites (a) and interpolations (b) in CMONOC |
造成CMONOC与PBO插值结果间差异的原因包括：(1)CMONOC测站坐标时间序列中粗差、阶跃与缺值情况严重；(2)CMONOC的测站密度远小于PBO，对于中国境内复杂的地壳活动细节描述仍需加密数据分布.尽管如此，仍有88.6%的插值结果在3倍不确定性以内，对于变化不显著区域，稀疏的测站分布仍可以对研究区域运动趋势进行整体研判，但细节刻画上还需进一步优化算法.
3.2 棋盘格检测结果本节主要通过生成虚拟速度场，开展棋盘格检测，以进一步验证GNSS成像算法的正确性.首先，对棋盘格进行填充与滤波，生成虚拟速度场；随后，在PBO和CMONOC站网测站位置进行采样；使用采样的虚拟速度生成GNSS影像；最后，通过定性及定量分析评估GNSS成像方法的效果.对比原始虚拟速度场与GNSS影像结果，不仅能够验证MSF对无测站分布位置速度插值结果的准确程度，还能够检验GNSS影像对研究区域内最小形变的空间响应.
3.2.1 棋盘格检测定性分析PBO的测站数量大，分布密集，本文将研究区域的划分间隔设定为3°×3°，交替填入±3 mm·a-1并进行高斯平滑，棋盘格状虚拟速度场如图 6a所示.测站速度采样结果如图 6b所示，对应生成的GNSS影像如图 7所示.定性地对比虚拟速度场与GNSS影像，大部分区域都较好地还原了原始虚拟速度场中各3°×3°格网中的速度，其中以测站密度最大的美国西部海岸沿线地区与加利福尼亚州地区最佳，有效地反映了棋盘格之间速度的变化趋势，很好地恢复了3°×3°间隔的虚拟速度场；但测站分布较稀疏区域(如爱达荷州与蒙大拿州地区、亚利桑那州地区)的GNSS影像中，棋盘格不仅面积缩小，还发生了变形，边界也出现了混淆.因此，基于测站空间分布，选择尺度适中的棋盘格，对有效检测GNSS影像对形变区域的空间响应具有重要的参考价值.
|
图 6 美国西部地区虚拟速度场(a)与采样结果图(b) Fig. 6 Synthetic velocity field (a) and sampling result (b) in western US |
|
图 7 美国西部地区虚拟速度场的GNSS影像 Fig. 7 GNSS image of the synthetic velocity field in western US |
CMONOC覆盖范围广，却较为稀疏，进行棋盘格检测时，基于测站的分布情况，初步设定间隔为7°×7°.该区域高斯平滑后的棋盘格状虚拟速度场如图 8a所示.相似地，在测站的坐标位置上进行采样，结果如图 8b所示.将生成的GNSS影像(图 9)与虚拟速度场进行比较后发现，即使将划分的间隔设定为7°×7°，CMONOC的测站网对于原始虚拟速度场的还原情况在某些区域(如青藏高原地区、东北地区)仍有不足.本文认为测站的疏密程度对于GNSS影像的生成具有较大的影响，不同密度的测站网能够反映的地壳形变空间上的波长也有差异.
|
图 8 中国大陆区域虚拟速度场(a)与采样结果图(b) Fig. 8 Synthetic velocity field (a) and sampling result (b) in Chinese mainland |
|
图 9 中国大陆虚拟速度场的GNSS影像 Fig. 9 GNSS image of the synthetic velocity field in Chinese mainland |
本节进一步将中国大陆分为四个区域进行细化定性分析：在96°E—120°E的范围内(图 9中区域Ⅱ)，各棋盘格内速度的整体趋势得到了有效的反映，棋盘格的变形程度较低；96°E以西、40°N以北地区(图 9中区域Ⅳ)，虽然棋盘格的面积发生了缩小，但是仍然还原了棋盘格边缘较丰富的信息；在120°E以东、40°N以北地区(图 9中区域Ⅲ)棋盘格出现了较严重的变形，推测该区域测站网反映7°×7°范围的地壳垂向形变的能力较弱；而74°E—96°E，28°N—40°N范围内(图 9中区域Ⅰ)仅有22个测站，因此棋盘格发生了严重的缩小、变形，边界也发生混淆，几乎无法还原虚拟速度场的空间分布.因此，本文认为该区域测站仅能反映较长波长的形变.为了更加细致地描绘该区形变特征，有必要新增观测信息，或进一步研究反映较短波长形变的方法.
3.2.2 棋盘格检测定量分析对美国西部和中国大陆两区域的棋盘格检测进行定量分析，依据CMONOC和PBO测站速度不确定性的均值、中位数等统计信息，设置了大小不同的阈值，得到相应的总体精度如表 1所示.美国加利福尼亚州在阈值设置为0.59 mm·a-1时，总体精度仍然保持在90%以上；中国大陆在图 9区域Ⅱ的总体精度也能够达到70%以上；结果表明，在测站相对密集区域，成像效果能够反映研究区域的整体趋势及空间分布情况，达到预期水平.接下来分析美国西部地区与中国大陆的总体精度：阈值设定为0.77 mm·a-1时，美国西部地区的虚拟速度场根据PBO部分测站还原的总体精度达到68.32%，中国大陆的虚拟速度场根据CMONOC及周边部分测站还原的总体精度也达到60.42%；即使随着阈值的降低，美国西部地区还原虚拟速度场的总体精度仍然能够保持在60%以上，中国大陆的总体精度也能够保持在55%以上.根据上述结果，本文认为PBO在美国西部地区能够较好地反映地壳垂向形变的时空变化特征，且具有一定的可靠性；而中国大陆由于区域面积大，测站密度不足，总体精度较美国西部地区低，但是其仍能够在一定程度上反映地壳垂向形变的时空变化特征；结合定性分析结果可知，在测站密度稀疏的青藏高原、中国东北地区，GNSS影像的可靠性会下降.测站空间分布及密度对GNSS成像可靠性的影响需进一步研究.
|
表 1 不同阈值下两区域棋盘格检测的总体精度比较 Table 1 General accuracy of checkerboard test in two regions with difference thresholds |
我们对两区域内格网点在虚拟速度场与GNSS影像中的速度差值进行了统计(如图 10).速度差值反映了GNSS影像恢复虚拟速度场的效果；差值越小，说明GNSS影像越准确地还原了虚拟速度场；差值越大，则说明GNSS影像还原虚拟速度场的效果越不理想.该统计结果表明，美国西部地区速度差值普遍比较小，主要分布于0.2 mm·a-1以内，处于1倍中误差范围之内，且格网点个数随差值的增大迅速减小，也说明美国西部地区的虚拟速度场较好地得到了还原；而中国大陆的统计结果显示，其速度差值主要分布于0.4 mm·a-1以内，同样随差值增大，格网点个数减小；然而，其落在0.2 mm·a-1以内的个数远小于美国西部地区，再次说明CMONOC对原始速度场的还原情况稍差.由此推断，棋盘格检测时不仅要顾及测站的空间分布，也要考虑速度在空间上的梯度变化.
|
图 10 两区域格网点在虚拟速度场与GNSS影像中速度差值分布图 (a)美国西部地区的统计结果; (b)中国大陆区域的统计结果. Fig. 10 Distribution of the velocity difference between grid points from synthetic velocity field and GNSS image (a) The statistical result of western US; (b) The statistical result of Chinese mainland. |
GNSS测站速度及其不确定性能否反映研究区域内的地壳运动共性，是影响GNSS影像效果的关键因素之一.因此，本文分别使用测站网的原始及MSF滤波后的速度与不确定性估值，生成美国西部地区和中国大陆的GNSS影像(如图 11与图 12).
|
图 11 PBO站网的测站速度(a)与不确定性分布(b)及生成的GNSS影像(c和d，其中，c为利用测站网的原始速度及其不确定性估值生成的GNSS影像，d为利用MSF以后的速度及其不确定性生成的GNSS影像) Fig. 11 Distribution of PBO stations velocities (a) and uncertainties (b), and GNSS images of western US ((c) is the GNSS image generated with the initial estimated velocities and (d) is the GNSS image generated with the velocities after MSF) |
|
图 12 CMONOC站网的测站速度(a)与不确定性分布(b)及生成的GNSS影像(c和d，其中，c为使用测站网的原始速度及其不确定性估值生成的GNSS影像，d为使用MSF以后的速度及其不确定性生成的GNSS影像) Fig. 12 Distribution of CMONOC stations velocities (a) and uncertainties (b), and GNSS images of Chinese mainland ((a) is GNSS image generated with the initial estimated velocities and (b) is GNSS image generated with the velocities after MSF) |
通过定性比较GNSS影像发现，PBO站网MSF前生成的GNSS影像(图 11c)中具有部分的零碎图斑，其中包括区域内的速度峰值(如加利福尼亚州的海岸山脉地区)与突变值(如加利福尼亚州的中央谷地)以及速度估值中粗差在GNSS影像中的影响.然而，在PBO站网MSF后生成的GNSS影像(图 11d)中，仅在加利福尼亚州北部存在一处极为细小的图斑.因此，MSF有效地消除了速度估值粗差与区域内速度峰值、突变的影响，增强了区域的共性，达到了滤波的效果，能够体现研究区域整体的地壳垂向运动趋势，有利于GNSS影像的定性分析，是非常有必要的.
在CMONOC站网MSF前生成的GNSS影像(图 12c)中，同样具有由于速度估值的粗差、峰值与突变产生的零碎图斑.此外，由于CMONOC测站密度较小，Delaunay三角网的空圆特性(Renka，1997)造成GNSS影像中出现了许多规则圆弧形的突变边缘，会影响GNSS影像中各区域地壳垂向运动趋势的分析.然而，在CMONOC站网MSF后生成的GNSS影像(图 12d)中，除了测站稀疏的青藏高原地区与大兴安岭地区外，其规则的圆弧状边界的状况得到了很大程度的减弱.同时，原有的零碎图斑也有效地通过滤波进行了消除.因此，MSF滤波器使GNSS影像的整体趋势更加显著，更具区域代表性，有助于对不同区域地壳垂向运动空间特征进行整体评价.这说明了在稀疏的测站网中，MSF具有更加显著的效果，再次体现了MSF在GNSS影像生成中的重要性.
综上所述，MSF是在GNSS影像生成过程中不可或缺的步骤，因此下文主要选用各测站经过MSF后生成的GNSS影像与已有成果进行对比分析.
3.3.2 MSF存在的问题MSF滤波后的GNSS影像也呈现出部分问题.例如，美国西部GNSS影像结果(图 11)表明，华盛顿州的西北海岸地区的山地在经过滤波后，其速度的峰值减小，区域面积也轻微地发生了缩小；海岸山脉(Coast Ranges)在加州与俄勒冈州交界处的部分以及加州中南部与中央谷地平行的部分，其速度的峰值发生了减小，在滤波前复杂的速度变化情况被简化；中央谷地(Central Valley)本就由于测站稀疏导致反映出的沉降范围小且不连续，在经过滤波后不仅速度大小发生了减小，其峰值区域也被平滑，无法有效地对该区域的地壳垂向运动趋势进行描述；在怀俄明州、蒙大拿州与爱达荷州交界的赫布根湖(Hebgen Lake)区域，其隆升的速度与范围都由于滤波发生了缩小.以上这些现象都是由于在使用顾及空间结构的中位数滤波时，各测站原始的速度估值中的粗差与峰值、突变发生了混叠，而该滤波方法无法有效地对其进行区分，造成过度平滑现象.因此，在空间特征突变显著的重点区域需要人为干预，确定MSF滤波器的适用范围.
对中国大陆(图 12)进行分析时发现，其存在的问题仍然是对速度峰值的过度平滑现象，在青藏高原东南部、天山山脉、华北地区、长江中下游地区都有体现.此外，在测站网稀疏的西部地区、东北地区，由于Delaunay三角网的空圆特性产生的规则圆弧形边界也是在GNSS影像的结果分析造成影响的重要因素之一.
综上所述，MSF目前存在的问题包括：(1)无法有效区分粗差与区域内的速度峰值，导致速度峰值区域发生过度平滑现象；(2)会导致部分整体趋势发生扩大或缩小；(3)在测站极为稀疏的区域无法解决规则圆弧形边界的问题.因此，本文认为在进行GNSS影像结果分析时可以结合MSF前后生成的GNSS影响对地壳垂向形变的速度场进行分析，即利用MSF后GNSS影像的整体趋势，结合MSF前GNSS影像中的有用信息，对GNSS影像结果进行后续分析.针对上述问题，在今后的研究中，可以将已有的大地测量学结果与地球动力学模型作为先验信息，进一步对MSF进行优化.
3.3.3 交叉验证棋盘格定性与定量分析结果表明，本文使用美国西部地区PBO测站网构建的GNSS影像能够通过3°×3°棋盘格检测，总体精度较高，可靠地反映该地区较小规模的地壳垂向形变趋势.将美国西部GNSS影像(图 11c和d)与已有的大地测量学与地球动力学成果进行对比分析得出：海岸山脉(Coast Ranges)总体呈现出隆升的趋势，并且在华盛顿州西北部、俄勒冈州与加利福尼亚州交界处以及加利福尼亚州中南部速度最大，达到2 mm·a-1；内华达山脉(Sierra Nevada)以较大的速度隆升，其速度峰值达到3 mm·a-1，与Hammond等(2016)结果接近；加利福尼亚州东南的莫哈韦沙漠(Mojave Desert)南部以1.5 mm·a-1的速度隆升，莫哈韦沙漠北部与内华达州南部的黑岩沙漠(Black Rock Desert)较为稳定(速度为0 mm·a-1左右)，与Hammond等(2016)、Hamilton和Myers(1966)结果一致；洛杉矶(LA)地区呈现小范围的沉降趋势，速度最大达到-2.5 mm·a-1；维多利亚河(Victoria River)及其河谷地区呈现大范围的沉降，速度为-1.5 mm·a-1，与Hamilton和Myers(1966)的结果具有一致性；蒙大拿州、爱达荷州与怀俄明州交界的赫布根湖(Hebgen Lake)地区以较大速度(1.5 mm·a-1)呈现隆升运动，而其周边的蛇河平原(Snake River Plain)东部、蒙大拿州西南部呈现沉降运动，速度在-1 mm·a-1左右，该结果与Reilinger等(1979)一致.这些结果都能够和已有结果相互解释，验证了GNSS影像具有合理性与正确性.
然而，即使将棋盘格的划分间隔定为7°×7°，CMONOC在测站稀疏区域还原虚拟速度场的情况依旧不理想.因此，本文在使用已有的精密水准测量与GPS成果进行交叉验证时，一方面仅挑选了几个本文认为具有代表性的区域(张培震等，2003)；另一方面，通过较大范围的定性分析结合定量描述的方法进行对比分析.
图 12d展示的沉降区域分别呈现出：(1)柴达木、塔里木、准噶尔盆地的下降运动，速度分别为-0.6 mm·a-1、-1.2 mm·a-1、-0.9 mm·a-1，与董鸿闻等(2002)结果一致；(2)甘青地块的下沉，速度为-0.3~-0.9 mm·a-1，且西部下降较大，与精密水准测量结果(董鸿闻等，2002)相同；(3)东北地区的下沉趋势，速度为-0.6~-1.5 mm·a-1，该结果与应绍奋等(1988)、郭良迁(1990)、董鸿闻(2002)与赖锡安等(2004)都具有一致性；其中三江平原下沉速度最大，达到-1.5 mm·a-1，该结果与赖锡安等(2004)吻合；(4)长江下游地区的沉降运动以速度为-2.1 mm·a-1的苏锡常地区为中心，向外逐渐减小为-0.6 mm·a-1，与胡惠民等(1992)、王伟(2013)的研究结果一致；(5)华南地区的总体沉降趋势，速度为-0.6~-1.5 mm·a-1，与王敏(2009)的结果一致；(6)尽管天津及其以南地区的速度在图 12b中呈现为0 mm·a-1，与王伟(2013)得到的总体沉降趋势不尽相符，但结合MSF前GNSS影像(图 12c)中该区域的沉降趋势，仍可认为该区域为下沉区.
图 12d的隆升区域中，可以观察得到：(1)天山地块在新疆地区沉降的背景下隆升，呈现的速度为0.3~0.6 mm·a-1，与赖锡安等(2004)的水准测量结果一致；(2)青藏高原南部、喜马拉雅山脉的隆升运动，隆升速度较大，达到1.0 mm·a-1，总体趋势与董鸿闻等(2002)、赖锡安等(2004)、顾国华(2005)、王敏(2009)、王伟(2015)的研究结果一致，量级的差异由3.3.2节中MSF相关问题导致；(3)大兴安岭的隆升运动，速度为0.3 mm·a-1，与王伟(2013)一致；(4)华北地块呈现较强烈的隆升运动，速度达到1.5 mm·a-1，同样与王伟(2013)一致.华北地块的隆升速度大于喜马拉雅山脉的现象在以往的GPS观测结果(王伟，2013；顾国华，2005)中已有体现，其原因有待进一步研究.
在与已有成果进行对比分析中，图 12在青藏高原东南缘呈现的运动趋势与已有成果(Liang et al., 2013；王伟，2013；Cao et al，2015)存在差异，其原因可能为该地区的测站受到地震同震、震后，接收机、天线设备的更换，天线高量取等影响(王敏，2009).因此本文并没有对该区域进行深入分析，有待后续使用更高精度、更大密度的测站速度估值生成GNSS影像进行研究.
综合上述交叉验证结果分析，本文认为在中国大陆，GNSS影像能够正确地反映大范围的总体趋势，有助于地壳垂向形变的时空分布特征研究.
4 结论本文首先通过对美国西部地区以及中国大陆进行MSF验证实验与棋盘格检测，探求了两个测站网反映地壳垂向形变的空间分辨率及其可靠性.其次，使用两个区域内各测站未经过与经过顾及空间结构的中位数滤波的速度与不确定性估值，分别生成了GNSS影像；对比滤波前后的速度与不确定性生成的GNSS影像，说明了顾及空间结构的中位数滤波能够有效地剔除区域内粗差的影响，避免GNSS影像中出现零碎图斑，削弱稀疏测站网GNSS影像中规则圆弧形突变边缘问题的影响，是GNSS影像生成中一个重要的步骤；但同时也发现，MSF无法有效分辨测站速度估值中存在的速度峰值、突变等有效信息与粗差，导致滤波后生成的GNSS影像存在过度平滑的问题；针对该问题，本文结合滤波前后生成的GNSS影像，充分利用有效信息，初步提取了研究区域地壳形变的时空分布特征，如美国西部地区的海岸山脉与内华达山脉，中国大陆的天山地块、华北地区的隆升运动；美国西部地区维多利亚河及其河谷地区、洛杉矶地区，中国大陆华南地区、长江中下游地区、东北地区的沉降运动.这些结果与已有的大地测量学与地球动力学结果具有一致性，再次验证了GNSS影像能够正确反映地壳垂向形变的时空分布特征，GNSS成像方法具有合理性与正确性，有助于研究地壳垂向形变的时空分布特征，具有应用价值.
致谢 感谢美国内华达大学的William C. Hammond和Corné Kreemer教授在论文的实验期间多次与我们进行交流与讨论.感谢Renka R J编写构建球面Delaunay三角网的开源代码.本文中使用的GNSS坐标时间序列数据由UNAVCO与中国地震局发布；测站分布与GNSS影像成图使用的软件为GMT，在此表示感谢.
Blewitt G, Altamimi Z, Davis J, et al. 2010. Geodetic observations and global reference frame contributions to understanding sea-level rise and variability.//Church J A, Woodworth P L, Wilson W S eds. Understanding Sea-Level Rise and Variability. Chichester, UK: Wiley-Blackwell.
|
Blewitt G, Kreemer C, Hammond W C, et al. 2016. MIDAS robust trend estimator for accurate GPS station velocities without step detection. Journal of Geophysical Research:Solid Earth, 121(3): 2054-2068. DOI:10.1002/2015JB012552 |
Cao J L, Wang H, Wu Y Q, et al. 2015. Crustal deformation derived from GPS in the Chinese mainland from 1991-2004. Earthquake Research in China, 29(2): 225-236. |
Cormen T H, Leiserson C E, Rivest R L, et al. 2009. Introduction to Algorithms. 3rd ed. Cambridge, USA: The MIT Press.
|
Dong H. 2002. Research on Vertical Recent Crustal Movement of the Mainland of China (in Chinese). Xi'an: Xi'an Cartographic Publishing House..
|
Gu G H. 2005. Vertical crustal movement obtained from GPS observation in China's mainland. Earthquake (in Chinese), 25(3): 1-8. |
Guo L Q. 1990. Modern crustal vertical deformation in Northeast China block region and meaning of tectonic activity. Northeastern Seismological Research (in Chinese), 6(3): 15-20. |
Hamilton W, Myers W B. 1966. Cenozoic tectonics of the western United States. Reviews of Geophysics, 4(4): 509-549. |
Hammond W C, Blewitt G, Kreemer C. 2016. GPS imaging of vertical land motion in California and Nevada:implications for Sierra Nevada uplift. Journal of Geophysical Research:Solid Earth, 121(10): 7681-7703. DOI:10.1002/2016JB013458 |
Hu H M, Huang L R, Yang G H. 1992. Recent crustal vertical movement in the Chang-Jiang river delta and its adjacent area. Acta Geographica Sinica (in Chinese), 47(1): 22-30. |
Jiang G Y, Xu X W, Chen G H, et al. 2015. Geodetic imaging of potential Seismogenic Asperities on the Xianshuihe-Anninghe-Zemuhe fault system, Southwest China, with a new 3-D viscoelastic interseismic coupling model. Journal of Geophysical Research:Solid Earth, 120(3): 1855-1873. DOI:10.1002/2014JB011492 |
Jiang T C, Li T, Liu J N, et al. 2013. Deformation analysis of Wenchuan earthquake based on ScanSAR interferometry. Science of Surveying and Mapping (in Chinese), 38(6): 7-9. |
Lai X A, Huang L R, Xu J S. 2004. Present-Day Crustal Movement in China Constrained (in Chinese). Beijing: Seismological Press.
|
Liang S M, Gan W J, Shen C Z, et al. 2013. Three-dimensional velocity field of present-day crustal motion of the Tibetan plateau derived from GPS measurements. Journal of Geophysical Research:Solid Earth, 118(10): 5722-5732. DOI:10.1002/2013JB010503 |
Liu Z, Dong D N, Lundgren P. 2011. Constraints on time-dependent volcanic source models at long valley Caldera from 1996 to 2009 using InSAR and geodetic measurements. Geophysical Journal International, 187(3): 1283-1300. DOI:10.1111/j.1365-246X.2011.05214.x |
Reilinger R E, Citron G P, Brown L D. 1979. Recent vertical crustal movements from precise leveling data in Southwestern Montana, western Yellowstone national Park and the Snake River plain. Developments in Geotectonics, 13: 191-192. DOI:10.1016/B978-0-444-41783-1.50036-2 |
Renka R J. 1997. Algorithm 772:STRIPACK:Delaunay triangulation and Voronoi diagram on the surface of a sphere. ACM Transactions on Mathematical Software (TOMS), 23(3): 416-434. DOI:10.1145/275323.275329 |
Sen P K. 1968. Estimates of the regression coefficient based on Kendall's tau. Journal of the American Statistical Association, 63(324): 1379-1389. DOI:10.1080/01621459.1968.10480934 |
Shan X J, Zhang G H, Wang C S, et al. 2015. Joint inversion for the spatial fault slip distribution of the 2015 Nepal MW7.9 earthquake based on InSAR and GPS observation. Chinese Journal of Geophysics (in Chinese), 58(11): 4266-4276. DOI:10.6038/cjg20151131 |
Sun J B. 2009. Principles and Applications of Remote Sensing (in Chinese). Wuhan: Wuhan University Press.
|
Tan K, Zhao B, Zhang C, et al. 2016. Rupture models of the Nepal MW7.9 earthquake and MW7.3 aftershock constrained by GPS and InSAR coseismic deformation. Chinese Journal of Geophysics (in Chinese), 59(6): 2080-2093. DOI:10.6038/cjg20160614 |
Theil H. 1950. A rank-invariant method of linear and polynomial regression analysis.//Raj B, Koerts J eds. Henri Theil's Contributions to Economics and Econometrics. Dordrecht: Springer, 12: 85-91. https: //link.springer.com/chapter/10.1007/978-94-011-2546-8_20
|
Wang M. 2009. Analysis of GPS data with high precision & study on present-day crustal deformation in China[Ph. D. thesis] (in Chinese). Beijing: Institute of Geology, China Earthquake Administration.
|
Wang W. 2013. Study on present-day crustal movement with GPS data and modeling the active deformation in the Chinese mainland[Ph. D. thesis] (in Chinese). Beijing: Institute of Geophysics, China Earthquake Administration.
|
Wilcox R R. 2013. Introduction to Robust Estimation and Hypothesis Testing. 3rd ed. New York: Academic Press.
|
Xu C J, He P, Wen Y M, et al. 2015. Recent advances InSAR interferometry and its applications. Journal of Geomatics (in Chinese), 40(2): 1-9. |
Ying S F, Zhang Z S, Geng S C, et al. 1988. Basic characteristics of recent vertical crustal movement of China mainland. Earthquake Research in China (in Chinese), 4(4): 3-10. |
Zhang P Z, Deng Q D, Zhang G M, et al. 2003. Active tectonic blocks and strong earthquakes in the continent of China. Science in China Series D:Earth Sciences, 46(2): 13-24. |
Zhu J J, Li Z W, Hu J. 2017. Research progress and methods of InSAR for deformation monitoring. Acta Geodaetica et Cartographica Sinica (in Chinese), 46(10): 1717-1733. |
董鸿闻. 2002. 中国大陆现今地壳垂直运动研究. 西安: 西安地图出版社.
|
顾国华. 2005. GPS观测得到的中国大陆地壳垂直运动. 地震, 25(3): 1-8. |
郭良迁. 1990. 东北断块区的现代地壳垂直形变及其构造活动的意义. 东北地震研究, 6(3): 15-20. |
胡惠民, 黄立人, 杨国华. 1992. 长江三角洲及其邻近地区的现代地壳垂直运动. 地理学报, 47(1): 22-30. DOI:10.3321/j.issn:0375-5444.1992.01.004 |
蒋廷臣, 李陶, 刘经南, 等. 2013. ScanSAR干涉测量汶川地震形变场. 测绘科学, 38(6): 7-9. |
赖锡安, 黄立人, 徐菊生. 2004. 中国大陆现今地壳运动. 北京: 地震出版社.
|
单新建, 张国宏, 汪驰升, 等. 2015. 基于InSAR和GPS观测数据的尼泊尔地震发震断层特征参数联合反演研究. 地球物理学报, 58(11): 4266-4276. DOI:10.6038/cjg20151131 |
孙家抦. 2009. 遥感原理与应用. 武汉: 武汉大学出版社.
|
谭凯, 赵斌, 张彩红, 等. 2016. GPS和InSAR同震形变约束的尼泊尔MW7.9和MW7.3地震破裂滑动分布. 地球物理学报, 59(6): 2080-2093. DOI:10.6038/cjg20160614 |
王敏. 2009. GPS观测结果的精化分析与中国大陆现今地壳形变场研究[博士论文].北京: 中国地震局地质研究所. http://www.cnki.com.cn/Article/CJFDTotal-GJZT201003013.htm
|
王伟. 2013.中国大陆现今地壳运动的GPS分析与构造变形模拟[博士论文].北京: 中国地震局地球物理研究所. http://www.cnki.com.cn/Article/CJFDTotal-GJZT201402013.htm
|
许才军, 何平, 温扬茂, 等. 2015. InSAR技术及应用研究进展. 测绘地理信息, 40(2): 1-9. |
应绍奋, 张祖胜, 耿士昌, 等. 1988. 中国大陆垂直向现代地壳运动基本特征. 中国地震, 4(4): 3-10. |
张培震, 邓起东, 张国民, 等. 2003. 中国大陆的强震活动与活动地块. 中国科学(D辑), 33(S1): 12-20. |
朱建军, 李志伟, 胡俊. 2017. InSAR变形监测方法与研究进展. 测绘学报, 46(10): 1717-1733. DOI:10.11947/j.AGCS.2017.20170350 |