卡爾曼濾波算法在長白山天池火山水準(zhǔn)測量中的應(yīng)用

劉俊清，武成智，肖輝鋒

(1.吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，吉林長春130026;2.吉林省地震局，吉林長春130117;3.長白山火山監(jiān)測站，吉林安圖133600;4.北京望邦天鑫科技發(fā)展有限公司，北京100043)

一、引 言

地殼運(yùn)動(dòng)過程可以看作是一個(gè)動(dòng)態(tài)系統(tǒng)，因而觀測點(diǎn)每個(gè)時(shí)刻的運(yùn)動(dòng)狀態(tài)是動(dòng)態(tài)(速度和加速度)變化的，而大地測量學(xué)通過多期重復(fù)測量通常獲得的是觀測點(diǎn)的靜態(tài)值，為獲得觀測點(diǎn)的動(dòng)態(tài)變化，需要在一個(gè)動(dòng)態(tài)系統(tǒng)中進(jìn)行數(shù)據(jù)處理。復(fù)測水準(zhǔn)測量作為研究地殼垂直形變手段已成為地震監(jiān)測預(yù)報(bào)的重要方法之一，網(wǎng)(路線)型穩(wěn)定、觀測周期短，復(fù)測期數(shù)多，數(shù)據(jù)量大為其主要特點(diǎn)［1-2］，快速準(zhǔn)確地進(jìn)行數(shù)據(jù)處理從而獲得地殼動(dòng)態(tài)垂直變化速率是現(xiàn)代地震、火山預(yù)報(bào)的基本要求［3］。

選擇一個(gè)系統(tǒng)合適的數(shù)學(xué)模型及觀測量的隨機(jī)模型，并進(jìn)行動(dòng)態(tài)系統(tǒng)的數(shù)據(jù)處理，在大地測量學(xué)中稱這一過程為“動(dòng)態(tài)測量平差”。通過動(dòng)態(tài)測量平差來獲得相對真實(shí)的地殼運(yùn)動(dòng)狀態(tài)，其關(guān)鍵之處是選擇一個(gè)能夠準(zhǔn)確描述動(dòng)態(tài)系統(tǒng)的數(shù)學(xué)模型［4］。卡爾曼濾波自從1960年被提出以后［5］，廣泛應(yīng)用于動(dòng)態(tài)系統(tǒng)中，已被證明是描述動(dòng)態(tài)系統(tǒng)最佳的算法。在地球物理研究中也應(yīng)用廣泛，如在監(jiān)測地殼水平運(yùn)動(dòng)的 VLBI、GPS 等測量手段中［6-7］，著名的 GPS數(shù)據(jù)后處理軟件GAMIT/GLOBK就是采用卡爾曼濾波做平面點(diǎn)的動(dòng)態(tài)平差，通過這種算法可以得到平面點(diǎn)相對真實(shí)的動(dòng)態(tài)運(yùn)動(dòng)狀態(tài)(速度和加速度)。

本文對長白山天池火山多期的水準(zhǔn)測量結(jié)果進(jìn)行真正意義上的動(dòng)態(tài)平差，選擇描述動(dòng)態(tài)系統(tǒng)最佳的卡爾曼濾波方法，獲得地殼垂直運(yùn)動(dòng)的速度及速度的變化率，對火山活動(dòng)過程利用這種動(dòng)態(tài)的測量結(jié)果進(jìn)行深入研究，同時(shí)可以最大限度地提高資料的利用率。

二、測量平差中的卡爾曼濾波方程

1.卡爾曼濾波原理

水準(zhǔn)測量觀測方程模型為

求矩陣X的過程，即是統(tǒng)計(jì)學(xué)中的參數(shù)估計(jì)，在測量學(xué)中叫測量平差，通過L的觀測值ln×1向量求定信號X的估值的方法稱為濾波。在測量平差中，將含有誤差的觀測值求解參數(shù)的最佳估值的方法分為兩類:一類是經(jīng)典的最小二乘平差法，另一類就是濾波。由于濾波考慮了參數(shù)的隨機(jī)性質(zhì)，因此，當(dāng)已知參數(shù)的先驗(yàn)特性時(shí)，它所得到的估值比經(jīng)典的最小二乘估值具有更高的精度。在卡爾曼濾波中，當(dāng)觀測向量存在粗差時(shí)，平差結(jié)果將受到粗差影響。根據(jù)穩(wěn)健M估計(jì)等價(jià)權(quán)原理［8］，通過分析增益矩陣，可以選取適當(dāng)?shù)臋?quán)函數(shù)代替觀測噪聲協(xié)方差陣，以減小或消除粗差對估計(jì)結(jié)果的影響。水準(zhǔn)測量動(dòng)態(tài)平差過程通常分兩步進(jìn)行:首先根據(jù)抗差卡爾曼濾波模型建立水準(zhǔn)平差模型，然后根據(jù)模型算法進(jìn)行程序設(shè)計(jì)。

應(yīng)用卡爾曼濾波方法對動(dòng)態(tài)系統(tǒng)進(jìn)行數(shù)據(jù)處理，就是通過觀測向量L(t)估計(jì)隨時(shí)間t不斷變化的狀態(tài)向量X(t)。嚴(yán)格地說，動(dòng)態(tài)系統(tǒng)的狀態(tài)是指全面確定動(dòng)態(tài)系統(tǒng)運(yùn)動(dòng)狀態(tài)的最少的一組參數(shù)。對于連續(xù)時(shí)間系統(tǒng)，在一般情況下，狀態(tài)隨時(shí)間t變化的規(guī)律，可用一個(gè)具有隨機(jī)初始狀態(tài)的向量微分方程來描述，稱之為動(dòng)態(tài)方程，而觀測向量L(t)與狀態(tài)向量X(t)之間的函數(shù)關(guān)系就是人們熟知的觀測方程。因此動(dòng)態(tài)系統(tǒng)的函數(shù)模型應(yīng)包括狀態(tài)方程和觀測方程。

卡爾曼濾波基本狀態(tài)方程和觀測方程如下

隨機(jī)模型為

2.動(dòng)態(tài)水準(zhǔn)測量平差模型

在動(dòng)態(tài)水準(zhǔn)監(jiān)測系統(tǒng)中，若以點(diǎn)的高程及其速率為狀態(tài)，根據(jù)上述基本方程推導(dǎo)復(fù)測水準(zhǔn)測量狀態(tài)方程和觀測方程為［9］

這兩個(gè)方程就是動(dòng)態(tài)水準(zhǔn)測量平差的抗差卡爾曼濾波算法遞推方程，對其進(jìn)行程序設(shè)計(jì)后，即可實(shí)現(xiàn)電算化。

三、長白山天池火山一等水準(zhǔn)測量動(dòng)態(tài)平差

1.長白山天池火山一等水準(zhǔn)測量介紹

為有效地監(jiān)測和研究長白山天池火山的地殼活動(dòng)，深入分析火山區(qū)的地殼變形特征，以及為火山噴發(fā)危險(xiǎn)性預(yù)測和評價(jià)提供基礎(chǔ)的科學(xué)數(shù)據(jù)，從1999年開始，吉林省地震局陸續(xù)在天池火山錐體北側(cè)和西側(cè)建設(shè)了兩條精密水準(zhǔn)測量路線(如圖1所示)，為監(jiān)測火山錐體不同位置的垂直變形量，沿路線按一等精密水準(zhǔn)測量規(guī)范設(shè)立固定點(diǎn)，如圖北、西路線分別以字母N、W+序號表示。北側(cè)水準(zhǔn)路線位于天池和二道白河鎮(zhèn)之間，北起景區(qū)山門以南4.5 km的黃松甸黃松蒲，南至天池瀑布，全長約25 km，相對高差901 m，共設(shè)立12個(gè)水準(zhǔn)點(diǎn)，全部按一等水準(zhǔn)測量水準(zhǔn)標(biāo)志建造。西側(cè)水準(zhǔn)路線2006年新建，位于天池至西山門之間，以天池西山門以西8 km起，東距天池2 km為止，全長約30 km，相對高差1084 m，共設(shè)立15個(gè)水準(zhǔn)點(diǎn)，全部按一等水準(zhǔn)測量水準(zhǔn)標(biāo)志建造，兩條路線如圖1所示。從2002年開始，每年8—9月間進(jìn)行兩條水準(zhǔn)路線的測量，測量儀器使用Leica DNA03精密水準(zhǔn)儀。

圖1 長白山天池火山水準(zhǔn)監(jiān)測路線

2.數(shù)據(jù)處理

天池北側(cè)水準(zhǔn)路線截至2012年進(jìn)行了11期測量，西側(cè)水準(zhǔn)路線進(jìn)行了7期測量，獲得了長白山天池火山區(qū)珍貴的垂直形變監(jiān)測資料。但是，火山區(qū)附近沒有找到任何等級的水準(zhǔn)測量標(biāo)志，即天池兩條水準(zhǔn)路線沒有控制點(diǎn)與起算點(diǎn)，無法應(yīng)用嚴(yán)格的經(jīng)典平差方法進(jìn)行數(shù)據(jù)處理。在以往的數(shù)據(jù)處理中利用了經(jīng)典測量平差進(jìn)行數(shù)據(jù)處理，但需要假定遠(yuǎn)離火山口的較穩(wěn)定的山下起始點(diǎn)作為固定點(diǎn)。實(shí)際上在地球動(dòng)力學(xué)范疇內(nèi)，這個(gè)假定的起始點(diǎn)也是運(yùn)動(dòng)的，經(jīng)典的靜態(tài)平差處理這種超一等的水準(zhǔn)測量顯然不嚴(yán)密，另外也不能獲得復(fù)測時(shí)間段內(nèi)地表的動(dòng)態(tài)運(yùn)動(dòng)信息。使用卡爾曼濾波算法進(jìn)行動(dòng)態(tài)平差，可以獲得真正意義上的地殼動(dòng)態(tài)運(yùn)動(dòng)信息。表1列出了2006—2012年利用卡爾曼濾波算法的平差結(jié)果(為制表方便，數(shù)據(jù)從“整米位”開始)。其中，點(diǎn)號行表示各條水準(zhǔn)路線靠近火山口的3個(gè)水準(zhǔn)點(diǎn)的高程平差值;λ值行表示水準(zhǔn)點(diǎn)的變化速率。

將兩條水準(zhǔn)路線近火山口的3個(gè)水準(zhǔn)點(diǎn)高程的平差值分別畫出時(shí)間序列圖(如圖2、圖3所示)。由平差值曲線圖可知，北側(cè)線路在多期的測量中一直呈上升趨勢，上升速率最大時(shí)間段為2002—2003年，且上升速率不一致，N_11、N_10、N_9水準(zhǔn)點(diǎn)沿火山口向下速率依次減小，2003年之后上升速率在逐年減小，北側(cè)線路N_0點(diǎn)至N_8水準(zhǔn)點(diǎn)沒有明顯的運(yùn)動(dòng)速率。這種變形特征完全符合火山活動(dòng)時(shí)的點(diǎn)源火山膨脹模型［10］，這種火山模型是把地表以下一定深度的巖漿囊作為點(diǎn)源，由于巖漿囊的膨脹，會(huì)使地表隆起，隨著地表位置的變化，垂直變形量具有一定規(guī)律。根據(jù)長白山電磁測深資料［11-12］，地下10 km處有低速體的存在，根據(jù)點(diǎn)源膨脹模型計(jì)算，地表變形最明顯范圍是以天池火山口為中心，半徑為8 km的范圍內(nèi)，N_11、N_10、N_9恰好位于這一區(qū)域。

表1 長白山天池火山水準(zhǔn)高差卡爾曼濾波結(jié)果 mm

圖2 北側(cè)線路水準(zhǔn)點(diǎn)時(shí)間序列

圖3 西側(cè)線路水準(zhǔn)點(diǎn)時(shí)間序列

2002年天池火山經(jīng)歷一次明顯的火山擾動(dòng)事件［13］，由火山地震震級與時(shí)間的關(guān)系圖及地震頻度圖(如圖4所示)顯示，2002年年底火山地震頻度很高，進(jìn)入2003年之后火山地震頻次逐年減少，地震頻度逐漸衰減的特征表明天池火山的確經(jīng)歷過一次短時(shí)間的擾動(dòng)過程。根據(jù)2002—2007年天池火山區(qū)的流動(dòng)GPS監(jiān)測結(jié)果，2002—2003年火山區(qū)地表存在顯著的面膨脹［14-16］。顯然，通過卡爾曼濾波算法獲得的天池火山北側(cè)水準(zhǔn)路線的動(dòng)態(tài)平差結(jié)果，較好地反映了這次火山擾動(dòng)事件，由于火山錐體的垂直變形特征符合點(diǎn)源模型，也說明的確是地表下小尺度巖漿囊的小規(guī)模膨脹。

西側(cè)線路2006年為第一期觀測，2007年之后的平差結(jié)果呈緩慢下降趨勢，第二期平差結(jié)果各點(diǎn)高程上升可能與路線水準(zhǔn)點(diǎn)沉降有關(guān)。靠近火山口的3個(gè)水準(zhǔn)點(diǎn)在各期觀測中運(yùn)動(dòng)差異在誤差范圍內(nèi)，可以認(rèn)為速率一致，實(shí)際上西側(cè)線路其他水準(zhǔn)點(diǎn)的平差值變化速率與這3個(gè)點(diǎn)基本一致，推測可能與東北區(qū)域地殼構(gòu)造垂直形變速率一致［17］。火山口兩側(cè)垂直形變不一致，初步認(rèn)為2002年的巖漿擾動(dòng)事件發(fā)生在錐體北側(cè)地下，而整個(gè)火山區(qū)地殼垂直運(yùn)動(dòng)背景與東北地殼運(yùn)動(dòng)速率一致。

圖4 天池火山地震M-T圖及頻度圖

四、結(jié)束語

本文在無起算點(diǎn)的長白山天池火山區(qū)水準(zhǔn)路線平差中應(yīng)用了卡爾曼濾波算法作水準(zhǔn)測量的動(dòng)態(tài)平差，把水準(zhǔn)線路上的所有水準(zhǔn)點(diǎn)看作垂直方向自由運(yùn)動(dòng)的點(diǎn)，沒有假定固定點(diǎn)，可以計(jì)算出所有點(diǎn)的運(yùn)動(dòng)狀態(tài)。利用長白山天池火山水準(zhǔn)測量數(shù)據(jù)的平差結(jié)果修正了原始觀測值，計(jì)算出更為合理的點(diǎn)運(yùn)動(dòng)速率，其時(shí)間序列曲線的變化很好地反映了2002年的天池火山活動(dòng)，水準(zhǔn)路線上設(shè)立的固定點(diǎn)垂直變化量完全符合火山區(qū)點(diǎn)源理論模型的變形特征，表明卡爾曼濾波算法可用于火山區(qū)的復(fù)測水準(zhǔn)測量動(dòng)態(tài)平差。

［1］薄志鵬，劉國輝，汪曉慶.對全國復(fù)測水準(zhǔn)網(wǎng)平差的幾點(diǎn)探討［J］.武漢測繪科技大學(xué)學(xué)報(bào)，1991，16(3):31-37.

［2］劉曉云，張世娟，程傳錄.精密水準(zhǔn)測量數(shù)據(jù)處理自動(dòng)化系統(tǒng)的研究與實(shí)現(xiàn)［J］.測繪通報(bào)，2013(10):67-69，86.

［3］劉俊清，孫繼財(cái)，武成智，等.長白山天池火山潰湖洪水最大流量的初步估算及影響分析［J］.地震地質(zhì)，2013，40(1):85-91.

［4］郝明.基于精密水準(zhǔn)數(shù)據(jù)的青藏高原東緣現(xiàn)今地殼垂直運(yùn)動(dòng)與典型地震同震及震后垂直形變研究［D］.北京:中國地震局地質(zhì)研究所，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動(dòng)態(tài)濾波的新方法［J］.中國慣性技術(shù)學(xué)報(bào)，1997，5(2):4-10.

［7］劉站科，付文舉，黃觀文，等.GPS/GLONASS組合精密單點(diǎn)定位［J］.測繪通報(bào)，2014(7):6-10.

［8］楊元喜.自適應(yīng)抗差濾波理論及應(yīng)用的主要進(jìn)展［C］∥中國測繪學(xué)會(huì)第九次全國會(huì)員代表大會(huì)暨學(xué)會(huì)成立 50 周年紀(jì)念大會(huì)論文集.北京:［s.n.］，2009.

［9］于正林.多期復(fù)測水準(zhǔn)網(wǎng)的動(dòng)態(tài)平差［J］.武漢測繪學(xué)院學(xué)報(bào)，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］張先康，張成科，嘉世旭，等.天山地殼結(jié)構(gòu)及其動(dòng)力學(xué)意義［C］∥中國科協(xié)2000年學(xué)術(shù)年會(huì)論文集.北京:［s.n.］，2000.

［12］湯吉，鄧前輝，趙國澤，等.長白山天池火山區(qū)電性結(jié)構(gòu)和巖漿系統(tǒng)［J］.地震地質(zhì)，2001，32(2):191-200.

［13］吳建平，明躍紅，張恒榮，等.長白山天池火山區(qū)的震群活動(dòng)研究［J］.地球物理學(xué)報(bào)，2007，S0(4):1089-1096.

［14］劉俊清.GPS測量技術(shù)在長白山火山形變監(jiān)測中的應(yīng)用［D］.長春:吉林大學(xué)，2010.

［15］胡亞軒，王慶良，崔篤信，等.長白山天池火山區(qū)形變監(jiān)測及火山活動(dòng)狀態(tài)分析［J］.大地測量與地球動(dòng)力學(xué)，2007，32(5):22-5.

［16］崔篤信，王慶良，李克，等.長白山天池火山近期形變場演化過程分析［J］.地球物理學(xué)報(bào)，2007，S0(6):1731-1739.

［17］胡亞軒，王慶良，王雄.利用垂直形變資料分析龍崗火山的活動(dòng)性［J］.地震研究，2009，32(3):289-294.