劉俊清，武成智，肖輝鋒

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

一、引 言

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

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

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

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

1.卡爾曼濾波原理

水準測量觀測方程模型為

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

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

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

隨機模型為

2.動態(tài)水準測量平差模型

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

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

三、長白山天池火山一等水準測量動態(tài)平差

1.長白山天池火山一等水準測量介紹

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

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

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

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

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

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

圖2 北側(cè)線路水準點時間序列

圖3 西側(cè)線路水準點時間序列

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

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

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

四、結(jié)束語

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

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

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

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

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

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

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

［9］于正林.多期復(fù)測水準網(wǎng)的動態(tài)平差［J］.武漢測繪學(xué)院學(xué)報，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)及其動力學(xué)意義［C］∥中國科協(xié)2000年學(xué)術(shù)年會論文集.北京:［s.n.］，2000.

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

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

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

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

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

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