卡爾曼濾波算法在長白山天池火山水準測量中的應用

劉俊清，武成智，肖輝鋒

(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.