陳 祥 楊志強(qiáng) 楊 兵 楊偉華

1 長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，西安市雁塔路126號(hào)，710054

全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system, GNSS)作為空間大地測(cè)量的主要技術(shù)手段之一，已被廣泛應(yīng)用于地球內(nèi)部動(dòng)力學(xué)機(jī)制、板塊運(yùn)動(dòng)和地殼形變監(jiān)測(cè)等研究[1]。目前，中國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)(陸態(tài)網(wǎng)絡(luò))已建成由近300個(gè)GNSS基準(zhǔn)站和2 000多個(gè)GNSS流動(dòng)站構(gòu)成的觀測(cè)網(wǎng)絡(luò)[2]。由于受數(shù)據(jù)處理策略、鐘差、電離層延遲、對(duì)流層延遲和多路徑效應(yīng)等因素的影響，解算得到的GNSS垂向坐標(biāo)時(shí)間序列存在誤差，從而影響GNSS定位的精度和可靠性，導(dǎo)致某些地球物理現(xiàn)象出現(xiàn)錯(cuò)誤解釋[3]。因此，對(duì)GNSS垂向坐標(biāo)時(shí)間序列進(jìn)行降噪研究具有非常重要的科學(xué)意義和現(xiàn)實(shí)意義[4]。

目前用于GNSS坐標(biāo)時(shí)間序列分解的方法主要有小波[5]、經(jīng)驗(yàn)?zāi)B(tài)分解[6]及STL(seasonal-trend decomposition procedure based on LOESS)[7]等，其中小波和經(jīng)驗(yàn)?zāi)B(tài)分解方法側(cè)重于去噪分析，而STL方法主要用于季節(jié)性變化的提取，較少涉及去噪效果分析[8]。本文采用LOESS方法對(duì)GNSS垂向坐標(biāo)時(shí)間序列進(jìn)行降噪分析，并對(duì)其有效性及降噪效果進(jìn)行驗(yàn)證和評(píng)價(jià)。

1 LOESS原理及降噪方法

1.1 LOESS方法基本原理

LOESS是一種基于局部回歸分析的非參數(shù)降噪方法[7]，與線性回歸、多項(xiàng)式回歸等方法相比，該方法未對(duì)預(yù)測(cè)點(diǎn)進(jìn)行限制，而是直接基于數(shù)據(jù)進(jìn)行分析。該方法需要預(yù)先設(shè)置一個(gè)窗口，根據(jù)窗口內(nèi)的點(diǎn)計(jì)算擬合點(diǎn)的擬合值，具體過(guò)程如下：

2)設(shè)(xk,yk)為窗口在最左端時(shí)的中心點(diǎn)，重復(fù)步驟1)，對(duì)(x1,y1)到(xk,yk)進(jìn)行擬合，此時(shí)窗口位置固定，但從點(diǎn)(xk+1,yk+1)開始需要以擬合點(diǎn)為中心移動(dòng)窗口位置，并逐個(gè)計(jì)算各點(diǎn)的擬合值，參與擬合值計(jì)算的仍為窗口內(nèi)所有點(diǎn)。

3)當(dāng)窗口移動(dòng)到坐標(biāo)時(shí)間序列的最右端時(shí)，與最左端的情況類似，此時(shí)窗口中心點(diǎn)(xm-k+1,ym-k+1)到右端點(diǎn)(xm,ym)(m為坐標(biāo)時(shí)間序列的總長(zhǎng)度)均采用同一窗口內(nèi)的所有點(diǎn)進(jìn)行擬合。

4)通過(guò)上述步驟對(duì)坐標(biāo)時(shí)間序列中所有點(diǎn)進(jìn)行擬合，將擬合后的點(diǎn)連接成完整的LOESS回歸曲線，該曲線即為降噪后的坐標(biāo)時(shí)間序列。

該方法的關(guān)鍵在于設(shè)置窗口寬度，經(jīng)過(guò)多次實(shí)驗(yàn)，本文將窗寬設(shè)置為75。關(guān)于窗口內(nèi)所有點(diǎn)權(quán)重的計(jì)算，可通過(guò)窗口內(nèi)各點(diǎn)與擬合點(diǎn)的距離進(jìn)行確定，該距離指X軸方向上的距離。隨著距離的增加，權(quán)重減小，因此需要采用三次權(quán)重函數(shù)進(jìn)行轉(zhuǎn)化：

(1)

此時(shí)還需要確定距離擬合點(diǎn)最遠(yuǎn)的點(diǎn)，并計(jì)算最大距離，對(duì)其他距離進(jìn)行歸一化處理，從而確定各點(diǎn)的權(quán)重：

(2)

式中，x為擬合點(diǎn)橫坐標(biāo)，xi為窗口內(nèi)某點(diǎn)的橫坐標(biāo)，Δ(x)為最遠(yuǎn)點(diǎn)與擬合點(diǎn)之間的距離，υi(x)為該點(diǎn)相對(duì)于x的權(quán)重。

當(dāng)目標(biāo)模型為非線性模型時(shí)，局部加權(quán)回歸算法相對(duì)于線性回歸算法更合適，該算法選擇窗口內(nèi)的點(diǎn)而不是全部點(diǎn)進(jìn)行回歸。加權(quán)最小二乘法的目標(biāo)函數(shù)為：

(3)

式中，J(θ)為損失函數(shù)，(xi,yi)為窗口內(nèi)某點(diǎn)，θ為最佳回歸系數(shù)。

1.2 技術(shù)路線

用LOESS方法對(duì)陸態(tài)網(wǎng)絡(luò)289個(gè)測(cè)站的垂向坐標(biāo)時(shí)間序列進(jìn)行降噪分析，具體流程如下：1)降噪前先剔除坐標(biāo)時(shí)間序列中的粗差；2)用LOESS方法對(duì)剔除粗差后的序列進(jìn)行降噪處理，得到降噪后的坐標(biāo)時(shí)間序列及噪聲序列；3)確定噪聲模型，并利用極大似然估計(jì)(MLE)方法解算GNSS坐標(biāo)時(shí)間序列諧波模型對(duì)應(yīng)的各項(xiàng)運(yùn)動(dòng)參數(shù)；4)檢驗(yàn)降噪后噪聲序列的自相關(guān)性和降噪前后各指標(biāo)的顯著性；5)綜合數(shù)學(xué)指標(biāo)和GNSS測(cè)站運(yùn)動(dòng)模型參數(shù)對(duì)降噪結(jié)果進(jìn)行評(píng)價(jià)。技術(shù)路線如圖1所示。

圖1 技術(shù)路線

2 數(shù) 據(jù)

2.1 數(shù)據(jù)來(lái)源

考慮到坐標(biāo)時(shí)間序列數(shù)據(jù)缺失的問(wèn)題，本文采用陸態(tài)網(wǎng)絡(luò)中289個(gè)GNSS基準(zhǔn)站原始的垂向坐標(biāo)時(shí)間序列進(jìn)行降噪分析。GNSS坐標(biāo)時(shí)間序列數(shù)據(jù)來(lái)源于中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)(http:∥www.cgps.ac.cn)，數(shù)據(jù)處理方法與精度指標(biāo)可參考該平臺(tái)提供的說(shuō)明文檔(ftp:∥ftp.cgps.ac.cn/doc/processing_manual.pdf)。

2.2 數(shù)據(jù)預(yù)處理

由于受各種外界因素的影響，原始數(shù)據(jù)中含有一定的粗差，在分析前需要對(duì)序列中的粗差進(jìn)行探測(cè)與剔除。本文采用IQR方法[9]消除GNSS原始坐標(biāo)時(shí)間序列中的異常值。

3 算例及結(jié)果分析

3.1 算例分析

采用§1.1的方法對(duì)選取的289個(gè)GNSS站進(jìn)行降噪處理。為驗(yàn)證LOESS方法對(duì)GNSS坐標(biāo)時(shí)間序列降噪的可靠性，以XJHT站為例進(jìn)行說(shuō)明，圖2為XJHT站原始坐標(biāo)時(shí)間序列、降噪后時(shí)間序列及噪聲序列。由表1可知，XJHT站垂向坐標(biāo)時(shí)間序列降噪后的信噪比達(dá)48.47 dB，降噪效果較好，標(biāo)準(zhǔn)差減小了1.11 mm，白噪聲和閃爍噪聲分別從1.81 mm、11.27 mm·a-0.25降至0.01 mm、2.40 mm·a-0.25，表明序列中包含的噪聲明顯降低，速度不確定度從0.48 mm·a-1降至0.10 mm·a-1。

圖2 XJHT站垂向GPS原始坐標(biāo)時(shí)間序列、降噪后時(shí)間序列和噪聲項(xiàng)

表1 XJHT站降噪前后評(píng)價(jià)指標(biāo)

3.2 自相關(guān)性檢驗(yàn)

為進(jìn)一步分析降噪后噪聲序列是否存在自相關(guān)性，本文采用Durbin-Watson(DW)檢驗(yàn)對(duì)降噪后的噪聲序列進(jìn)行自相關(guān)性檢驗(yàn)。DW檢驗(yàn)是一種檢驗(yàn)序列自相關(guān)性的方法[10]，采用該方法對(duì)各測(cè)站降噪后的噪聲序列進(jìn)行檢驗(yàn)可得到檢驗(yàn)統(tǒng)計(jì)量。由判斷準(zhǔn)則可知，若DW值約為2，即可推斷序列完全不相關(guān)。

采用§1.1的方法對(duì)289個(gè)GNSS站的垂向坐標(biāo)時(shí)間序列進(jìn)行降噪處理后，對(duì)噪聲序列進(jìn)行DW自相關(guān)性檢驗(yàn)，結(jié)果如圖3所示。結(jié)果表明，DW平均值為1.64±0.19(n=289)，其中78%介于1.50～2.10之間，表明采用LOESS方法進(jìn)行降噪處理后，各測(cè)站的噪聲序列不存在自相關(guān)性，可進(jìn)行后續(xù)的評(píng)價(jià)工作。

圖3 DW檢驗(yàn)結(jié)果

3.3 評(píng)價(jià)指標(biāo)

本文采用以下4種指標(biāo)來(lái)定量評(píng)價(jià)LOESS方法對(duì)GNSS坐標(biāo)時(shí)間序列降噪的可靠性和效果：1)信噪比(signal-noise ratio, SNR)；2)標(biāo)準(zhǔn)差(standard deviation, STD)；3)白噪聲(white noise, WN)和閃爍噪聲(flicker noise, FN)；4)速度不確定度。

3.3.1 信噪比

SNR計(jì)算公式可表示為：

(4)

3.3.2 GNSS測(cè)站運(yùn)動(dòng)模型

研究表明[11]，白噪聲(WN)和閃爍噪聲(FN)組成的噪聲模型(WN+FN)是GNSS垂向坐標(biāo)時(shí)間序列的最優(yōu)噪聲模型，因此本文采用WN+FN噪聲模型。

作為常用的GNSS坐標(biāo)時(shí)間序列分析軟件，Hector軟件采用極大似然估計(jì)方法求解GNSS測(cè)站的運(yùn)動(dòng)模型參數(shù)：

(5)

3.3.3 指標(biāo)改正率

參照文獻(xiàn)[12-14]使用的指標(biāo)，本文計(jì)算評(píng)價(jià)指標(biāo)的改正率(correction rate, CR)，以進(jìn)一步分析評(píng)價(jià)LOESS方法對(duì)GNSS垂向坐標(biāo)時(shí)間序列的降噪效果：

(6)

式中，Ibefore和Iafter分別為各評(píng)價(jià)指標(biāo)在時(shí)間序列降噪前后的指標(biāo)值，CR為正表示LOESS方法能夠有效削弱GNSS垂向坐標(biāo)時(shí)間序列中包含的噪聲，CR值越大則表示降噪效果越明顯。

3.4 降噪效果評(píng)價(jià)

采用LOESS方法對(duì)各測(cè)站進(jìn)行降噪，得到降噪后的時(shí)間序列及噪聲部分，并計(jì)算降噪前后SNR、STD、噪聲和速度不確定度等4個(gè)評(píng)價(jià)指標(biāo)值，分析其變化情況。為驗(yàn)證降噪前后序列的各項(xiàng)指標(biāo)是否存在顯著性差異，本文采用非參數(shù)Wilcoxon秩和檢驗(yàn)[15]對(duì)各項(xiàng)指標(biāo)在降噪前后進(jìn)行顯著性檢驗(yàn)，若檢驗(yàn)結(jié)果中P值均小于0.05，則表明在95%的置信區(qū)間內(nèi)，指標(biāo)在降噪前后存在顯著性差異。

3.4.1 SNR與STD

由式(4)可知，原始序列的SNR均為0。圖4為各測(cè)站坐標(biāo)時(shí)間序列降噪后的信噪比，測(cè)站分布直方圖如圖5所示。由圖4和5可知，垂直分量坐標(biāo)時(shí)間序列降噪后的信噪比均較高，最高達(dá)48.47 dB，均值為8.14 dB，80%以上測(cè)站介于0～15 dB之間；部分測(cè)站的信噪比為負(fù)，其原因可能為使用LOESS方法降噪時(shí)設(shè)置的窗寬較小。結(jié)果表明，采用LOESS方法進(jìn)行降噪后信噪比總體有顯著提高。

圖4 各測(cè)站坐標(biāo)時(shí)間序列降噪后的信噪比

圖5 降噪后信噪比的測(cè)站分布

圖6為各測(cè)站垂向坐標(biāo)時(shí)間序列降噪前后的標(biāo)準(zhǔn)差變化，標(biāo)準(zhǔn)差改正率的測(cè)站分布情況如圖7所示。由圖6可知，采用LOESS方法降噪后的STD均有所減小，平均減小21.35%，最高減小48.77%，且Wilcoxon檢驗(yàn)表明，降噪前后標(biāo)準(zhǔn)差具有顯著性差異(p<0.05,n=289)。由圖7可知，22%的測(cè)站改正率介于30%～50%之間，34%的測(cè)站改正率介于20%～30%之間，因此采用LOESS方法能夠顯著降低時(shí)間序列的標(biāo)準(zhǔn)差。

圖6 降噪前后各測(cè)站序列的標(biāo)準(zhǔn)差

圖7 標(biāo)準(zhǔn)差改正率的測(cè)站分布

3.4.2 運(yùn)動(dòng)模型參數(shù)

圖8(a)和8(b)分別為各測(cè)站垂向坐標(biāo)時(shí)間序列降噪前后白噪聲和閃爍噪聲變化情況，2種噪聲改正率的測(cè)站分布情況分別見圖9(a)和9(b)。由圖8(a)可知，降噪后序列的噪聲均明顯降低，特別是白噪聲均降低至0.03 mm以下，且Wilcoxon檢驗(yàn)表明，降噪前后序列的白噪聲具有顯著性差異(p<0.05,n=289)。由圖9(a)可知，94.81%的測(cè)站白噪聲改正率在95%以上。從圖8(b)可以看出，閃爍噪聲降低至0～5 mm·a-0.25范圍內(nèi)，改正率最低為68.35%，最高為85.28%，且Wilcoxon檢驗(yàn)表明，降噪前后序列的閃爍噪聲具有顯著性差異(p<0.05,n=289)。由圖9(b)可知，約34%的測(cè)站閃爍噪聲改正率在80%～86%之間，約64%的測(cè)站在70%～80%之間。結(jié)果表明，LOESS方法可有效降低時(shí)間序列中的噪聲，其中HNCS站FN由0變?yōu)?1.20 mm·a-0.25，其原因可能是該測(cè)站的噪聲并不符合WN+FN噪聲模型，而降噪會(huì)將非整數(shù)譜指數(shù)冪律噪聲PL變成FN。

圖8 降噪前后各測(cè)站坐標(biāo)時(shí)間序列的白噪聲和閃爍噪聲

圖9 白噪聲和閃爍噪聲改正率的測(cè)站分布

圖10為各測(cè)站垂向坐標(biāo)時(shí)間序列降噪前后的速度不確定度，其改正率的測(cè)站分布情況如圖11所示。由圖10可知，降噪后速度不確定度也明顯降低，平均改正率為78.68%，且Wilcoxon檢驗(yàn)表明，降噪前后序列的速度不確定度具有顯著性差異(p<0.05,n=289)。由圖11可知，35%的測(cè)站改正率在80%以上，63%的測(cè)站改正率在70%～80%之間。分析表明，LOESS方法可顯著降低序列的速度不確定度。

圖10 降噪前后各測(cè)站坐標(biāo)時(shí)間序列的速度不確定度

圖11 速度不確定度改正率的測(cè)站分布

4 結(jié) 語(yǔ)

本文將LOESS方法引入GNSS垂向坐標(biāo)時(shí)間序列的降噪研究中，并使用數(shù)學(xué)指標(biāo)、統(tǒng)計(jì)檢驗(yàn)和GNSS測(cè)站運(yùn)動(dòng)模型參數(shù)進(jìn)行效果分析，結(jié)果表明，LOESS方法可有效提高GNSS原始坐標(biāo)時(shí)間序列的精度。

分析表明，降噪后信噪比均值提高至8.14 dB，最高達(dá)48.47 dB；標(biāo)準(zhǔn)差平均減小1.38 mm，平均改正率為21.35%，最高達(dá)48.77%。對(duì)于降噪后的序列，白噪聲和閃爍噪聲均值分別由2.685 mm、12.415 mm·a-0.25降至0.006 mm、2.723 mm·a-0.25，平均改正率為78.63%；同時(shí)，速度不確定度得到明顯改善，平均改正率達(dá)78.68%。

本文在以下方面仍需進(jìn)一步研究：1)未將本文方法與其他降噪方法進(jìn)行對(duì)比，以分析本文方法的優(yōu)勢(shì)和不足；2)本文直接采用了普遍認(rèn)為的最佳噪聲模型WN+FN，未考慮GNSS基準(zhǔn)站復(fù)雜的噪聲特性。