卡尔曼滤波算法在长白山天池火山水准测量中的应用

刘俊清，武成智，肖辉锋

(1.吉林大学地球探测科学与技术学院，吉林长春130026;2.吉林省地震局，吉林长春130117;3.长白山火山监测站，吉林安图133600;4.北京望邦天鑫科技发展有限公司，北京100043)

一、引 言

地壳运动过程可以看作是一个动态系统，因而观测点每个时刻的运动状态是动态(速度和加速度)变化的，而大地测量学通过多期重复测量通常获得的是观测点的静态值，为获得观测点的动态变化，需要在一个动态系统中进行数据处理。复测水准测量作为研究地壳垂直形变手段已成为地震监测预报的重要方法之一，网(路线)型稳定、观测周期短，复测期数多，数据量大为其主要特点［1-2］，快速准确地进行数据处理从而获得地壳动态垂直变化速率是现代地震、火山预报的基本要求［3］。

选择一个系统合适的数学模型及观测量的随机模型，并进行动态系统的数据处理，在大地测量学中称这一过程为“动态测量平差”。通过动态测量平差来获得相对真实的地壳运动状态，其关键之处是选择一个能够准确描述动态系统的数学模型［4］。卡尔曼滤波自从1960年被提出以后［5］，广泛应用于动态系统中，已被证明是描述动态系统最佳的算法。在地球物理研究中也应用广泛，如在监测地壳水平运动的 VLBI、GPS 等测量手段中［6-7］，著名的 GPS数据后处理软件GAMIT/GLOBK就是采用卡尔曼滤波做平面点的动态平差，通过这种算法可以得到平面点相对真实的动态运动状态(速度和加速度)。

本文对长白山天池火山多期的水准测量结果进行真正意义上的动态平差，选择描述动态系统最佳的卡尔曼滤波方法，获得地壳垂直运动的速度及速度的变化率，对火山活动过程利用这种动态的测量结果进行深入研究，同时可以最大限度地提高资料的利用率。

二、测量平差中的卡尔曼滤波方程

1.卡尔曼滤波原理

水准测量观测方程模型为

求矩阵X的过程，即是统计学中的参数估计，在测量学中叫测量平差，通过L的观测值ln×1向量求定信号X的估值的方法称为滤波。在测量平差中，将含有误差的观测值求解参数的最佳估值的方法分为两类:一类是经典的最小二乘平差法，另一类就是滤波。由于滤波考虑了参数的随机性质，因此，当已知参数的先验特性时，它所得到的估值比经典的最小二乘估值具有更高的精度。在卡尔曼滤波中，当观测向量存在粗差时，平差结果将受到粗差影响。根据稳健M估计等价权原理［8］，通过分析增益矩阵，可以选取适当的权函数代替观测噪声协方差阵，以减小或消除粗差对估计结果的影响。水准测量动态平差过程通常分两步进行:首先根据抗差卡尔曼滤波模型建立水准平差模型，然后根据模型算法进行程序设计。

应用卡尔曼滤波方法对动态系统进行数据处理，就是通过观测向量L(t)估计随时间t不断变化的状态向量X(t)。严格地说，动态系统的状态是指全面确定动态系统运动状态的最少的一组参数。对于连续时间系统，在一般情况下，状态随时间t变化的规律，可用一个具有随机初始状态的向量微分方程来描述，称之为动态方程，而观测向量L(t)与状态向量X(t)之间的函数关系就是人们熟知的观测方程。因此动态系统的函数模型应包括状态方程和观测方程。

卡尔曼滤波基本状态方程和观测方程如下

随机模型为

2.动态水准测量平差模型

在动态水准监测系统中，若以点的高程及其速率为状态，根据上述基本方程推导复测水准测量状态方程和观测方程为［9］

这两个方程就是动态水准测量平差的抗差卡尔曼滤波算法递推方程，对其进行程序设计后，即可实现电算化。

三、长白山天池火山一等水准测量动态平差

1.长白山天池火山一等水准测量介绍

为有效地监测和研究长白山天池火山的地壳活动，深入分析火山区的地壳变形特征，以及为火山喷发危险性预测和评价提供基础的科学数据，从1999年开始，吉林省地震局陆续在天池火山锥体北侧和西侧建设了两条精密水准测量路线(如图1所示)，为监测火山锥体不同位置的垂直变形量，沿路线按一等精密水准测量规范设立固定点，如图北、西路线分别以字母N、W+序号表示。北侧水准路线位于天池和二道白河镇之间，北起景区山门以南4.5 km的黄松甸黄松蒲，南至天池瀑布，全长约25 km，相对高差901 m，共设立12个水准点，全部按一等水准测量水准标志建造。西侧水准路线2006年新建，位于天池至西山门之间，以天池西山门以西8 km起，东距天池2 km为止，全长约30 km，相对高差1084 m，共设立15个水准点，全部按一等水准测量水准标志建造，两条路线如图1所示。从2002年开始，每年8—9月间进行两条水准路线的测量，测量仪器使用Leica DNA03精密水准仪。

图1 长白山天池火山水准监测路线

2.数据处理

天池北侧水准路线截至2012年进行了11期测量，西侧水准路线进行了7期测量，获得了长白山天池火山区珍贵的垂直形变监测资料。但是，火山区附近没有找到任何等级的水准测量标志，即天池两条水准路线没有控制点与起算点，无法应用严格的经典平差方法进行数据处理。在以往的数据处理中利用了经典测量平差进行数据处理，但需要假定远离火山口的较稳定的山下起始点作为固定点。实际上在地球动力学范畴内，这个假定的起始点也是运动的，经典的静态平差处理这种超一等的水准测量显然不严密，另外也不能获得复测时间段内地表的动态运动信息。使用卡尔曼滤波算法进行动态平差，可以获得真正意义上的地壳动态运动信息。表1列出了2006—2012年利用卡尔曼滤波算法的平差结果(为制表方便，数据从“整米位”开始)。其中，点号行表示各条水准路线靠近火山口的3个水准点的高程平差值;λ值行表示水准点的变化速率。

将两条水准路线近火山口的3个水准点高程的平差值分别画出时间序列图(如图2、图3所示)。由平差值曲线图可知，北侧线路在多期的测量中一直呈上升趋势，上升速率最大时间段为2002—2003年，且上升速率不一致，N_11、N_10、N_9水准点沿火山口向下速率依次减小，2003年之后上升速率在逐年减小，北侧线路N_0点至N_8水准点没有明显的运动速率。这种变形特征完全符合火山活动时的点源火山膨胀模型［10］，这种火山模型是把地表以下一定深度的岩浆囊作为点源，由于岩浆囊的膨胀，会使地表隆起，随着地表位置的变化，垂直变形量具有一定规律。根据长白山电磁测深资料［11-12］，地下10 km处有低速体的存在，根据点源膨胀模型计算，地表变形最明显范围是以天池火山口为中心，半径为8 km的范围内，N_11、N_10、N_9恰好位于这一区域。

表1 长白山天池火山水准高差卡尔曼滤波结果 mm

图2 北侧线路水准点时间序列

图3 西侧线路水准点时间序列

2002年天池火山经历一次明显的火山扰动事件［13］，由火山地震震级与时间的关系图及地震频度图(如图4所示)显示，2002年年底火山地震频度很高，进入2003年之后火山地震频次逐年减少，地震频度逐渐衰减的特征表明天池火山的确经历过一次短时间的扰动过程。根据2002—2007年天池火山区的流动GPS监测结果，2002—2003年火山区地表存在显著的面膨胀［14-16］。显然，通过卡尔曼滤波算法获得的天池火山北侧水准路线的动态平差结果，较好地反映了这次火山扰动事件，由于火山锥体的垂直变形特征符合点源模型，也说明的确是地表下小尺度岩浆囊的小规模膨胀。

西侧线路2006年为第一期观测，2007年之后的平差结果呈缓慢下降趋势，第二期平差结果各点高程上升可能与路线水准点沉降有关。靠近火山口的3个水准点在各期观测中运动差异在误差范围内，可以认为速率一致，实际上西侧线路其他水准点的平差值变化速率与这3个点基本一致，推测可能与东北区域地壳构造垂直形变速率一致［17］。火山口两侧垂直形变不一致，初步认为2002年的岩浆扰动事件发生在锥体北侧地下，而整个火山区地壳垂直运动背景与东北地壳运动速率一致。

图4 天池火山地震M-T图及频度图

四、结束语

本文在无起算点的长白山天池火山区水准路线平差中应用了卡尔曼滤波算法作水准测量的动态平差，把水准线路上的所有水准点看作垂直方向自由运动的点，没有假定固定点，可以计算出所有点的运动状态。利用长白山天池火山水准测量数据的平差结果修正了原始观测值，计算出更为合理的点运动速率，其时间序列曲线的变化很好地反映了2002年的天池火山活动，水准路线上设立的固定点垂直变化量完全符合火山区点源理论模型的变形特征，表明卡尔曼滤波算法可用于火山区的复测水准测量动态平差。

［1］薄志鹏，刘国辉，汪晓庆.对全国复测水准网平差的几点探讨［J］.武汉测绘科技大学学报，1991，16(3):31-37.

［2］刘晓云，张世娟，程传录.精密水准测量数据处理自动化系统的研究与实现［J］.测绘通报，2013(10):67-69，86.

［3］刘俊清，孙继财，武成智，等.长白山天池火山溃湖洪水最大流量的初步估算及影响分析［J］.地震地质，2013，40(1):85-91.

［4］郝明.基于精密水准数据的青藏高原东缘现今地壳垂直运动与典型地震同震及震后垂直形变研究［D］.北京:中国地震局地质研究所，2012.

［5］KALMAN R E.A New Approach to Linear Filtering and Prediction Problems［J］.Journal of Basic Engineering，1960，82(1):35-45.

［6］房建成，万德钧，吴秋平，等.GPS动态滤波的新方法［J］.中国惯性技术学报，1997，5(2):4-10.

［7］刘站科，付文举，黄观文，等.GPS/GLONASS组合精密单点定位［J］.测绘通报，2014(7):6-10.

［8］杨元喜.自适应抗差滤波理论及应用的主要进展［C］∥中国测绘学会第九次全国会员代表大会暨学会成立 50 周年纪念大会论文集.北京:［s.n.］，2009.

［9］于正林.多期复测水准网的动态平差［J］.武汉测绘学院学报，1981，6(2):66-76.

［10］MOGI K，茂木清夫.Flow and Fracture of Rocks under General Triaxial Compression［J］.Applied Mathematics and Mechanics:English Edition，1981，111(6):635-651.

［11］张先康，张成科，嘉世旭，等.天山地壳结构及其动力学意义［C］∥中国科协2000年学术年会论文集.北京:［s.n.］，2000.

［12］汤吉，邓前辉，赵国泽，等.长白山天池火山区电性结构和岩浆系统［J］.地震地质，2001，32(2):191-200.

［13］吴建平，明跃红，张恒荣，等.长白山天池火山区的震群活动研究［J］.地球物理学报，2007，S0(4):1089-1096.

［14］刘俊清.GPS测量技术在长白山火山形变监测中的应用［D］.长春:吉林大学，2010.

［15］胡亚轩，王庆良，崔笃信，等.长白山天池火山区形变监测及火山活动状态分析［J］.大地测量与地球动力学，2007，32(5):22-5.

［16］崔笃信，王庆良，李克，等.长白山天池火山近期形变场演化过程分析［J］.地球物理学报，2007，S0(6):1731-1739.

［17］胡亚轩，王庆良，王雄.利用垂直形变资料分析龙岗火山的活动性［J］.地震研究，2009，32(3):289-294.