本文内容來源于《測繪學報》2024年第1期（審圖号GS京(2024)0107号）
基于變分模态分解的GNSS高程時間序列時變信号提取武曙光1
, 邊少鋒2, 李厚樸1
, 李昭3, 歐陽華11. 海軍工程大學電氣工程學院, 湖北 武漢 430034;
2. 中國地質大學(武漢)地質探測與評估教育部重點實驗室, 湖北 武漢 430074;
3. 武漢大學衛星導航定位技術研究中心, 湖北 武漢 430079基金項目：國家自然科學基金(42122025;42174030;41771487);湖北珞珈實驗室專項基金(220100020);資源與環境資訊系統國家重點實驗室開放基金摘要：針對GNSS坐标時間序列中的時變信号難以由現有最小二乘、最大似然估計(MLE)等參數化方法準确提取的問題, 本文采用變分模态分解(VMD)方法将中國大陸構造環境監測網絡(CMONOC)測站的高程時間序列分解為一系列本征模态函數(IMF), 進而重構出測站位置時間序列中含有的時變信号。結果表明, 相對于MLE方法, VMD方法在97.9%的測站上均方根誤差(RMSE)改進率為正值, 是以該方法有助于絕大多數測站精确提取出時變信号, 減弱高程時間序列中的非線性形變。另外, 從相關系數和信噪比的角度來看, VMD方法得到的重構序列與原始序列之間的相關系數更高, 信噪比也更大, 表明降噪效果較好。通過特定測站的分析表明, VMD方法能有效探測出GNSS高程時間序列預進行中包含遺漏的階躍信号的測站, 表現為較大的RMSE改進率, 這在大批量測站的階躍信号探測中具有一定的實用價值。VMD方法相對于小波分解(WD)經驗模态分解(EMD)具有更好的自适應性, 但IMF分量個數仍然需要針對具體測站進行逐一确定, 當分解個數和重構分量選取恰當時, VMD方法在GNSS高程時間序列中的應用效果可進一步提高。關鍵詞：GNSS高程時間序列 變分模态分解 CMONOC測站 RMSE改進率
引文格式：武曙光, 邊少鋒, 李厚樸, 等. 基于變分模态分解的GNSS高程時間序列時變信号提取[J]. 測繪學報，2024，53(1)：79-90. DOI: 10.11947/j.AGCS.2024.20220673
WU Shuguang, BIAN Shaofeng, LI Houpu, et al. Extraction of time-varying signals from GNSS height time series by variational mode decomposition[J]. Acta Geodaetica et Cartographica Sinica, 2024, 53(1): 79-90. DOI: 10.11947/j.AGCS.2024.20220673 閱讀全文：http://xb.chinasmp.com/article/2024/1001-1595/20240107.htm
引 言
20多年來，中國大陸構造環境監測網絡(Crustal Movement Observation Network of China, CMONOC, 簡稱“陸态網”)測站累積了大量的GNSS觀測檔案，通過GNSS精密資料處理，得到測站長時期的坐标時間序列，為大地測量學、地球動力學等學科的研究提供了有效的資料保障[1-3]。GNSS坐标時間序列一般包含了由地殼運動造成的長期線性變化，同時也包含各類空間尺度的地球實體效應(如大氣、水文、潮汐等)導緻的非線性變化，主要表現為測站位置的季節性震蕩。此外，測站位置還包含了由于地震、火山活動等構造運動，以及儀器更換等因素引起的階躍和震後形變。另外，GNSS坐标時間序列還可能含有因資料處理模型模組化不完善或錯誤模組化導緻的虛假(非線性)變化。近20年來，國内外學者開展了衆多非線性形變來源的研究工作[4-6]。研究表明，通過建立地表負載模型(surface loading model, SLM)，并計算GNSS測站位置的彈性形變位移序列，可以有效研究測站位置處的地球實體影響機制。目前，地表負載模型最理想情況下能夠解釋約50%的測站垂向周年振幅，水準方向則不到20%[7-8]。除了環境負載，GNSS測站還受到熱膨脹效應的影響，表現出一種周期性溫度變化引起的天線觀測墩及其所在基岩的熱彈性形變，也被認為是GNSS坐标時間序列非線性信号的來源之一[9-10]。除此之外，區域性的地球實體因素也會導緻GNSS測站産生非線性形變。例如，區域内地下水、石油等的開采、地表強降水[11]、高原凍土周期性當機與消融[12]等。研究表明，大陸的GNSS坐标時間序列在水準方向主要以線性運動為主[13-14]，而高程方向則表現為明顯的季節性振蕩，且不同測站的振蕩特性并不相同。對于坐标高程分量中的季節性信号，常将其視為年周期信号及其整數階諧波分量的組合[15-16]，并用最小二乘、最大似然估計(maximum likelihood estimation, MLE)等方法估計出常數振幅和相位。這種參數化的諧波模型對于平穩序列是合适的，而對于非平穩的測站，年際間的測站環境變化及其響應存在差異時，不變的振幅和相位就不再适合描述測站的運動特征，這時就需要将季節性信号的振幅和相位視為随時間變化的。目前，GNSS坐标時間序列分析中常用的時變信号提取方法有小波分解(wavelet decomposition, WD)和經驗模态分解(empirical mode decomposition, EMD)[17-19]。研究表明，WD的效果受小波基的選擇及分解層數的影響較大，當依靠經驗選擇不當時，容易使結果出現偏差[20]。相對而言，EMD按照資料自身的時間尺度特征将信号分解為不同的本征模态函數(intrinsic mode functions, IMF)，是一種自适應的方法。自文獻[21]提出至今，已在大地測量各類資料進行中有了廣泛的應用[22-24]。然而研究表明，對于較為複雜的信号，EMD分解過程中往往會出現模态混疊現象[25]，造成不同的模态在同一個IMF分量中共存或相同的模态分成了多個頻率相近的模态，因而無法根據時間特征尺度有效區分不同模态。為解決模态混疊問題，文獻[26]提出了變分模态分解(variational mode decomposition, VMD)，該方法也是一種自适應信号處理方法，該方法假設每個IMF分量是具有不同中心頻率的有限帶寬，為使得每個IMF的估計帶寬之和最小，通過轉換解決變分問題，将各IMF解調到相應的基頻帶，最終提取各個IMF及相應的中心頻率。自VMD方法提出以來，在機械故障檢測[27-29]、電力負荷預測[30-31]及語音去噪[32-33]等領域廣泛應用。但是，變分模态分解的模态個數目前一般是經驗值，同時鮮有學者将其應用在GNSS坐标時間序列分析之中[34]。考慮到CMONOC測站遍布大陸大陸，測站所在地的地形地貌、氣候、土壤及水文條件等地理環境因素差異巨大，VMD方法在提取測站時變信号中的有效性和普适性需要進行全面分析。本文将利用VMD和MLE兩種方法探測CMONOC測站高程時間序列中的時變信号，驗證VMD在GNSS坐标時間序列時變信号提取中的效果，這對于更加真實地反映大陸GNSS測站的運動特征，進而研究引起測站運動的地球實體起源具有十分重要的意義。
1 理論方法
1.1 GNSS坐标時間序列分析方法GNSS測站坐标時間序列的函數模型表示為[16]
(1)式中，x(t)表示GNSS測站坐标時間序列；等式右邊第1項表示多項式，系數和階次為pi和i，對于大多數測站，一階線性項是由地殼闆塊運動引起的；第2項表示階躍，是指測站坐标時序中的突變，H表示Heaviside階梯函數，bj和tj表示階躍的大小和發生的時刻；第3項代表周期信号，Ak、ωk和φk分别表示第k個周期信号的振幅、角頻率和初相，nA指周期信号的數目；第4、5項用于描述測站的震後形變位移，al、Tl和tl表示對數函數的振幅、弛豫時間和震後形變開始時刻，am、Tm和tm表示指數函數的振幅、弛豫時間和震後形變開始時刻，ε為拟合殘差。在CMONOC測站的資料預處理過程中，已剔除了含有明顯階躍信号和震後形變的測站，是以式(1)中第2、4、5項的影響可暫不予考慮。在正常進行中，多項式項的階數采用常用的一階項，而在各類周期信号中，周年和半周年信号最為顯著，是以式(1)可簡化為
(2)式中，a為線性項的截距；b表示線性速度。
1.2 最大似然估計法設單站、單分量GNSS坐标時間序列滿足
(3)式中，p(1)、p(2)分别表示測站的初始坐标和線性速率；q-1為周期信号的數量；p(2k－1)、p(2k)為調和函數的系數，表示周期信号的振幅；ε表示拟合殘差。矩陣形式表示為
(4)若僅考慮周年和半周年信号，則q=3。矩陣x、A和p的具體形式表示為
(5)假設式(4)中的坐标殘差ε是由白噪聲α和有色噪聲β(t)的線性組合
(6)式中，系數a、bk分别表示白噪聲和有色噪聲的振幅大小。通過最小二乘估計法可得到最小二乘解為
(7)協方差矩陣Cp可以表達若幹随機噪聲過程，如白噪聲、可變白噪聲、閃爍噪聲、随機遊走噪聲、幂律噪聲、一階高斯馬爾可夫噪聲，以及它們之間的各類組合。對不同的噪聲模型組合，Cp可以表示為不同形式。通過調整Cp使得似然函數取得最大值，即可得到與該時間序列最相近的噪聲模型。假設觀測值服從高斯正态分布，極大似然函數可表示為
(8)式中，det是矩陣的行列式；N為GNSS坐标時間序列的曆元數目。按照MLE原理，通過選擇不同的噪聲模型，使得殘差
與其協方差陣Cp的聯合機率密度達到最大，數值越大，結果越可靠。聯合機率密度函數的對數形式為
(9)MLE可以處理非均勻采樣，以及資料缺失較多的時間序列，而且可以同時估計測站的線性速度、周期項、階躍及噪聲振幅等參數。
1.3 變分模态分解不同于EMD方法，在VMD方法中定義IMF為一個調幅調頻信号uk(t)，其表達式為
(10)式中，Ak(t)為瞬時幅值；φk(t)為瞬時相位，對其求導得瞬時頻率ωk(t)
(11)VMD算法将變分問題描述為在限制條件為各IMF之和等于輸入信号x下尋求k個IMF，使得所有IMF的估計帶寬之和最小。具體構造步驟如下。(1) 對每個分量uk(t)進行Hilbert變換，得到單側頻譜
(12)(2) 通過加入指數項調整各IMF的中心頻率，将各IMF的頻譜調制到相應的基頻帶
(13)(3) 計算上述解調信号梯度的L2範數，估計出各IMF的帶寬，假定将原始信号x(t)分解為N個IMF分量，則對應的限制變分模型表達式為
(14)式中，∂/∂t表示求偏導；δ(t)為沖激函數；*表示卷積運算；uk={u1, u2, …, uN}代表原始信号x(t)分解得到的N個IMF分量；ωk={ω1, ω2, …, ωN}表示各IMF分量的中心頻率。為求取限制變分模型的最優解，VMD方法通過引入拉格朗日乘子λ(t)和二次懲罰因子α，将待求解的限制性變分問題轉變為非限制性變分問題。增廣拉格朗日表達式為
(15)式中，〈·〉表示點積運算，各模态分量uk和對應的中心頻率ωk采用乘法算子交替方向法求解[35]。IMF分量個數k的确定在VMD方法中十分重要，k值過小導緻模态欠分解，即重構信号的分量個數不足，會影響後續重構過程的精度；k值過大會導緻模态重複或産生額外的噪聲。在機械故障診斷的應用中，一般選擇VMD分量的個數為3~8，利用中心頻率觀察法确定。首先分别設定k值為3~8依次進行試驗，估計不同k值時每個分量的中心頻率，标記第一次出現中心頻率相近的k值，則确定k-1為分解模态個數。
1.4 時變信号的提取精度為了說明VMD方法提取GNSS坐标時間序列中時變信号的精度，使用RMSE改進率進行判定
(16)式中，RMSEimprove為RMSE的改進率；RMSEreduct1為VMD方法的RMSE減小率；RMSEreduct2為參數化方法求得的RMSE減小率；x表示原始序列；ε表示MLE方法得到的殘差序列，與式(2)中含義一緻；ξ表示VMD方法得到的殘差序列，其中RMSE表達式為
(17)式中，s為坐标序列列向量；n為序列的曆元數目。RMSEimprove展現了VMD方法相對于MLE方法在降低原始信号RMSE值的有效性方面的差異，值越大，改進效果越好。除了上述定義的RMSE改進率，還可以利用原始序列x(t)與重構序列(t)之間的相關系數ρ和信噪比SNR定量評價時變信号的提取效果[36-37]，它們定義如下。(1) 原始序列與重構序列之間的相關系數ρ表示為
(18)式中，Cov[x(t), (t)]為原始序列x(t)與重構序列(t)之間協方差；D[x(t)]、D[(t)]分别為x(t)與(t)的方差；n為序列長度。ρ值越接近1，表明原始信号與重構信号的相似度越高，時變信号的提取精度越高。(2) 信噪比SNR的計算公式表示為
(19)SNR主要是展現噪聲信号在原始信号中的比重，其值越大，降噪效果越好。
1.5 陸态網高程時間序列本文使用的GNSS坐标時序資料由中國地震局(China Earthquake Administration, CEA)提供，包括所有CMONOC測站的原始的(raw)和去趨勢項的(detrend)GNSS坐标時間序列，時間跨度為1999—2019年。中國地震局采用GAMIT/GLOBK10.4軟體[38]對所有CMONOC基準站進行了統一解算，得到了基準站的單日松弛解結果，然後經過七參數相似變換将松弛解轉換至ITRF參考架構中[30]，具體的資料處理政策及說明參考 ftp://ftp.cgps.ac.cn/doc/processing_manual.pdf。GNSS坐标時間序列通常包含了一定量的異常值，以及由于地震或裝置更換引起的坐标階躍，是以應先采用一系列的預處理手段擷取相對幹淨的坐标時間序列。首先，剔除形式誤差在水準分量(北向、東向)大于10 mm，同時在高程分量大于20 mm的曆元。然後，采用CEA識别出的階躍值來校正測站坐标時間序列，計算方法是将完整的時間序列在階躍時刻分為前後兩段分别進行函數拟合，前後兩段序列在階躍時刻的估計值之差即為階躍值。對于某些無法識别的不明原因的階躍，有必要對測站逐一檢查，以避免未扣除的階躍信号對測站長期趨勢和周期項參數估值的不利影響。此外，有些測站還表現出了異常的非線性形變。如果這種異常非線性信号産生的原因是已知的，并且異常僅持續了較短的時間，就簡單地剔除測站的異常時間段；如果異常未知并持續時間較長，則删除該站，在下一步的計算分析中不予考慮。在完成資料預處理并删除有效時間跨度少于3 a的測站或存在未知的異常非線性形變的測站之後，剩下的243個CMONOC測站是本文研究的對象。
2 結果與分析
在利用VMD對目标信号進行分解的過程中，IMF分量個數的選擇至關重要。是以，需要首先讨論IMF分量個數的選取，在此基礎上，對IMF分量進行重構，分析VMD方法相較于MLE方法在提取測站時變信号精度方面的差異。對于一些名額的極大、極小值對應的測站，有必要單獨進行研究分析。
2.1 IMF分量個數的确定考慮到GNSS測站位置受到多種因素的影響，不僅包含來自衛星端、傳播路徑、接收端的定位誤差，也包含了環境負載效應、熱膨脹效應及測站位置相關的局部地球實體因素等。影響GNSS測站位置的信号源衆多且存在互相耦合的情況，是以，本文選擇VMD分量的個數為3~10，另外懲罰因子α采用預設值1000。确定IMF分量個數一般選擇中心頻率觀察法，具體步驟為：首先，分别設定k=3、4、5、6、7、8、9、10進行試驗；然後，估計不同k值時每個分量的中心頻率，标記第1次出現中心頻率相近的k值；最後，确定k-1為分解模态個數。以測站TJWQ(天津武清)為例，不同k值下(3~10)的中心頻率見表 1，機關：次每年(cpy)。
表 1 測站TJWQ在不同分解層數下IMF分量的中心頻率Tab. 1 Center frequency of IMF components at TJWQ under different decomposition layers cpy
表選項
可以看出，對于k值取3~10的過程中，沒有出現中心頻率相近的模态，是以從避免出現模态混疊的角度看，取IMF分量個數取10是可行的；然而，對于IMF2(主要表現為測站季節性信号)，頻率約為1 cpy，而k≤5時，并沒有出現該頻率的模态；k=6時，第一次出現該模态，是以當顧及實際信号的季節性特征，且分解模态數盡可能少的原則，測站TJWQ應該取的模态個數為6。k=6時，該站原始序列、各IMF分量、殘差序列及相應的頻譜如圖 1所示。
| |
| 圖 1 測站TJWQ高程式列、IMF分量、殘差序列及相應的頻譜Fig. 1 Height time series, IMF components, residual series and corresponding power spectrum of station TJWQ |
| 圖選項 |
可以看出，信号分量從IMF1到IMF6由低頻到高頻變化，IMF1包含了測站的長期線性項，IMF2則反映了測站的季節性運動特征，IMF3—IMF6表示頻率更高、頻譜強度更弱的某種信号。當選擇IMF1和IMF2之和作為重構信号，根據式(16)計算出的RMSE改進率為1.6%。是以，相對于最大似然估計法，VMD具有更優越的減弱測站高程時間序列非線性形變的性能。可見，VMD對複雜信号的處理能力，能夠将多分量信号一次性分解成多個單分量調幅調頻信号，是一種高效的處理非線性信号的工具。然而，VMD算法本身也存在一些不足之處。對于不同測站，受到的地球實體信号源的影響不同，VMD分解層數也不同。對于CMONOC測站而言，對所有測站逐一使用中心頻率觀察法進行參數确定任務量較大。同時對于懲罰因子本文采用了常用值，并未給予過多讨論分析。是以，後續研究還需要針對具體測站進行VMD分解層數和懲罰因子的确定，計劃采用鲸魚優化算法(whale optimization algorithm, WOA)自适應确定IMF分量個數。
2.2 CMONOC測站時變信号提取精度VMD分解過程中，IMF分量個數一般選擇為3~8。信号重構時選擇的IMF分量越多，越能還原出原始信号，與原始信号的相似性就越強，相較于最大似然估計法的RMSE改進率越大，但同時重構信号的組成也會越複雜。是以選擇合适的IMF分量至關重要。本節将從RMSE改進率、相關系數和信噪比這3個方面分析VMD方法提取CMONOC測站高程時變信号的精度，重構分量的個數選擇3個。
2.2.1 RMSE改進率利用式(9)，計算得到VMD方法相對于MLE方法的RMSE改進率，CMONOC測站的RMSE改進率的空間分布與分布直方圖如圖 2、圖 3所示。全部243個測站的RMSE改進率均為正值，表明VMD方法有助于所有測站提取時變信号，減弱高程時間序列中的非線性形變。全部測站RMSE改進率的均值為69.8%，RMSE改進率在60%以上的測站數為193個，占比79.4%，超過100%的測站有2個，分别是YONG(104.6%)和GXBH(100.2%)，最低的測站是TJWQ(7.8%)。這些測站将在2.3節中詳細分析。RMSE改進率并無明顯的地理分布特征。
| |
| 圖 2 CMONOC測站RMSE改進率空間分布Fig. 2 Spatial distribution of RMSE improvement rate at CMONOC stations |
| 圖選項 |
| |
| 圖 3 RMSE改進率分布直方圖Fig. 3 Distribution histogram of RMSE improvement rates |
| 圖選項 |
2.2.2 相關系數利用式(10)，計算得到VMD提取出的重構信号與原始信号之間的相關系數，選取前3個IMF分量進行信号重構時，原始信号與重構信号之間的相關系數達0.99~1，表明原始信号與重構信号的相似度極強，測站的時變信号能夠通過VMD方法完美地提取出來。由于所有IMF分量之和即為原始信号，是以重構信号包含的IMF越多，其與原始信号之間的相關性越強。若隻選取IMF1進行信号重構時，重構信号與原始信号之間的相關系數降低，此時相關系數超過0.9的測站數為40個。圖 4表示CMONOC測站的VMD重構信号與原始信号之間的相關系數，可以看出，這些測站較為均勻地分布在大陸大陸，并未表現出明顯的地理分布特征。位于天津、河北等沿海地區的測站相關系數較高。相對于VMD方法，MLE方法得到的拟合信号與原始信号之間的相關系數明顯較低，全部測站相關系數的均值為0.62，兩種方法得到的相關系數分布直方圖如圖 5所示。
| |
| 圖 4 CMONOC測站VMD重構信号與原始信号的相關系數空間分布Fig. 4 Spatial distribution of correlation coefficients between reconstructed signals derived from VMD method and original signals at CMONOC stations |
| 圖選項 |
| |
| 圖 5 相關系數分布直方圖Fig. 5 Distribution histogram of correlation coefficients |
| 圖選項 |
2.2.3 信噪比利用式(11)，計算得到VMD方法相對于MLE方法的信噪比，信噪比的空間分布與分布直方圖如圖 6、圖 7所示。由式(11)可知，原始信号的信噪比均為0，經過VMD提取出時變信号以後，信噪比提高，全部243個測站的信噪比均值為28.0 dB，有7個測站的信噪比達到40 dB以上，其中信噪比最大的測站是HECX(53.8%)。相對于VMD方法，MLE方法得到的拟合信号的信噪比較低，全部測站的信噪比均值僅為3.1 dB。
| |
| 圖 6 CMONOC測站VMD方法信噪比空間分布Fig. 6 Spatial distribution of SNR derived from VMD method at CMONOC stations |
| 圖選項 |
| |
| 圖 7 信噪比分布直方圖Fig. 7 Distribution histogram of SNR |
| 圖選項 |
2.3 特殊測站分析上節提到，全部CMONOC測站中有18個測站的RMSE改進率大于90%，其中測站GXBH、HBES的RMSE改進率分别為100.2%和91.9%，它們的原始信号、MLE拟合信号與VMD重構信号如圖 8所示。
| |
| 圖 8 RMSE改進率大于90%的2個CMONOC測站Fig. 8 Two CMONOC stations with RMSE improvement rates greater than 90% |
| 圖選項 |
測站GXBH(廣西北海)經過VMD方法提取時變信号，RMSE改進率達100.2%。從原始序列可以看出，該站并未表現出明顯的季節性信号，即整年周期運動，是以基于“周年+半周年”假設的MLE方法得到的拟合序列并未出現年度振蕩(綠點)，拟合效果不佳。相對而言，VMD方法首先将原始信号分解為不同中心頻率的IMF分量，然後選取前3個IMF分量進行信号重構，重構序列能夠較好拟合出原始序列中的高頻信号。測站HBES(湖北恩施)RMSE改進率為91.9%。從原始序列(灰點)可以看出，該測站在2014—2015年間出現一段異于其他年份的下沉現象，VMD方法得到的重構序列可以較好反映這一情況，而MLE方法得到固定的周年振幅和相位，與該時間段内的實際位置變化差異較大。是以，通過VMD方法的分解與重構過程，測站非線性形變可以得到較好的提取。GNSS坐标時間序列的預處理是一項重要的工作，然而在測站數目衆多時，階躍信号的準确探測和識别會變得相當複雜。對于已知原因的階躍信号，可以根據地震或儀器更換發生時刻估算階躍值的大小，進而對原始序列進行改正；對于原因未知或發生時刻難以确定的階躍信号，目前最有效的方法仍然是人工目視檢查[40]，但不可避免會出現遺漏情況。如圖 9所示，測站TASH位于新疆塔什庫爾幹塔吉克自治縣，該測站存在較長時間的資料缺失和較大的階躍，當遺漏該階躍信号時，MLE方法會錯誤估計原始序列，而VMD方法仍然能夠較準确地得到反映測站實際變化的重構序列，因而使得RMSE改進率達到79.8%。同樣地，位于南海永興島的測站YONG，存在較長時間的資料缺失和階躍現象，MLE方法不能較好反映原始序列中的高頻信号，使得RMSE改進率在所有CMONOC測站中最大，達到了104.6%。
| |
| 圖 9 存在階躍信号的2個CMONOC測站Fig. 9 Two CMONOC stations with offset signal |
| 圖選項 |
全部CMONOC測站中，有10個測站的RMSE改進率低于20%，其中最低的3個測站為TJWQ(7.8%)、TJBH(11.6%)和HECX(11.7%)，它們的原始信号、拟合信号與重構信号分别如圖 10所示。分析發現，這些測站均存在顯著的趨勢項。
| |
| 圖 10 RMSE改進率小于20%的3個CMONOC測站Fig. 10 Three CMONOC stations with RMSE improvement rates smaller than 20% |
| 圖選項 |
測站TJWQ(天津武清)由于近十年來華北平原地區地下水開采嚴重，導緻地下水位逐年下降，目前該地區已成為世界上最大的地下水漏鬥區[41-43]，該測站位于漏鬥中心區域。可以發現，該測站地下水下降造成了地表位置的逐年下降，2011—2014年，測站位置呈線性下降，而2014—2019年，測站位置在每年約6、7月出現明顯的回升。對于這種整體呈線性趨勢的原始序列，MLE方法與VMD方法均能夠較好提取出測站中的線性項，相對而言，VMD能夠更多地顧及測站的年度時變信号，是以RMSE改進率為7.8%。測站TJBH(天津濱海)、HECX(河北滄縣)與測站TJWQ情況是一緻的，其中HECX的信噪比在全部CMONOC測站中最大。
3 結語
針對現有參數化方法不能有效拟合GNSS高程時間序列，以及一些非參數化方法(WD、EMD)在信号分解中存在的自适應性較差、模态混疊等缺點，本文将VMD方法應用于中國CMONOC測站高程時間序列的分析中，得出以下結論。(1) 相對于MLE方法，VMD方法全部測站的RMSE改進率為正值，全部測站RMSE改進率的均值為69.8%，該方法有助于GNSS測站提取時變信号，減弱高程時間序列中的非線性形變，驗證了VMD方法在CMONOC測站高程時間序列時變信号提取中的有效性。(2) 當重構分量個數取3時，VMD方法得到的重構序列與原始序列之間的相關系數在全部CMONOC測站的均值近似為1，比MLE方法的相應結果提高了0.38；VMD方法得到的重構序列的信噪比在全部測站的均值為28.04 dB，比MLE方法的相應結果提高了24.9 dB。(3) 雖然VMD方法具有較好的自适應性，但是IMF分量個數仍然與具體的每一個測站有關，當分解個數和重構分量選取恰當時，VMD方法在GNSS高程時間序列中的應用效果可進一步提高。
作者簡介第一作者簡介：武曙光(1992-), 男, 博士, 講師, 研究方向為地學資料處理與信号分析。E-mail: [email protected]通信作者：李厚樸 E-mail：[email protected]
初審：張豔玲複審：宋啟凡
終審：金 君
資訊