GNSS垂向形变场信息提取的尝试

杨 博，占 伟，高艳龙，韩月萍

（中国地震局第一监测中心，天津 300180）

GNSS垂向形变场信息提取的尝试

杨 博，占 伟，高艳龙，韩月萍

（中国地震局第一监测中心，天津 300180）

一定条件下GNSS流动复测资料可为我们揭示地形变提供必要的依托，但由于垂向分量观测误差较大而没有被恰当地利用。利用华北地区1999—2007年的4期资料，从形变场的连续变化角度出发，提出了利用多核函数进行解析与滤波的方法，理论与结果表明：①可获得垂直运动场的数学解析；②可实现对垂向运动场不同层次的滤波与信息分离；③当流动测站具有一定的空间密度、一定数量的复测结果和时间长度则可获得较可靠的形变信息。最后利用华北地区GNSS流动观测资料给出了实际算例与初步分析结果。

GNSS；垂向运动；多核函数；滤波与解析

0 引言

在地球动力学研究领域里，以GNSS观测资料为依托的研究已是其重要的一个分支；具体到地壳形变研究而言，利用GNSS资料人们有望从空间与时间的角度观察地壳三维的运动及形变的详细过程。因此，多角度的研究也相应地展开［1－12］，但目前基本都是处在二维水平形变的研究层面上。这主要是因为流动GNSS垂向观测误差较大、影响的因素也较多，所以垂向观测分量在形变分析中一直没有得到应用。然而，随着复测资料的积累，垂向运动处理的信噪比也在不断的提高，这是一个不争的事实。为了更加较有效地利用垂向形变信息，本文试图从运动“观测值”入手，在形变场连续变化的前提下，通过一定的数理方法，实现垂直形变场的空间解析与滤波，以更快更好地利用垂向运动信息并服务于地壳形变研究、地表监测等领域。另一方面，大规模地开展精密水准测量在较短的时段内已不太可能，就华北而言，已近20年没有进行过全面的复测。所以，开发利用其他手段的监测已是非常必要了，GNSS或许是当今最切实际的一种选择。

1 垂直运动场的解析描述与滤波计算

1.1 多核函数与垂直运动场的解析表达

多核函数法是Hardy在1971年提出的一种数

式中sj（x，y，xj，yj）为核函数，（xj，yj）为核点的坐标，cj为待定系数。

一般说核函数是一个非常简单的曲面形态，通常为“钵”形或“钟”形曲面函数。尽管哪一种形态的核函数更优尚无非常明确的定论，但在此后的应用研究中，该法已被广泛地用于地壳形变研究［13－14］、DEM 内插［15］和地球重力场元内插［16］，并被证明是非常成功和有效的。

然而，多核函数在具体应用中的效果如何仍存在选择问题。对于不同的描述对象可能会有不同的选择，即使对相同的对象也并非是唯一选择，甚至很难区分它们之间的优劣。所以，在实际应用中主要取决于使用者对多核函数及描述对象的理解。一般取决于3个方面：①核函数的选择；②核函数中的参数选择；③核函数的核点选择。由于我们所描述的对象是较大尺度的地壳形变场，而地壳介质由于其各向非同性、活动构造带的存在及应力大小与方向的空间变化等，使得形变场具有广谱性。这就是说，我们所描述的对象是很复杂的。据此，这里所选用的核函数为“钵”形函数，具体为：值逼近方法，其基本思想是，任何数学表面和任何不规则的圆滑表面，总可以通过一系列有规则的数学表面的合成，以任意精度逼近。其数学表达式为：

式中：dj为球面上两点间的大地线长度（以km为单位），ST＝ （s1，s2，…，snx），CT＝ （c1，c2，…，cnx）。

在GNSS计算中通常是以ITRF（或WGS84）作为参考框架获取速度场，因此测站j的位置（λj，φj）、运动速率νu（λj，φj）及其相应的误差mu（λj，φj）都是已知的。在用多核函数进行数值逼近时应将所有测站位置作为核函数的核点位置，因此可列出n个方程（测站的个数），故可求解n个未知数（C）。这时即可获得垂向运动的解析式fu（λ，φ），同理也可获得误差解析式fm（λ，φ）。即

1.2 垂向运动的滤波计算

由于上述数值逼近所得到的是垂向运动观测值的解析结果，它不仅包含运动信息同时也包含误差干扰（运动观测值中所带来的），所以由此得到的运动结果实际上是各种成分的综合结果。因此，进行滤波与信息分离是必要的。多核函数不但可以进行数值逼近，而且也具有滤波的功能，故这里进行空间滤波同样利用多核函数法。进行滤波计算时，滤波效果或滤波层次则取决于核点的选用与分布。从性质上看，核点的间隔越小（不得小于测站间的平均间隔），由此所获得的运动信息所包含的频谱带宽越宽，可含盖从低频至高频的范围，反之带宽越窄，也越趋于低频。现假定垂向运动分量的滤波函数为：

式中：AT＝（a1，…，anx）为待定系数矩阵；ST＝ （s1，…，snx）为 核 函 数 矩 阵，si（λ，φ，λi，φi）＝ exp（－（di／780）2.8）（di为球面上两点间的大地线长度，以km 为单位），（λi，φi）为核点位置坐标。

在实际计算时，依据解析式（4）首先以相同的网格化形式（等步长）获取垂向分量和相应误差的数据，并以垂向分量作为“观测值”。然后，根据最小二乘法求解（5）式中的待定系数，即通过误差方程：

其中N＝DTPL，W＝DTPD，L为相应分量观测值的矩阵，P 为相应分量观测值的权矩阵（pi＝1／m2i）。当X求得后（即（5）式中的A），（5）式则变为经滤波后垂向运动计算的函数解析式，据此可获得任意位置垂向运动滤波值。

根据误差传播律，可进行滤波值的误差估计。由（6）式等可计算单位权方差估值：

式中nx为核点的个数，n为测站数量，且大于nx（否则无法给出精度评定）。垂向运动方差则为：

1.3 梯度的计算

考虑到信息的全面性，还可进一步给出梯度在空间上的变化结果。由于以上已经获得了运动场的解析式，所以梯度计算则变得非常方便了。根据梯度的定义，有：

据此就可很容易地给出梯度及其在空间的变化结果，并服务于垂直形变的全面分析。

2 算例与分析

华北上地壳最突出的现今动力学表现是新生代盆地的发育，而且盆地基底面的起伏与莫霍面的起伏呈上凹下凸型的镜像关系。这实际上表明深部活动对浅部构造的制约活动，其垂直活动特征为山区隆升，平原下降，并称之为常态活动。华北地形地貌的空间特征如图1所示：由北京－石家庄－郑州连线以西以北地区为山区或海拔较高的地区；其次是山东地区，其中较大一部分山区或丘陵地区，在其南部的安徽地区也有一些具有一定海拔高度的丘陵地貌；除此之外，基本上为海拔较低的平原地区（京津及以南的河北、河南和江苏等地区）。由于华北地区是我国地震活动较为强烈的地区之一，因此中国地壳运动网络在此布设了较密集的GNSS流动测站（测站分布为图2中的黑点所示），以用于监测该地区地壳形变及其动态变化。

图1 华北地形地貌空间特征略图

GNSS前期数据处理主要用较新版本GAMIT／GLOBK／QOCA软件完成，其步骤可概括为：①利用GAMIT软件获得测站坐标和卫星轨道的单日松弛解；②利用GLOBK软件将获得的每天单日松弛解和IGS数据中心SOPAC产出的全球IGS跟踪站有关的单日松弛解合并；③利用QOCA软件综合所有的单日松弛解计算各测站的三维运动速度。在利用QOCA软件求解测站速度时，采用ITRF05地球框架的速度解作为约束，联合解算了1999—2007年（1999、2001、2004和2007年4期观测资料）这一时段各测站包括垂向的三维运动速度。在由GNSS所获得的垂向运动速度的基础上，利用文中所述的方法进行进一步处理。其基本步骤是，首先依据（3）、（4）式对研究区相对运动速度及误差结果进行数据逼近，以获得解析式；然后以此基础以经、纬步长分别为35km间隔计算网格点上的垂向运动结果；接下来再利用（5）～（8）式进行空间滤波计算。滤波计算时选用的核点其经、纬间隔均为150km，结果见图2。

图2所显示的是华北地区1999—2007年垂向运动的趋势性结果，由这样的结果我们可获得第一印象是华北地区近十来年的运动形态在性质上与地形地貌形态是一致的，零值线的展布基本上体现了山区与平原的分界线；它展现了山区、丘陵等海拔较高的地区为相对隆升区（图2中的蓝色区域），海拔较低的平原地区为相对下沉区（图2中的红色区域）的特征。换句话说，华北地区当前的垂向运动方式仍处于继承性运动方式。但就空间的细部变化来看，不同的地区仍有较大的差别。华北地区西部的鄂尔多斯及周边地区整体虽表现为上升，平均上升接近3mm／a，但内部的差异变化并不大，变化范围0～5mm／a；此外，由于等值线分布宽缓，这表明其梯度变化是很小的，体现了内部运动具有很好的一致性；从地震的孕育角度看，内部差异运动较小的地域或块体难以储存较高的能量，因此也难以发生较大的地震，该区的地震历史也表明了这一点，同时也说明这一状态还将持续下去。鄂尔多斯东边缘的山西地域是华北地区隆升最突出的地域，除总体表现为上升外，最大上升量也位于其中，达10mm／a，其地点在太原以东地区。在山西境内隆升垂向运动变化范围为0～10mm／a，这表明山西地域的垂直形变是比较突出的，但并非均匀，主要集中在太原的周围地区。由于最小隆升位于太原以西地区，因此与其东侧最大隆升相比形成了鲜明的对照，说明该较小区域是能量相对积累最快的区域。除此之外，山西地域垂向运动的另一个特征是山西断陷带相对较小，其东部相对较大，具有区域上的继承性质。华北北端燕山带的整体隆升表现为“西小东大”，变化范围为2～6mm／a，形变相对均匀。山东的山区和丘陵地域基本上为隆升运动，但变化不均匀，大体以郯庐带为界，两侧上升，具有一定的对称形态，西侧最大上升为6mm／a，东侧为8mm／a，并存在着变形。其南的安徽地区也以隆升运动为主，最大达6mm／a，但等值线分布较为宽缓。以上所述的隆升运动，基本上体现的是近十年来的构造运动。然而，华北平原、江苏等地区运动发生了很大的变化，它们均处在下沉的状态之中。就性质而言，更多的属于地面下沉，尤其是石家庄以东、北京以南、济南以西等所围地区，明显表现为一个巨型的沉降“漏斗”，其中心沉降接近50mm／a。就其形成机理而言，基本上是过量开采地下水所致，并非属于构造运动的反映。河南和江苏地区等下沉量较小，一般在5mm／a以内，形变相对也较小。由于其“漏斗”的特征不明显，所以构造运动的比重相对较大些。

图2 华北地区1999—2007年垂向运动的趋势性结果（单位：mm／a）

如果不考虑沉降漏斗这一空间范围，图2所示结果基本上表征了华北近十年来具有构造活动含义的垂向形变场。因此，它告诉我们当前华北地区的垂向运动是继承性的。但在继承性运动中，或许存在着变形偏大的运动空间，如太原周围地区，故应给予关注。事实上，华北地区十几年来并未发生强烈地震；根据历史地震，一方面说明华北存在缺震现象；另一方面也可能说明缺震是有一定背景的。该区的继承性运动表明能量还处在积累之中，但对形变较大的部位，及无形变且历史上又是多震的部位不能掉以轻心。

3 结语

综上所述，在一定条件下GNSS的垂向运动信息是可以利用的，如果在GNSS的前期与后期数据处理进一步开展工作，可以获得服务于不同领域而又有效的形变信息。本文主要围绕后期数据方法而进行一些探讨和尝试。通过多核函数法实现了运动场的解析和滤波计算，同时也给出了运动场空间滤波结果的严密精度评定算式。实际结果表明，由此所获得的形变信息是有参考价值的。这在当前利用精密水准进行较大空间尺度垂直形变场研究难以维系的情况下，可以说是希望。然而，本方法更加恰当的使用还需结合实际情况、研究目标、及测站的分布等。核函数点的选择与滤波和信息分离有密切的关系，一般地说，核函数点间隔越大则信息越趋于低频，反之亦然。

通过对华北地区GNSS复测资料的处理与垂向形变信息的提取，使我们获得的最基本看法是，华北地区的构造运动在1999—2007年期间以继承性运动为主，即山区或海拔较高的地区为隆升，隆升量级绝大部分在4mm／a以内，极个别区域达8～10mm／a，平原地区为下降，可辨别的构造运动约在5mm／a以内，华北平原的下降主要属于过量开采地下流体造成的，其下沉量最大达50mm／a。

致谢：本项研究得到了杨国华研究员的具体指导，在此表示感谢。

［1］ 杨国华，张风霜，武艳强，等.利用 GPS连续观测资料进行强震危险性预测的探讨［J］.地震，2008，28（1）：33-39.

［2］ 黄立人，王敏.中国大陆构造块体的现今活动和变形［J］.地震地质，2003，25（1）：23-32.

［3］ 杨国华，李延兴，韩月萍，等.由 GPS观测结果推导中国大陆现今水平应变场［J］.地震学报，2002，24（4）：337-347.

［4］ 江在森，马宗晋，张希.GPS初步结果揭示的中国大陆水平应变场与构造变形［J］.地球物理学报，2003，46（3）：352-358.

［5］ 杨国华，江在森，王敏.印尼地震对我国川滇地区地壳水平活动的影响［J］.大地测量与地球动力学，2006，26（1）：9-14.

［6］ 杨博，周伟，陈阜超，等.GPS连续站水平位置时间序列共模白噪声识别与估计的欧拉－滤波法［J］.山东科技大学学报（自然科学版），2010，29（3）：26-31.

［7］ 杨国华，杨博，张风霜，等.汶川地震对华北地区水平形变场影响及有关含义的讨论［J］.地震，2009，29（1）：77-83.

［8］ 杨国华，韩月萍，杨博，等.华北地区最近几年水平形变场的含义与讨论［J］.山东科技大学学报（自然科学版），2009，28（6）：1-8.

［9］ 王敏，沈正康，牛之俊，等.现今中国大陆地壳运动与活动块体模型［J］.中国科学（D辑），2003，33（增）：21-32.

［10］ 王琪，丁国瑜，乔学军，等.天山现今地壳快速缩短与南北地块的相对运动［J］.科学通报，2000，45（14）：1543-1547.

［11］ 杨博，陈阜超，周伟，等.GPS连续站水平位置坐标时间序列插值的一种新方法［J］.华北地震科学，2009，28（2）：22-27.

［12］ 杨国华，江在森，武艳强，等.中国大陆整体无净旋转基准及其应用［J］.大地测量与地球动力学，2005，29（4）：6-10.

［13］ 黄立人，刘天奎.用速度面拟合法研究华北部分地区的现今地壳垂直运动［J］.地壳形变与地震，1989，9（3）：6-12.

［14］ 杨国华，黄立人.速度面拟合法中多面函数几个特性的初步数值研究［J］.地壳形变与地震，1990，10（4）：70-82.

［15］ 周光文，黄莜容.一种改进的多面函数法［J］.测绘通报，1996，（6）：6-8.

［16］ 陶本藻.GPS水准似大地水准面拟合和正常高计算［J］.测绘通报，1992，（4）：14-18.

［17］ 陈阜超，纪静，塔拉，等.京津水准复测与垂直形变特征［J］.华北地震科学，2011，（2）：31-34.

Some Attempt about Distilling Information of Vertical Deformation Field from GNSS Data

YANG Bo，ZHAN Wei，GAO Yan-long，HAN Yue-ping
（First Crust Monitoring and Application Center，CEA，Tianjin 300180，China）

In certain condition，GNSS mobile observations can provide essential support to reveal crustal deformation，but they were not used properly because of large observation error in vertical component.From the perspective of continuous variation of deformation field，a method by using multi-kernel function to filtering and analysis is proposed with 4observations（1999—2007）in the GPS network of North China，theory and results show as follows.Firstly，mathematical analysis of vertical motion field can be obtained.Secondly，different levels of filtering and information separation of vertical deformation field can be realized.Thirdly，if the space density of mobile stations，quantity and time length of repetition measurement can reach a certain level，reliable deformation information can be achieved.Finally，apractical example and preliminary analysis results with GNSS mobile observations in North China were given.

GNSS；vertical movement；multi-kernel function；filtering and analysis

P315.72

A

1003－1375（2012）01－0001－05

2011-07-15

地震行业科研专项重大项目（200908029）；地震行业科研专项（201008007）

杨博（1985－），男（汉族），天津人，助理工程师，主要从事GNSS数据处理与分析工作.E-mail：boyyang1234＠qq.com