基于奇異譜分析的地基GPS 可降水汽含量變化規(guī)律分析

李豐翔，林 琳

（1.山東省國土測繪院，山東 濟(jì)南 250013；2.山東省水利科學(xué)研究院，山東 濟(jì)南 250014）

水汽是水循環(huán)中最活躍的成分，可降水量是天氣預(yù)報和氣候研究中的一個重要參數(shù)?？山邓浚≒WV）可反映大氣水汽變化，是大氣水分和熱量傳遞的基礎(chǔ)，是降水、蒸發(fā)和水汽輻合平衡的結(jié)果[1]。南亞季風(fēng)和西風(fēng)帶來的水汽在日本地區(qū)上空匯聚[2-5]。水汽在固體、液體、氣體3 種形態(tài)變化的過程中會產(chǎn)生能量的吸收或釋放，對維持大氣穩(wěn)定性和形成強(qiáng)對流天氣具有重要影響。水汽的移動和變化也對地區(qū)的氣候干旱具有重要影響[6-8]。此外，水汽還是臺風(fēng)等極端災(zāi)害性天氣形成和發(fā)展的重要因素[9]。因此，準(zhǔn)確獲取PWV 及其時空變化特征對于研究氣候變化、天氣預(yù)警等具有重要意義。

無線電探空儀、微波輻射法、太陽光度法和衛(wèi)星遙感[10]等PWV探測技術(shù)的空間分辨率、時間分辨率低且精度不高，其應(yīng)用受到了一定的限制[11-13]。GPS 具有分布廣泛、觀測精度高、全天候觀測、觀測成本低和觀測方便等優(yōu)點，已廣泛應(yīng)用于車輛和個人導(dǎo)航、工程和大地測量、地球物理和地球動力學(xué)研究、姿態(tài)確定等領(lǐng)域[14]。GPS 信號在通過對流層和電離層時會產(chǎn)生延遲，采用雙頻載波相位線性組合可估計或顯著降低電離層延遲的影響，通過雙差解算或精確點定位技術(shù)可實現(xiàn)對天頂總延遲的精確估計。Bevis M[15]等提出了GPS氣象學(xué)的概念，使PWV能在全天候條件下從GNSS 觀測數(shù)據(jù)中高精度、高時間分辨率檢索。一些研究[15-16]表明，GPS反演PWV的精度是可靠的，RMSE值為1～3 mm；且被廣泛應(yīng)用于極端短期災(zāi)害監(jiān)測和預(yù)報[17-18]等領(lǐng)域。Choy S[19]等基于GPS 導(dǎo)出的PWV 監(jiān)測了2010年發(fā)生在墨爾本的風(fēng)暴；Tang X[20]等研究發(fā)現(xiàn)，GNSS反演的PWV在2016年臺風(fēng)逼近浙江省時呈現(xiàn)持續(xù)增加趨勢。我國學(xué)者從20 世紀(jì)90 年代開始接觸學(xué)習(xí)GPS/MET原理方法，并開始將這一理論知識運用到實踐中，如將地基GPS反演技術(shù)應(yīng)用于對PWV的探測實驗中[21-22]。1999 年黨亞民[23]等首次選取位于北京、上海、武漢的3 個GPS 跟蹤站進(jìn)行實驗研究，得到了監(jiān)測站完整的PWV 時間序列，誤差結(jié)果小于1 mm；劉焱雄[24]等基于香港地區(qū)GPS 站的觀測數(shù)據(jù)，采用雙差網(wǎng)解方法計算測站的PWV 時間序列并進(jìn)行精度評定，結(jié)果表明GPS反演的PWV達(dá)到了較高的精度；鄒友峰[25]等通過計算GNSS 測站的PWV，分析了不同區(qū)域PWV的分布特點與規(guī)律；謝娜[26]等在深入分析GPS/PWV 理論的基礎(chǔ)上，結(jié)合多普勒天氣雷達(dá)在自然界觀測實驗中測得的液態(tài)水信號數(shù)據(jù)資料，初步研究得到GPS/PWV 與自然界平均降雨量的相互關(guān)系，為研究人工降雨提供了參考資料。因此，為了研究基于GPS技術(shù)反演PWV的時空分布特點與規(guī)律，本文以日本地區(qū)為例，基于GAMIT/GLOBK 10.7 軟件對數(shù)據(jù)進(jìn)行解算，得到高精度PWV 時間序列，并利用奇異譜分析（SSA）和快速傅里葉變換（FFT）方法研究其時空變化規(guī)律。

1 研究數(shù)據(jù)和方法

1.1 研究區(qū)概況與數(shù)據(jù)

本文選取日本作為研究區(qū)域。日本位于太平洋西北部海域，亞洲大陸東部地區(qū)，主體由九州、四國、本州、北海道四大島組成；屬溫帶海洋性季風(fēng)氣候，終年溫和濕潤，降雨較多；6 月多梅雨，夏、秋季多臺風(fēng)。本文選取分布在日本的20 個GPS 測站實測數(shù)據(jù)，數(shù)據(jù)的時間跨度為2010—2019年，數(shù)據(jù)采樣間隔為30 s。實驗數(shù)據(jù)來源于SONEL 網(wǎng)站（https://www.so?nel.org）。

1.2 研究方法

1.2.1 PWV反演方法

本文利用GAMIT/GLOBK 10.7 對GPS 數(shù)據(jù)進(jìn)行解算得到PWV。其原理為利用GNSS 載波相位觀測方程將對流層總延遲作為未知參數(shù)求解[27-28]，即

式中，i為接收機(jī)；j為衛(wèi)星；為接收機(jī)到衛(wèi)星的幾何距離；c 為光速；dti為接收機(jī)鐘差；dtj為衛(wèi)星鐘差；M為對流層映射函數(shù)；為 電離層延遲；λf為相位觀測波長；為相位模糊度；為其他誤差。

利用Saastamoinen模型計算對流層干延遲[29]，即

式中，P為大氣壓；L為GNSS站點緯度；H為GNSS站點高程。

對流層總延遲減去干延遲得到對流層濕延遲，再根據(jù)式（3）可求得GNSS反演PWV值[15]。

式中，gs為氣體常數(shù)（461±0.003 Jkg-1K-1）；ρw為液體水密度；k2、k3為折射常數(shù)，k2=24±11 Kpa-1，k3=3.75±0.03 K2Kpa-1；Tm為大氣垂直面的加權(quán)平均溫度。

1.2.2 PWV時間序列預(yù)處理方法

本文先對PWV 時間序列進(jìn)行粗差探測和剔除，再對剔除的數(shù)據(jù)和原本缺失的數(shù)據(jù)進(jìn)行插值處理。本文采用的粗差探測和剔除方法為3 倍標(biāo)準(zhǔn)差方法，首先對時間序列進(jìn)行扣除趨勢項和中心化處理，然后求取序列的標(biāo)準(zhǔn)差，將超過3 倍標(biāo)準(zhǔn)差的值認(rèn)定為粗差并剔除，統(tǒng)計分析發(fā)現(xiàn)，粗差剔除率為3.96%。本文采用SSA 插值方法。SSA[30]是一種對非線性、非平穩(wěn)時間序列進(jìn)行分析和預(yù)測的方法，可從具有未知先驗信息的時間序列中提取出變化緩慢的趨勢、周期、噪聲等有效信息，已被應(yīng)用于氣候?qū)W、測量學(xué)[31]和海洋學(xué)[32]。SSA 迭代插值[33-36]的主要步驟為：①采用符號“99999”對PWV 時間序列中的缺失數(shù)據(jù)進(jìn)行補全并進(jìn)行中心化處理；②利用SSA對補全后的PWV時間序列進(jìn)行分解；③采用加權(quán)相關(guān)分析方法確定前n階為序列的主成分，將第1個重構(gòu)成分中的值代替序列中缺失值，重復(fù)該步驟，直至前后兩次插補數(shù)據(jù)殘差的RMSE小于0.01 mm；④對第2個重構(gòu)成分重復(fù)步驟③，直至完成第n個重構(gòu)成分的插值。

1.2.3 利用SSA和FFT探測趨勢項和周期項的方法

PWV變化具有明顯的周期項和趨勢項，為探測日本地區(qū)PWV序列的周期項和趨勢項，本文引入SSA和FFT 進(jìn)行時間序列分析，具體步驟為：①輸入預(yù)處理后的PWV 時間序列；②基于PWV 時間序列構(gòu)建時滯矩陣；③奇異值分解，步驟②和③稱為SSA分解；④分組；⑤對角平均化，步驟④和⑤稱為SSA重構(gòu)；⑥利用FFT探測子序列的周期信號。

2 研究結(jié)果與分析

2.1 PWV解算和精度評定

本文采用高精度數(shù)據(jù)處理軟件GAMIT/GLOBK 10.7，在Linux系統(tǒng)下運行，具體解算策略見表1。逐年平均的PWV 解算結(jié)果見圖1，可以看出，2010 年起PWV 均在100 mm 以上，其中2010 年、2014 年、2015 年和2016 年的PWV 均高于平均水平，這可能與全球變化的氣候大環(huán)境有一定關(guān)系；PWV 含量呈現(xiàn)先減后增再減的趨勢，具有明顯的周期現(xiàn)象，2015年達(dá)到最大值，約為110.2 mm，2012 年達(dá)到最小值，約為100.6 mm。

圖1 2010—2019年逐年平均PWV

表1 GAMIT解算策略

本文采用RMSE 來評價PWV 的精度，對10 年的PWV 時間序列進(jìn)行精度統(tǒng)計，結(jié)果見表2，可以看出，日本地區(qū)20 個GPS 站10 年數(shù)據(jù)反演的PWV 的RMSE 的平均值、最大值、最小值和標(biāo)準(zhǔn)差分別為3.41 mm、4.42 mm、2.07 mm和0.58 mm，PWV的解算精度較高，滿足時空分布規(guī)律研究的精度需求。

表2 2010—2019年P(guān)WV精度統(tǒng)計結(jié)果/mm

式中，yi為第i年的PWV；?為2010—2019 年P(guān)WV均值；I為待分析的PWV總年數(shù)，本文取10。

2.2 PWV時空規(guī)律分析

本文首先利用SSA方法對2010—2019年日本地區(qū)20個測站的PWV年平均值序列進(jìn)行分解，并利用FFT方法對分解所得的子序列進(jìn)行周期探測，提取序列中存在的趨勢項和周期項；然后利用SSA和最小二乘擬合相結(jié)合的方法擬合線性趨勢，該方法能避免周期項和噪聲對線性趨勢的影響，盡可能高精度地估計PWV年變化速率。分析結(jié)果見圖2，可以看出，日本地區(qū)PWV存在明顯的上升趨勢，上升速率約為0.13 mm/a；日本地區(qū)PWV 存在明顯的周期現(xiàn)象，較突出的周期為周年周期（b1、b2）和半年周期（c1、c2），周年周期的振幅約為2 mm，振幅較大，半年周期的振幅存在明顯變化，2014 年之前振幅約為0.5 mm，2014—2016 年振幅約為0.8 mm，2017—2018 年振幅約為0.3mm，2018年以后振幅出現(xiàn)明顯增大趨勢，最大達(dá)到1 mm，半年周期的振幅總體比周年周期小。

圖2 2010—2019年日本地區(qū)PWV的SSA和FFT分析結(jié)果

同時，本文分析了2010—2019 年日本地區(qū)PWV均值的空間分布特點，結(jié)果見圖3，可以看出，日本地區(qū)PWV 變化隨緯度變化呈先減后增的趨勢，且PWV 的最低點出現(xiàn)在中部地區(qū)，北部PWV 總體最高；且沿海區(qū)域PWV比內(nèi)陸地區(qū)高。

圖3 2010—2019年日本地區(qū)PWV均值的空間分布

綜上所述，從時間維度分析，2010—2019 年日本地區(qū)PWV 變化具有明顯的增長趨勢，且伴有明顯的周期現(xiàn)象；從空間維度分析，日本地區(qū)沿海區(qū)域PWV 比內(nèi)陸高，且南北地區(qū)PWV 含量高、中部地區(qū)含量低。

3 結(jié) 語

本文利用GAMIT/GLOBK 10.7 軟件對日本地區(qū)2010—2019 年GPS 數(shù)據(jù)進(jìn)行解算得到高精度PWV 數(shù)據(jù)。首先利用3 倍標(biāo)準(zhǔn)差和SSA 方法對數(shù)據(jù)進(jìn)行預(yù)處理；然后利用SSA和FFT方法提取PWV序列的趨勢項和周期項，并利用SSA與最小二乘擬合相結(jié)合的方法估計PWV 變化速率，有效避免了周期項和噪聲對PWV 變化序列的影響，得到高精度的PWV 變化速率。本文的主要結(jié)論為：①利用GAMIT/GLOBK 軟件對GPS 數(shù)據(jù)進(jìn)行解算，解算精度為3.35 mm，精度較高；②通過對PWV 時間序列的粗差探測發(fā)現(xiàn)，2010—2019 年日本地區(qū)PWV 數(shù)據(jù)的粗差剔除率為3.98%，數(shù)據(jù)質(zhì)量較好；③日本地區(qū)PWV存在明顯上升趨勢，上升速率為0.13 mm/a，且存在明顯的周年周期和半年周期信號，周年周期信號強(qiáng)度高于半年周期，空間上內(nèi)陸地區(qū)PWV低于沿海地區(qū)。