韋 云 王 迅 王 浩 李 黎 陳國(guó)棟 趙 偉

1 蘇州科技大學(xué)地理科學(xué)與測(cè)繪工程學(xué)院，江蘇省蘇州市學(xué)府路99號(hào)，215009 2 蘇州科技大學(xué)北斗導(dǎo)航與環(huán)境感知研究中心，江蘇省蘇州市學(xué)府路99號(hào)，215009

近年來(lái)，利用由GNSS技術(shù)獲取的PWV時(shí)空變化趨勢(shì)來(lái)預(yù)測(cè)降雨的方法成為國(guó)內(nèi)外學(xué)者關(guān)注的重點(diǎn)[1-3]。傳統(tǒng)的獲取PWV的方法是利用由GNSS技術(shù)得到的對(duì)流層延遲(ZTD)減去干延遲(ZHD)，得到濕延遲(ZWD)，再利用ZWD推算得到PWV。計(jì)算ZHD需要的參數(shù)有測(cè)站處地面大氣壓(Ps)、測(cè)站緯度(θ)、測(cè)站大地高(H)，ZWD與PWV之間的轉(zhuǎn)換系數(shù)也要用到加權(quán)平均溫度(Tm)、地面氣溫(Ts)等參數(shù)，參數(shù)越多，計(jì)算過(guò)程越復(fù)雜，會(huì)導(dǎo)致數(shù)據(jù)量大、計(jì)算效率低，且易產(chǎn)生誤差累積等問(wèn)題。因此，建立ZTD與GNSS-PWV之間的直接轉(zhuǎn)換模型，以計(jì)算特定時(shí)段的PWV時(shí)序信息十分必要[4-6]。

本文基于2017年?yáng)|南沿海地區(qū)18個(gè)GNSS站的觀測(cè)數(shù)據(jù)，采用多元函數(shù)線性回歸法建立該地區(qū)GNSS-PWV與ZTD、地面氣溫(Ts)、地面大氣壓(Ps)之間的直接轉(zhuǎn)換模型，并利用模型預(yù)報(bào)2018年GNSS-PWV，以各站實(shí)測(cè)PWV為參考值，驗(yàn)證本文模型的預(yù)報(bào)精度。

1 數(shù)據(jù)來(lái)源及其相關(guān)性分析

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

選取由江蘇省氣象局和中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)(http:∥www.cgps.ac.cn)提供的2017～2018年?yáng)|南沿海地區(qū)18個(gè)GNSS站資料，主要包括PWV、ZTD、地面大氣壓(Ps)和地面氣溫(Ts)等，其中常州、黃山、蚌埠、武夷山4個(gè)氣象站用于精度檢驗(yàn)，其余14個(gè)測(cè)站用于建模，測(cè)站空間位置信息見(jiàn)表1，選取的18個(gè)GNSS站均勻分布于108°～122°E、19°～33°N區(qū)域，涵蓋了東南沿海大部分地區(qū)。

表1 GNSS測(cè)站空間位置信息

中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)采用Saastamoinen模型計(jì)算ZHD[7]，其表達(dá)式為：

(1)

用ZTD減去ZHD獲得ZWD，ZWD乘以轉(zhuǎn)換系數(shù)Π便可得到GNSS-PWV：

(2)

式中，ρw為液態(tài)水密度，k′2、k3為大氣折射常數(shù)，Rv為水汽氣體常數(shù)，Tm為加權(quán)平均溫度。

1.2 影響因子相關(guān)性驗(yàn)證

皮爾森(Pearson)相關(guān)系數(shù)R可用來(lái)描述2個(gè)變量之間的線性相關(guān)程度，R的絕對(duì)值越大，表明其相關(guān)性越強(qiáng)，在0.8～1.0之間表示極強(qiáng)相關(guān)，0.6～0.8之間為強(qiáng)相關(guān)，0.4～0.6之間為中等程度相關(guān)，0.2～0.4之間為弱相關(guān)，0.0～0.2則無(wú)相關(guān)性[8]。

PWV與ZTD、Ts、Ps之間的相關(guān)程度如圖1所示，可以看出，ZTD與PWV的相關(guān)系數(shù)為0.98，為極強(qiáng)相關(guān)，表明ZTD對(duì)PWV值有極大影響；PWV與Ts的相關(guān)系數(shù)為0.78，為強(qiáng)相關(guān)；PWV與Ps的相關(guān)系數(shù)為-0.65，為負(fù)強(qiáng)相關(guān)。

圖1 相關(guān)性分析

表2為各變量間相關(guān)性統(tǒng)計(jì)結(jié)果，可以看出，各自變量之間也存在一定的線性關(guān)系。若建模選取的樣本自變量之間本身就存在高度共線性關(guān)系，會(huì)導(dǎo)致線性回歸模型不穩(wěn)定，難以區(qū)分各自變量對(duì)模型結(jié)果的影響，因此需要對(duì)自變量進(jìn)行共線性分析。

表2 相關(guān)性分析統(tǒng)計(jì)

本文利用方差膨脹因子(variance inflation factor，VIF)分析各自變量之間的共線性，其表達(dá)式為：

VIFi=1/(1-R2)

(3)

式中，R為自變量x與其余自變量作回歸分析的復(fù)相關(guān)系數(shù)。通常認(rèn)為，當(dāng)0

2 PWV直接轉(zhuǎn)換模型的建立及精度分析

2.1 建模方法

利用多元函數(shù)線性擬合方法建立PWV直接轉(zhuǎn)換模型。設(shè)多元線性回歸模型為：

y=β0+β1x1+…+βmxm+ε，ε～N(0,σ2)

(4)

式中，y為因變量；x1、x2、…、xm為自變量；β0為常數(shù)項(xiàng)；β1、…、βm為回歸系數(shù)，是與自變量無(wú)關(guān)的未知參數(shù)；ε為誤差項(xiàng)[10]，是均值為0、方差σ2>0的不可觀測(cè)隨機(jī)變量。若有n個(gè)獨(dú)立觀測(cè)數(shù)據(jù)，則可列出n個(gè)關(guān)于未知參數(shù)x1、x2、…、xm的方程：

yi=β0+β1xi1+…+βmxim+εi

εi～N(0,σ2),i=1,…,n

(5)

則有：

PWV轉(zhuǎn)換模型可表示為：

Y=Xβ+ε，ε～N(0,σ2En)

(6)

式中，En為n階單位矩陣。該線性回歸模型可用最小二乘法原理求解，其中自變量的系數(shù)βi為最小二乘估計(jì)值。

2.1.1 多因子模型

基于PWV與ZTD、Ts、Ps的線性關(guān)系，PWV與其他參數(shù)函數(shù)的關(guān)系為：

PWV=β0+β1ZTD+β2Ts+β3Ps

(7)

將2017年14個(gè)測(cè)站的PWV與ZTD、Ts、Ps代入式(7)，利用最小二乘法可得到多因子PWV的直接轉(zhuǎn)換模型為：

PWV=-13.380 5+0.162 9ZTD+

0.074 1Ts-0.359 7Ps

(8)

2.1.2 雙因子模型

將2017年14個(gè)測(cè)站的PWV與ZTD、Ts代入式(7)，得到基于ZTD和Ts的雙因子PWV直接轉(zhuǎn)換模型：

PWV=-377.315 7+0.161 0ZTD+

0.415 7Ts

(9)

將PWV與ZTD、Ps代入式(7)，得到基于ZTD和Ps的雙因子PWV直接轉(zhuǎn)換模型：

PWV=3.335 4+0.165 8ZTD-0.3821Ps

(10)

2.1.3 單因子模型

假設(shè)PWV和ZTD之間的線性關(guān)系為：

PWV=β0+β1ZTD

(11)

將2017年14個(gè)測(cè)站的PWV和ZTD代入式(11)，得到基于ZTD的單因子PWV直接轉(zhuǎn)換模型：

PWV=-423.991 2+0.182 7ZTD

(12)

2.2 精度檢驗(yàn)

利用2017年14個(gè)測(cè)站的數(shù)據(jù)建立PWV模型，分別將2017～2018年的數(shù)據(jù)代入所建模型中，計(jì)算得到PWV模型預(yù)測(cè)值，并與原始數(shù)據(jù)中的PWV實(shí)測(cè)值(作為真值)進(jìn)行比較，驗(yàn)證模型的精度。

本文將PWV實(shí)測(cè)值減去PWV模型預(yù)測(cè)值，計(jì)算得到真值與預(yù)測(cè)值之間的偏差(bias)，通過(guò)統(tǒng)計(jì)分析平均bias與RMS值來(lái)驗(yàn)證PWV模型的精度。圖2為各PWV模型的bias變化趨勢(shì)，可以看出，單因子PWV模型預(yù)測(cè)值與真值的bias基本在15 mm以內(nèi)；雙因子(無(wú)Ps)PWV模型預(yù)測(cè)值與真值的bias在10 mm以內(nèi)；雙因子(無(wú)Ts)PWV模型預(yù)測(cè)值與真值的bias在2 mm以內(nèi)；多因子PWV模型預(yù)測(cè)值與真值的bias在1 mm以內(nèi)，精度最高。2018年的數(shù)據(jù)未參與建模，PWV模型預(yù)測(cè)值及bias的變化趨勢(shì)與2017年情況基本一致。

圖2 2017～2018年不同PWV模型預(yù)測(cè)值偏差

圖3為未參與建模的常州、黃山、蚌埠和武夷山等4個(gè)測(cè)站2017～2018年各PWV模型預(yù)測(cè)值bias。由圖可知，常州站和蚌埠站基于雙因子(無(wú)Ps)模型和單因子模型的預(yù)測(cè)值大部分大于真值，而基于雙因子(無(wú)Ts)模型的預(yù)測(cè)值大部分小于真值，基于多因子模型的預(yù)測(cè)值則均勻分布于真值兩側(cè)；黃山站和武夷山站基于雙因子(無(wú)Ps)模型和單因子模型的預(yù)測(cè)值均小于真值，基于雙因子(無(wú)Ts)模型的預(yù)測(cè)值均大于真值，而基于多因子模型的預(yù)測(cè)值同樣均勻分布于真值兩側(cè)。從bias的變化范圍來(lái)看，4個(gè)測(cè)站的PWV偏差值與參與建模的PWV偏差值分布基本一致，各類模型的bias變化范圍也與圖2非常接近。

圖3 2017～2018年各測(cè)站PWV模型預(yù)測(cè)值偏差

表3(單位mm)為各PWV模型預(yù)測(cè)值的平均bias和RMS統(tǒng)計(jì)結(jié)果，可以看出，基于單因子模型的RMS值為4.66 mm，平均bias為3.73 mm；加入Ts的雙因子模型RMS值降至3.94 mm，平均bias為3.13mm，精度有小幅提升；加入Ps的雙因子模型RMS值降至0.50 mm，平均bias為0.40 mm，精度大幅提高；同時(shí)加入Ps和Ts的多因子模型RMS值降至0.33 mm，平均bias降為0.24 mm，精度達(dá)到最高。

表3 全年P(guān)WV模型精度統(tǒng)計(jì)

由式(1)和式(2)可知，ZHD主要通過(guò)Ps計(jì)算得到，ZWD主要通過(guò)Ts及Tm計(jì)算得到，而Tm也由Ts計(jì)算得出。由于ZHD約占對(duì)流層延遲的80%～90%，加入Ps的雙因子模型能提升ZHD的計(jì)算精度，而ZWD在ZTD中的占比較小，因此基于ZTD與Ps的雙因子PWV模型的精度明顯比基于ZTD和Ts的雙因子PWV模型的精度高。

與基于ZTD、Ts和Ps的多因子PWV模型相比，基于ZTD和Ps的雙因子PWV模型的精度稍低，但其RMS值有0.50 mm，可滿足絕大部分GNSS氣象學(xué)研究的需求，且其計(jì)算參數(shù)較少，轉(zhuǎn)換效率較高。

3 結(jié) 語(yǔ)

本文利用2017～2018年?yáng)|南沿海地區(qū)18個(gè)GNSS站的數(shù)據(jù)，通過(guò)分析PWV與ZTD、地面氣溫(Ts)、地面大氣壓(Ps)等參數(shù)的相關(guān)性，基于多元線性擬合方法構(gòu)建東南沿海地區(qū)的多因子PWV直接轉(zhuǎn)換模型，得出以下結(jié)論：

1)PWV與ZTD之間的相關(guān)性最強(qiáng)，相關(guān)系數(shù)達(dá)到0.98；PWV與Ps、Ts之間也具有較強(qiáng)的相關(guān)性，相關(guān)系數(shù)分別為-0.65和0.78。

2)基于ZTD的單因子PWV模型平均bias為3.73 mm，RMS為4.66 mm；基于ZTD和Ts的雙因子PWV模型平均bias為3.13 mm，RMS為3.94 mm；基于ZTD和Ps的雙因子PWV模型平均bias為0.40 mm，RMS為0.50 mm；基于ZTD、Ts和Ps的多因子PWV模型平均bias為0.24 mm，RMS為0.33 mm。多因子PWV直接轉(zhuǎn)換模型的精度最高，基于ZTD和Ps的雙因子模型次之，這兩種模型均可獲得精度優(yōu)于1 mm的PWV結(jié)果。

3)東南沿海地區(qū)的PWV直接轉(zhuǎn)換模型不僅可應(yīng)對(duì)傳統(tǒng)PWV計(jì)算過(guò)程數(shù)據(jù)量大、效率低且易產(chǎn)生誤差累積的問(wèn)題，也可在基本地面氣象數(shù)據(jù)缺失時(shí)，基于ZTD直接計(jì)算特定時(shí)段的PWV時(shí)序信息，以用于GNSS氣象學(xué)研究。