小波自回歸模型在天頂對流層延遲預測中的應用

李偉捷 ，劉根友

（1. 中科院測量與地球物理研究所 大地測量與地球動力學國家重點實驗室，武漢 430077；2. 中國科學院大學，北京 100049）

0 引言

全球衛星導航系統（global navigation satellite system， GNSS）信號在對流層的傳播過程中，會發生路徑彎曲以及傳播速率改變，導致信號接收端測距誤差，即對流層延遲。對流層延遲誤差與水汽、溫度等因素密切相關，隨衛星信號入射高度角減小而增大，在天頂方向約2.3 m，當衛星高度角降低到10°時，可達15～20 m[1-2]。對流層延遲是 GNSS衛星定位重要誤差之一，提高對流層改正精度有助于模糊度固定和定位精度的提高。一般情況下，對流層誤差可以用模型公式計算，傳統的 Hopfield模型及Saastamonien模型需要實測氣象參數的支持才能滿足實際需求[3-4]，為擺脫對實測氣象參數依賴，許多經驗模型被相繼提出[5-8]。其中，GPT2w 、IGGtrop_SH和GZTDS這 3種對流層延遲經驗模型精度較高，內符合精度分別為3.6、3.8、3.7 cm[9-11]。由于天頂對流層延遲（zenith tropospheric delay，ZTD）及相關氣象因子顯著的年及半年周期特性，這幾種模型的參數均是由頻率f1=1/365.25和f2=2/365.25的諧波函數表示，這種表示的方法簡單有效，有利于全球格網的建立。同時，文獻[12]基于小波分析發現，對流層延遲的噪聲主要集中在低階高頻部分，而通過小波提取的低頻部分能較好地保留有用信息。為了提高單站對流層誤差預測精度，本文在考慮周年和半年諧波項的同時，還將利用小波分析提取剩余殘差的低頻部分，并對相對平穩的低頻信號建立自回歸模型。

1 基本理論

1.1 諧波函數擬合

由于單站 ZTD序列具有顯著的年周期及半年周期，通常通過頻率為f1=1/365.25和f2=2/365.25的諧波函數對單站進行擬合[8]，其計算方法為

式中：t為簡化儒略日，即與1980-01-06零時即全球定位系統（global positioning system，GPS)時零點對應儒略日（2 444 244.5）的差值；X(t)為t時刻的對應ZTD的值；系數Ki通過最小二乘法確定。

1.2 離散小波變換

離散小波變換是通過低通和高通2個互補的濾波器將信號分解成a和d2個部分，a為信號的近似值（approximations），d為信號的細節值（detail）[13]。在地球物理鄰域，該法常用于分析含有多尺度特征、奇異值檢測以及瞬態現象的非平穩信號[14]。本文采用小波分解樹形式，只對低頻信號進行繼續分解，對應公式為原始信號，j為分解階數。其中，圖1為3階分解，則有S=a3+d3+ d2+d1。

圖1 離散小波變換原理

1.3 自回歸模型

自回歸模型（auto regressive model，AR）模型是將t時刻輸出表達成過去時刻輸出序列的線性組合[15]，是一種利用序列自相關特性建立的線性預測模型，模型階數n根據貝葉斯信息準則BIC及系數顯著性確定[16]，其計算方法為

式中：Xt為t時刻的值；?t為擬合殘差；系數Pj通過最小二乘估計；BIC為貝葉斯信息值；n為模型階數；б2為擬合殘差方差；N為樣本數。

1.4 小波及自回歸組合模型（wavelet-autoregressive model, WAMIX）實施流程

本文結合諧波函數擬合、離散小波變換以及自回歸模型3種手段，通過對已知數據的建模，實現對流層延遲的預測。具體步驟如下：

采用諧波式（1）擬合，提取ZTD的主要趨勢項（年周期及半年周期信號）；通過小波把上一步殘差結果分解成低頻及高頻分量，挖掘各分量的變化特性，選擇合適階數，將高頻分量作為噪聲剔除并保留低頻分量；對保留的低頻信號建立適宜的 AR模型，實現時序預測。具體流程如圖2所示。

圖2 WAMIX建模流程

2 實驗與結果分析

2.1 數據來源

對流層延遲數據采自國際 GNSS服務組織（International GNSS Service，IGS）數據共享平臺ftp：//cddis.gsfc.nasa.gov/，該平臺提供了高精度的天頂對流層延遲ZTD，采樣間隔為300 s。本文選取了2013年1月1日至2016年12月31日數據完整率達90 %的ARTU、BJNM、DRAO、GODE、LHAZ、 KOUR、YELL等7個測站作為研究對象。選用每天中午12時的數值為當日值，并通過樣條插值得到等間隔為1 d的數據，測站分布如圖3所示。

圖3 選用測站分布

2.2 方法與結果

首先利用2013—2016年全時段數據，通過快速傅里葉變換(fast Fourier transform，FFT)分析存在明顯的周年項和半年項，因此可以根據式（1）提取出各站ZTD序列中主要趨勢項年周期及半年周期，如圖4（鑒于篇幅，僅列出ARTU和BJNM站）所示。然后對觀測值減去諧波擬合后的殘差進行 db45小波 10階分解。文獻[12]中提到，ZTD小波變換前幾層高頻分量主要是由測量、儀器、天氣或人為原因等引起的誤差，因此，本文參考圖5的小波信號分解，保留相對平穩的a4分量進行AR建模，其中a1～a10表示 10 階分解得到的低頻分量，d1～d10對應的是高頻分量，a4中包含的為 16 d及更長時段對應低頻信號。建模對象為4組150 d連續樣本（2014和 2015年年積日第 1～150天及年積日第 51～200天），根據 BIC準則及系數顯著性，AR模型定階為11。最后，對樣本區數據進行一步預測檢驗并采用通過檢驗的AR[11]模型預測后30 d的結果，與觀測值進行比較，統計其精度，其結果見圖5至圖7。

圖4 ARTU及BJNM站ZTD時間序列及其諧波函數擬合圖

圖6 ARTU及BJNM站諧波擬合殘差及其傅里葉變換頻譜分析結果

圖7 ARTU及BJNM站低頻分量AR模型預測結果

圖 5 至圖 7 的（a）、（b）和(c)、(d)分別對應 ARTU和BJNM 2014年第1～150天a4低頻分量預測結果，（a）、（c）為樣本區內一步預測檢驗，（b）、（d）為樣本區外預測結果。

為驗證模型預測精度，本文將 WAMIX模型與目前認可精度較高的GPT2w(1°×1°)和IGGtrop_SH進行精度對比，以IGS站提供的ZTD為真值，計算各模型序列殘差，統計其均方根誤差(root mean square， RMS)和平均相對中誤差（average absolute relative error， AARD）。據預測結果統計（表 1至表4），GPT2w 的 RMS 為 36.7 mm， AARD 為 1.3 %；IGGtrop_SH的RMS為 36.6 mm， AARD為 1.3 %；WAMIX的 RMS為 30.5 mm， AARD為 1.1 %。前二者精度相當，WAXMIX較之精度有所提高，RMS平均改善6 mm，AARD改善0.2 %。4組預測數據中，最佳結果是2016年8月LHAZ預測，較 GPT2w模型 RMS改善 28.7 mm，AARD改善1.5 %，較 IGGtrop_SH模型RMS改善20.2 mm，AARD改善1.1 %。

對比 2014年及 2015年的預測結果時發現，2014年 WAMIX預測精度低于 2015年，尤其是在2014年6月，GODE站及KOUR站的RMS較GPT2w及IGGtrop_SH分別大了7和2 mm。為分析誤差來源，本文還納入了式（1）諧波擬合的精度結果，發現這2站WAMIX預測結果不及諧波，同時，在ARTU及BJNM站出現類似問題，但在圖6（b）、圖6（d）中預測值和樣本低頻分量a4基本吻合。因此，本文推測這不是 AR模型本身擬合的問題，而是由于建模時未考慮高頻分量造成的。

表1 2014年6月預測精度統計

表2 2014年8月預測精度統計

表3 2015年6月預測精度統計

表4 2015年8月預測精度統計

高頻信號中除噪聲外，還包含 ZTD真實信號受氣象影響產生的高頻分量，其主要影響因素是水汽，溫度以及氣壓[8]。據本文獲取的氣象資料記載，北京2014年6月氣候多變且晝夜溫差較大，陸續出現了16 d降雨，其中含15 d雷陣雨。由此產生的水汽和溫度的強烈變化會直接導致ZTD不穩定，產生大量高頻信號，對前文的推測予以佐證。相反，精度效果改善最佳的 LHAZ站，地處中國西部降雨偏少且海拔較高，高頻信號相對穩定。由于缺乏其他站的氣象數據，具體的產生原因還有待進一步驗證。

3 結束語

本文基于IGS 7個測站4 a（2013—2016年）的ZTD數據，采用一種小波-AR組合模型WAMIX實現單站建模預測，采用全時段數據進行諧波擬合和小波分析，在此基礎上利用150 d的低頻分量預測后 30 d的數據。WAMIX的預測精度為：RMS：3 cm， AARD：1.1 %，較 GPT2w和IGGtrop_SH模型整體精度為好，RMS平均改善6 mm，AARD改善0.2 %，RMS最優改善2 cm，AARD最優改善1 %。在氣象穩定區域，ZTD真實高頻分量較小，模型效果較好；在天氣多端甚至異常區域，高頻分量不能完全作為噪聲處理，模型效果可能不佳，需要綜合其他氣象資料進一步研究。