小波自回歸模型在天頂對(duì)流層延遲預(yù)測(cè)中的應(yīng)用

李偉捷 ，劉根友

（1. 中科院測(cè)量與地球物理研究所 大地測(cè)量與地球動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢 430077；2. 中國科學(xué)院大學(xué)，北京 100049）

0 引言

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

1 基本理論

1.1 諧波函數(shù)擬合

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

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

1.2 離散小波變換

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

圖1 離散小波變換原理

1.3 自回歸模型

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

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

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

本文結(jié)合諧波函數(shù)擬合、離散小波變換以及自回歸模型3種手段，通過對(duì)已知數(shù)據(jù)的建模，實(shí)現(xiàn)對(duì)流層延遲的預(yù)測(cè)。具體步驟如下：

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

圖2 WAMIX建模流程

2 實(shí)驗(yàn)與結(jié)果分析

2.1 數(shù)據(jù)來源

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

圖3 選用測(cè)站分布

2.2 方法與結(jié)果

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

圖4 ARTU及BJNM站ZTD時(shí)間序列及其諧波函數(shù)擬合圖

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

圖7 ARTU及BJNM站低頻分量AR模型預(yù)測(cè)結(jié)果

圖 5 至圖 7 的（a）、（b）和(c)、(d)分別對(duì)應(yīng) ARTU和BJNM 2014年第1～150天a4低頻分量預(yù)測(cè)結(jié)果，（a）、（c）為樣本區(qū)內(nèi)一步預(yù)測(cè)檢驗(yàn)，（b）、（d）為樣本區(qū)外預(yù)測(cè)結(jié)果。

為驗(yàn)證模型預(yù)測(cè)精度，本文將 WAMIX模型與目前認(rèn)可精度較高的GPT2w(1°×1°)和IGGtrop_SH進(jìn)行精度對(duì)比，以IGS站提供的ZTD為真值，計(jì)算各模型序列殘差，統(tǒng)計(jì)其均方根誤差(root mean square， RMS)和平均相對(duì)中誤差（average absolute relative error， AARD）。據(jù)預(yù)測(cè)結(jié)果統(tǒng)計(jì)（表 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 %。前二者精度相當(dāng)，WAXMIX較之精度有所提高，RMS平均改善6 mm，AARD改善0.2 %。4組預(yù)測(cè)數(shù)據(jù)中，最佳結(jié)果是2016年8月LHAZ預(yù)測(cè)，較 GPT2w模型 RMS改善 28.7 mm，AARD改善1.5 %，較 IGGtrop_SH模型RMS改善20.2 mm，AARD改善1.1 %。

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

表1 2014年6月預(yù)測(cè)精度統(tǒng)計(jì)

表2 2014年8月預(yù)測(cè)精度統(tǒng)計(jì)

表3 2015年6月預(yù)測(cè)精度統(tǒng)計(jì)

表4 2015年8月預(yù)測(cè)精度統(tǒng)計(jì)

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

3 結(jié)束語

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