鄒玉學(xué),詹金瑞,朱臘臘,上官明,歐陽光

(1.廣州市城市規(guī)劃勘測設(shè)計研究院，廣州 510060；2.中國地質(zhì)大學(xué)(武漢) 地理與信息工程學(xué)院，武漢 430074)

大氣可降水量(precipitable water vapor，PWV)的相位變化與降水直接相關(guān)，在天氣變化中起到了至關(guān)重要的作用[1-3]。全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system，GNSS)技術(shù)的快速發(fā)展為研究PWV提供全新的手段，克服水汽微波輻射計、無線電探空儀等方法的時空分辨率較低、成本高問題[4]。利用GNSS技術(shù)反演高精度PWV過程中，大氣加權(quán)平均溫度(weighted mean atmospheric temperature，Tm)是其關(guān)鍵參數(shù)之一。

國內(nèi)外學(xué)者對Tm模型研究主要集中在兩大部分：第一部分研究主要聚焦于基于無氣象要素建立的經(jīng)驗Tm模型，如Yao等[5]建立Tm與站點坐標(biāo)和年積日有關(guān)的GWMT模型，B?hm等[6]建立顧及年周期和半年周期的GPT2wTm模型，此模型顧及了季節(jié)性變化，僅根據(jù)測站經(jīng)緯度、高程和年積日等參數(shù)即可獲取Tm，因此精度較差，反演PWV時精度無法得到充分保障。第二部分研究主要集中在利用地表氣象數(shù)據(jù)建立的Tm模型，如Bevis等[1]利用美國地區(qū)探空數(shù)據(jù)建立Tm關(guān)于Ts(地表溫度)間的單因子線性回歸模型；多數(shù)學(xué)者研究發(fā)現(xiàn)Bevis模型應(yīng)用到其它地區(qū)存在明顯的系統(tǒng)性偏差，此后基于Bevis的線性回歸模型思想，諸多學(xué)者先后建立了適合極地、西非、歐洲和中國區(qū)域的本地線性回歸Tm模型[7-10]；姚宜斌等[11]推導(dǎo)Tm與Ts間的非線性關(guān)系，突破傳統(tǒng)的Tm與Ts間為一元線性回歸關(guān)系的認(rèn)識，建立Tm與Ts間的單因子非線性模型，并認(rèn)為擬合曲線在低溫地區(qū)擬合精度較低。Zhang等[12]研究線性Tm-Ts模型殘差的半年周期性時發(fā)現(xiàn)其隨緯度增加逐漸增大，因此在線性模型中加入諧波函數(shù)，合理解釋線性模型殘差周期性變化情況；臧建飛等[13]認(rèn)為多因子線性回歸Tm模型并不為最優(yōu)，基于誤差周期性分析，考慮年周期和半年周期性誤差修正現(xiàn)有的線性模型、指數(shù)模型和混合模型，而新模型在中高緯地區(qū)依然存在較大誤差。由此可見，現(xiàn)有Tm模型在我國東北地區(qū)普遍精度較低。為提高GNSS反演PWV精度，建立一個適合我國東北地區(qū)的Tm模型是當(dāng)前研究的首要任務(wù)。隨著計算機(jī)科學(xué)的發(fā)展，計算機(jī)的計算能力獲得急劇提升，許多數(shù)據(jù)驅(qū)動型模型得以實現(xiàn)并廣泛應(yīng)用于各行各業(yè)的生產(chǎn)研究工作中。Ding[14]以較少的地表氣象要素做模型輸入，構(gòu)建基于神經(jīng)網(wǎng)絡(luò)的Tm模型，結(jié)果表明該模型在全球尺度上優(yōu)于GPT2w、GTm、GTm-I和PTm-I這4類舊模型；龍鳳陽等[15]利用神經(jīng)網(wǎng)絡(luò)強(qiáng)大的非線性優(yōu)勢，基于補(bǔ)償方法構(gòu)建了適合中國區(qū)域的Tm融合模型。已有學(xué)者基于TFNN分別建立了NMFTm、SMFTm和MFMT模型，基于集成學(xué)習(xí)方式增強(qiáng)了模型的泛化能力[16]。多層感知器模型(multilayer perceptron,MLP)有著靈活的輸入和輸出且能解決非線性映射問題，具有很強(qiáng)的學(xué)習(xí)記憶功能。MLP模型的初始權(quán)值和閾值是隨機(jī)性的,難以獲得全局性最優(yōu)初始值，計算結(jié)果不穩(wěn)定，GA算法作為全局搜索算法恰能夠彌補(bǔ)誤差反傳過程中的缺陷，具有更好的非線性擬合能力[17]。

利用GA-MLP對Tm模型的建立卻鮮有報道。文中采用美國懷俄明大學(xué)(university of wyoming，UW)提供的東北地區(qū)7個國際探空站氣象數(shù)據(jù)(數(shù)據(jù)來自于http://weather.uwyo.edu/upperair/sounding.html)，通過GA-MLP模型，融合氣象數(shù)據(jù)，建立適合東北地區(qū)的Tm模型，驗證該模型在該區(qū)域的適用性。

1 計算原理

1.1 常用的Tm的計算方法

在GNSS氣象學(xué)中，PWV與對流層天頂濕延遲(zenith wet delay，ZWD)之間的轉(zhuǎn)換系數(shù)Π可由式(1)得到[1]：

(1)

(2)

式中：e，T分別為地表上的水氣壓(hPa)和溫度(K)；dz為積分高度；ei，Ti分別為第i層大氣層頂?shù)乃畾鈮汉徒^對溫度；Δhi為第i層的層高。

由于探空站在全球的分布數(shù)量有限，且僅提供每天的0時和12時的觀測數(shù)據(jù)，基于式(2)計算的Tm空間分辨率和時間分辨率較低，常用的Tm計算方式有以下幾種：

1)常數(shù)法：王小亞等發(fā)現(xiàn)轉(zhuǎn)換系數(shù)Π是隨季節(jié)和氣候變化的，近似值為0.15[18]，由式(1)可反求Tm=264.42 K，此計算的Tm無法反映其時空變化；

2)GPT2w模型：利用GPT2w模型結(jié)合參數(shù)格網(wǎng)文件可以求解全球任一點的Tm經(jīng)驗值[6]。

3)單/多因子線性回歸模型：多數(shù)學(xué)者根據(jù)Bevis模型(Tm=a+bTs)建立適合本地的單因子線性回歸模型，也有部分學(xué)者建立了基于氣象要素的多因子線性回歸模型，此類回歸分析方式計算的Tm值最為常用；

1.2 MLP模型構(gòu)建

MLP神經(jīng)網(wǎng)絡(luò)屬于人工智能常用的算法，對自變量與因變量間的擬合效果幾乎可以逼近任意非線性函數(shù)。它是一種多層前饋網(wǎng)絡(luò)，按照誤差反向傳播算法進(jìn)行訓(xùn)練，通過逆向傳播不間斷地修改網(wǎng)絡(luò)的閾值與權(quán)值，從而使網(wǎng)絡(luò)誤差達(dá)到最小[19]。

MLP由輸入層、隱含層和輸出層三部分構(gòu)成，常用的三層網(wǎng)絡(luò)模型結(jié)構(gòu)如圖1 所示，x1，x2，…，xn為輸入數(shù)據(jù)；y1，y2，…，y0為輸出數(shù)據(jù)；ωij為第i個輸入到第i個隱含層節(jié)點之間的連接權(quán)重；ωjk為第j個隱含層節(jié)點到第k個輸出節(jié)點的連接權(quán)重。輸入信號x經(jīng)過隱層被激活一次，后達(dá)到輸出層再次經(jīng)過激活函數(shù)得到輸出。

(3)

式中：zj為第j個隱節(jié)點的輸出；f為激活函數(shù)(Sigmoid、Tanh或者ReLU等)；ωij為第i個輸入到第j個隱節(jié)點之間的連接權(quán)重；xi為第i個輸入；bj為第j個隱節(jié)點的閾值；yk為輸出層的第k個輸出；ωjk為第j個隱節(jié)點到第k個輸出節(jié)點的連接權(quán)重；bk為第k個輸出節(jié)點的閾值。

圖1 MLP模型結(jié)構(gòu)

根據(jù)Tm的定義，氣壓與Tm無關(guān)，且二者間的Pearson相關(guān)系數(shù)較低，因此文中選取的輸入?yún)?shù)為與TmPearson相關(guān)性最高的Ts和Td(露點溫度)，輸出參數(shù)為單個輸出的Tm值，采用經(jīng)驗試算法最終確定隱含層神經(jīng)元個數(shù)為9個。隱含層和輸出層函數(shù)的確定對MLP模型預(yù)測精度影響較大，文中采用雙曲正切函數(shù)(Tanh)為隱含層節(jié)點傳遞函數(shù)，采用線性傳遞函數(shù)為輸出層節(jié)點傳遞函數(shù)，通過梯度下降算法對網(wǎng)絡(luò)進(jìn)行訓(xùn)練。

GA算法是模仿生物進(jìn)化理論機(jī)制的一種全局搜索方法，恰好能夠彌補(bǔ)誤差反傳過程中的缺陷[20]。GA算法使用生物進(jìn)化的思想，設(shè)定一定數(shù)量的初始個體，計算初始適應(yīng)度，選擇適應(yīng)度較大的個體進(jìn)行交叉、變異，并淘汰適應(yīng)度較小的個體，從而使整個群體產(chǎn)生進(jìn)化。在多層感知器中結(jié)合遺傳算法，具體過程如下：

1)MLP輸入：讀取數(shù)據(jù)并歸一化，確定MLP模型的拓?fù)浣Y(jié)構(gòu)；

2)編碼：將多層感知器的權(quán)值閾值視為基因，賦予多組不同的初始值，生成多個個體；

3)計算適應(yīng)度：將個體代入多層感知器計算誤差，以誤差的倒數(shù)作為個體適應(yīng)度；

4)選擇：根據(jù)適應(yīng)度的大小，采用輪盤賭的方式選擇個體；

5)交叉：對選擇的兩個個體的某一隨機(jī)位置進(jìn)行交叉操作；

6)變異：避免GA算法在訓(xùn)練中陷入局部最優(yōu)解，需將值在限制范圍內(nèi)進(jìn)行調(diào)整；

7)進(jìn)化：在個體基因經(jīng)過選擇、交叉和變異之后，新一代的個體產(chǎn)生，將新一代個體放入多層感知模型再次計算適應(yīng)度；

8)迭代：重復(fù)過程4)～7)直到最大迭代次數(shù)后輸出。

在經(jīng)過遺傳算法后，將最優(yōu)個體的基因作為初始權(quán)值和閾值帶入MLP模型參與訓(xùn)練，預(yù)測數(shù)據(jù)反歸一化并得到最終預(yù)測值。GA-MLP模型的流程圖如圖2所示。

圖2 算法流程圖

2 模型評估與分析

2.1 數(shù)據(jù)來源與精度評定指標(biāo)

圖3 探空站點分布

為評價模型的表現(xiàn)，文中以均方根誤差(root mean square error,RMSE)、平均偏差(Bias)和判定系數(shù)R2做模型精度評估標(biāo)準(zhǔn)。3種評估標(biāo)準(zhǔn)指標(biāo)如下：

(4)

(5)

(6)

2.2 結(jié)果分析

2.2.1Tm模型精度分析

為了比較文中所建基于GA-MLP的Tm模型和現(xiàn)有的Tm模型在東北地區(qū)與Tm真實值的擬合情況，分別計算了2018年7個探空站的Tm值，并按月取均值，模型化Tm值與真實值對比結(jié)果如圖4所示。

由圖4可以看出，各類模型化Tm值與真實值變化趨勢基本一致。無氣象要素參與的GPT2w模型值低于真實值，在7個探空站上具有明顯的系統(tǒng)性偏差，模型精度最差。具有氣象要素參與的各類模型化Tm值與真實值擬合度較高，但單因子線性和非線性模型在個別站點的部分時間段與真實值間的波動較為明顯，如嫩江、哈爾濱和通遼站的單因子線性和非線性模型值在4月和5月明顯高于真實值，嫩江、伊春和臨江站的單因子線性和非線性模型值在7月和8月明顯低于真實值，這可能由于僅有單一Td參與的模型無法更加準(zhǔn)確地反映Tm值。具有Ts和Td參與的多因子線性模型和基于GA-MLP的Tm模型與Tm真實值的擬合情況全年無明顯變化，Tm值變化趨勢一致，模型擬合情況最好。

為了進(jìn)一步定量分析不同Tm模型在2018年東北地區(qū)探空站的精度檢驗結(jié)果，計算了各類模型在2018年的Bias、RMSE和R2，結(jié)果如表1所示。

表1 不同模型年均精度檢驗結(jié)果

圖4 模型化Tm值與真實值對比

由表1可以看出，無氣象要素的GPT2w模型RMSE為 6.37 K，Bias為2.03 K，在東北地區(qū)精度最差，具有明顯的系統(tǒng)偏差，R2接近0.8，模型擬合優(yōu)度較差；具有氣象要素參與的各類模型的RMSE均低于GPT2w模型，Bias接近0 K，模型均無明顯系統(tǒng)偏差，判定系數(shù)均接近0.9，模型擬合優(yōu)度都較高，各類精度檢驗指標(biāo)明顯優(yōu)于無氣象要素參與的GPT2w模型，這說明氣象要素在Tm模型中具有至關(guān)重要的作用；僅有Ts參與的單因子線性模型和非線性模型的各項精度指標(biāo)無明顯差別，RMSE均接近5.0 K，R2=0.878，單因子線性模型和非線性模型在東北地區(qū)的年均精度相當(dāng)；多因子線性模型的各項精度指標(biāo)明顯優(yōu)于單因子線性模型，RMSE為4.13 K，R2=0.917，由此說明多因子較單因子更能反映模型預(yù)報的真實情況；GA-MLP的Tm模型的各項精度指標(biāo)均為最優(yōu)，較常用的單因子、多因子線性模型和MLP模型，RMSE分別減少了20%、3%和8%；對比MLP模型，采用GA優(yōu)化后的MLP模型精度提高明顯，平均偏差更接近0 K，判定系數(shù)更接近1，擬合優(yōu)度更優(yōu)。由此說明基于GA優(yōu)化MLP的Tm模型在東北地區(qū)精度最高，具有更優(yōu)的適用性。

2.2.2 GPS反演PWV精度分析

建立精確的區(qū)域大氣加權(quán)平均溫度模型最終目的是為了提高GNSS反演PWV的精度。利用東北地區(qū)IGS長春站2018年間的GPS觀測數(shù)據(jù)，由GAMIT軟件進(jìn)行高精度的基線解算得出對流層天頂總延遲(zenith total delay，ZTD)，從中減去由Saastamoinen模型[1]獲取的對流層天頂靜力學(xué)延遲(zenith hydrostatic delay，ZHD)，進(jìn)而得到與PWV直接相關(guān)的對流層天頂濕延遲(Zenith Wet Delay，ZWD)，ZWD需乘以轉(zhuǎn)換系數(shù)Π才能得到PWV[21]。分別通過最為常用且較易獲取Ts的單因子線性Tm模型、GPT2w模型、GA-MLP的Tm模型和真值Tm得到相應(yīng)的GPS PWV，以臨近的長春探空PWV為真值進(jìn)行對比分析(圖5)，分析結(jié)果見表2。

圖5 2018年長春PWV

表2 PWV精度比較

由圖5可以看出，3類模型反演的PWV變化趨勢一致。基于各類Tm模型反演的GPS PWV與探空PWV在冬季和春季基本吻合，兩類PWV數(shù)據(jù)一致性春冬明顯好于夏秋兩季，由于IGS站和探空站并非并置，可能因夏秋兩季局部天氣存在多變。由表2可知，基于GPT2w模型得到的GPS PWV Bias為1.36 mm，RMSE為3.68 mm；基于單因子線性模型得到的GPS PWV Bias為1.21 mm，RMSE為3.54 mm；基于本地化模型GA優(yōu)化的MLP模型得到的GPS PWV Bias為1.18 mm，RMSE為3.50 mm，符合國際上要求的3～4 mm的精度[21]。較普遍采用的GPT2w模型和單因子線性模型，GA算法優(yōu)化的MLP神經(jīng)網(wǎng)絡(luò)模型的RMSE降低了4.9%和1.1%，Bias降低了13.2%和2.5%，但三類模型Bias均接近1 mm，可能由于探空站與IGS站并非并置，存在一定系統(tǒng)偏差。由于GPT2w模型在東北地區(qū)本身存在系統(tǒng)性偏差，所以其反演的GPS PWV偏差稍大。簡而言之，基于GA優(yōu)化MLP的Tm模型在東北地區(qū)具有更優(yōu)的適用性，可以提高東北地區(qū)采用地基GPS技術(shù)反演大氣可降水量的精度。

3 結(jié)束語

文中采用東北地區(qū)7個國際探空站2014—2017年間的氣象數(shù)據(jù)，利用GA優(yōu)化的MLP模型，建立了適合東北地區(qū)的Tm模型。對2018年的Tm進(jìn)行預(yù)測，并驗證其應(yīng)用于GPS 反演PWV時精度，經(jīng)文中分析后得出如下結(jié)論：

1)GA-MLP的Tm模型在2018年東北地區(qū)的預(yù)測Bias為0.04 K，RMSE為4.06 K，R2為0.920，模型無系統(tǒng)性偏差，預(yù)測精度較高。較常用的Tm模型和GA-MLP的Tm模型各項精度評估指標(biāo)均為最優(yōu)。

2)文中建立的GA-MLP的Tm模型在GPS反演PWV的應(yīng)用中較常用的單因子線性模型和無氣象要素參與的GPT2w模型具有更高的精度，RMSE精度分別提升1.1%和4.9%，Bias精度分別提升2.5%和13.2%。

采用GA優(yōu)化MLP的Tm模型進(jìn)行GPS反演PWV是可行的。由于探空站與IGS站并非并置，以探空站PWV為真值檢驗，可能存在系統(tǒng)偏差影響，在后續(xù)實驗中有待進(jìn)一步驗證。