李豐翔，林 琳

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

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

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

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

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

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

1.2 研究方法

1.2.1 PWV反演方法

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

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

利用Saastamoinen模型計(jì)算對(duì)流層干延遲[29]，即

式中，P為大氣壓；L為GNSS站點(diǎn)緯度；H為GNSS站點(diǎn)高程。

對(duì)流層總延遲減去干延遲得到對(duì)流層濕延遲，再根據(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時(shí)間序列預(yù)處理方法

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

1.2.3 利用SSA和FFT探測(cè)趨勢(shì)項(xiàng)和周期項(xiàng)的方法

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

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

2.1 PWV解算和精度評(píng)定

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

圖1 2010—2019年逐年平均PWV

表1 GAMIT解算策略

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

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

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

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

本文首先利用SSA方法對(duì)2010—2019年日本地區(qū)20個(gè)測(cè)站的PWV年平均值序列進(jìn)行分解，并利用FFT方法對(duì)分解所得的子序列進(jìn)行周期探測(cè)，提取序列中存在的趨勢(shì)項(xiàng)和周期項(xiàng)；然后利用SSA和最小二乘擬合相結(jié)合的方法擬合線性趨勢(shì)，該方法能避免周期項(xiàng)和噪聲對(duì)線性趨勢(shì)的影響，盡可能高精度地估計(jì)PWV年變化速率。分析結(jié)果見圖2，可以看出，日本地區(qū)PWV存在明顯的上升趨勢(shì)，上升速率約為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)明顯增大趨勢(shì)，最大達(dá)到1 mm，半年周期的振幅總體比周年周期小。

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

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

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

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

3 結(jié) 語(yǔ)

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