劉俊清，武成智，肖輝鋒

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

一、引 言

地殼運(yùn)動(dòng)過(guò)程可以看作是一個(gè)動(dòng)態(tài)系統(tǒng)，因而觀(guān)測(cè)點(diǎn)每個(gè)時(shí)刻的運(yùn)動(dòng)狀態(tài)是動(dòng)態(tài)(速度和加速度)變化的，而大地測(cè)量學(xué)通過(guò)多期重復(fù)測(cè)量通常獲得的是觀(guān)測(cè)點(diǎn)的靜態(tài)值，為獲得觀(guān)測(cè)點(diǎn)的動(dòng)態(tài)變化，需要在一個(gè)動(dòng)態(tài)系統(tǒng)中進(jìn)行數(shù)據(jù)處理。復(fù)測(cè)水準(zhǔn)測(cè)量作為研究地殼垂直形變手段已成為地震監(jiān)測(cè)預(yù)報(bào)的重要方法之一，網(wǎng)(路線(xiàn))型穩(wěn)定、觀(guān)測(cè)周期短，復(fù)測(cè)期數(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é)模型及觀(guān)測(cè)量的隨機(jī)模型，并進(jìn)行動(dòng)態(tài)系統(tǒng)的數(shù)據(jù)處理，在大地測(cè)量學(xué)中稱(chēng)這一過(guò)程為“動(dòng)態(tài)測(cè)量平差”。通過(guò)動(dòng)態(tài)測(cè)量平差來(lái)獲得相對(duì)真實(shí)的地殼運(yùn)動(dòng)狀態(tài)，其關(guān)鍵之處是選擇一個(gè)能夠準(zhǔn)確描述動(dòng)態(tài)系統(tǒng)的數(shù)學(xué)模型［4］?？柭鼮V波自從1960年被提出以后［5］，廣泛應(yīng)用于動(dòng)態(tài)系統(tǒng)中，已被證明是描述動(dòng)態(tài)系統(tǒng)最佳的算法。在地球物理研究中也應(yīng)用廣泛，如在監(jiān)測(cè)地殼水平運(yùn)動(dòng)的 VLBI、GPS 等測(cè)量手段中［6-7］，著名的 GPS數(shù)據(jù)后處理軟件GAMIT/GLOBK就是采用卡爾曼濾波做平面點(diǎn)的動(dòng)態(tài)平差，通過(guò)這種算法可以得到平面點(diǎn)相對(duì)真實(shí)的動(dòng)態(tài)運(yùn)動(dòng)狀態(tài)(速度和加速度)。

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

二、測(cè)量平差中的卡爾曼濾波方程

1.卡爾曼濾波原理

水準(zhǔn)測(cè)量觀(guān)測(cè)方程模型為

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

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

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

隨機(jī)模型為

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

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

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

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

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

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

圖1 長(zhǎng)白山天池火山水準(zhǔn)監(jiān)測(cè)路線(xiàn)

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

天池北側(cè)水準(zhǔn)路線(xiàn)截至2012年進(jìn)行了11期測(cè)量，西側(cè)水準(zhǔn)路線(xiàn)進(jìn)行了7期測(cè)量，獲得了長(zhǎng)白山天池火山區(qū)珍貴的垂直形變監(jiān)測(cè)資料。但是，火山區(qū)附近沒(méi)有找到任何等級(jí)的水準(zhǔn)測(cè)量標(biāo)志，即天池兩條水準(zhǔn)路線(xiàn)沒(méi)有控制點(diǎn)與起算點(diǎn)，無(wú)法應(yīng)用嚴(yán)格的經(jīng)典平差方法進(jìn)行數(shù)據(jù)處理。在以往的數(shù)據(jù)處理中利用了經(jīng)典測(cè)量平差進(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)測(cè)量顯然不嚴(yán)密，另外也不能獲得復(fù)測(cè)時(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ù)從“整米位”開(kāi)始)。其中，點(diǎn)號(hào)行表示各條水準(zhǔn)路線(xiàn)靠近火山口的3個(gè)水準(zhǔn)點(diǎn)的高程平差值;λ值行表示水準(zhǔn)點(diǎn)的變化速率。

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

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

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

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

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

西側(cè)線(xiàn)路2006年為第一期觀(guān)測(cè)，2007年之后的平差結(jié)果呈緩慢下降趨勢(shì)，第二期平差結(jié)果各點(diǎn)高程上升可能與路線(xiàn)水準(zhǔn)點(diǎn)沉降有關(guān)?？拷鹕娇诘?個(gè)水準(zhǔn)點(diǎn)在各期觀(guān)測(cè)中運(yùn)動(dòng)差異在誤差范圍內(nèi)，可以認(rèn)為速率一致，實(shí)際上西側(cè)線(xiàn)路其他水準(zhǔn)點(diǎn)的平差值變化速率與這3個(gè)點(diǎn)基本一致，推測(cè)可能與東北區(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é)束語(yǔ)

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

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

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

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

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

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

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

［9］于正林.多期復(fù)測(cè)水準(zhǔn)網(wǎng)的動(dòng)態(tài)平差［J］.武漢測(cè)繪學(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］∥中國(guó)科協(xié)2000年學(xué)術(shù)年會(huì)論文集.北京:［s.n.］，2000.

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

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

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

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

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

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