基于中尺度氣象模式（MM5）計算GNSS對流層水汽延時①

宋愛虎，高興國，孔 梅，江峻毅，周 寧

（1.山東電力工程咨詢院有限公司，山東 濟南250013；2.山東省土地勘測規(guī)劃院，山東 濟南250014）

0 引 言

目前，基于GNSS衛(wèi)星信號反演水汽信息遙感方法形成的GNSS氣象學(xué)（GNSS MET）已經(jīng)得到了廣泛關(guān)注。我國香港、上海、北京等地區(qū)都陸續(xù)建立了GNSS水汽監(jiān)測網(wǎng)，國家海洋局也正在我國沿海地區(qū)建設(shè)GNSS業(yè)務(wù)化觀測系統(tǒng)。隨著國內(nèi)外地基GNSS水汽監(jiān)測系統(tǒng)的建立，GNSS水汽遙感技術(shù)逐步走向成熟，并實現(xiàn)業(yè)務(wù)化運行。

國內(nèi)外多種試驗表明：利用GNSS遙感測站上空可降水量（PWV）的精度可達1～2mm［1－2].采用GNSS探測水汽的空間三維分布，必須得到衛(wèi)星信號傳播方向的斜向路徑延時（STD），進而轉(zhuǎn)化為傾斜路徑的水汽總量（SWV）。求定衛(wèi)星信號斜向路徑方向上的SWV，以及層析大氣水汽的垂直結(jié)構(gòu)是當(dāng)前GNSS氣象學(xué)領(lǐng)域的研究前沿。國內(nèi)外學(xué)者在斜向路徑延時方面的研究已經(jīng)取得了一些進展。Ware等［3]人采用GNSS雙差解算的方法證實了獲取STD的可行性；2000年Alber等［4]提出了一種從雙差中獲得信號路徑相位延遲的方法；宋淑麗等［5]則基于非差算法提出一種以毫米級精度直接計算SWV的方法等。為獲得更加精確的斜向水汽延時，以評價采用各種方法得出的信號方向的斜向延時結(jié)果，采用世界上先進的中尺度氣象模式MM5模擬得到的三維氣象場數(shù)據(jù)，建立起渤海灣區(qū)域的空間折射率場，進而計算得到任意衛(wèi)星高度角方向的水汽延時。同時也為層析空間水汽和其他領(lǐng)域應(yīng)用研究提供較高精度的斜向水汽比對結(jié)果。

1 GNSS氣象學(xué)原理與斜向水汽延時的定義

由GNSS衛(wèi)星傳播到地面接收設(shè)備的衛(wèi)星信號，在穿越大氣層的電離層和對流層過程中會發(fā)生延時。在GNSS數(shù)據(jù)處理過程中，通常采用兩個不同的頻率，利用不同觀測值組合可以消除電離層延時。而對流層延時中包含大氣中溫度、氣壓和相對濕度等有用的氣象信息。GNSS信號在對流層中的延時包含兩部分：一部分主要是由水汽引起，稱為濕分量延時；另一部分由氮、氧等氣體引起，由于氮、氧等氣體遵守靜水力學(xué)規(guī)則，這部分延時稱為靜水力學(xué)延時（俗稱干分量延時），并可以用數(shù)學(xué)物理模型計算。如：Saastamoinen模型、Hopfield模型和Black模型［6]。因此，GNSS數(shù)據(jù)處理得到的對流層延時減去模型計算的靜水力學(xué)延時后，可以得到濕分量延時。濕分量延時與氣象學(xué)中的可降水份量之間具有對應(yīng)關(guān)系，其投影因子為

式中：μ為水的密度；Rv為水汽常量；k1、k2和k3為大氣折射因子；w為水汽分子與干空氣分子的質(zhì)量比；Tm為對流層加權(quán)平均溫度。通過物理變換，可以將濕分量延時轉(zhuǎn)化成可降水量。這就是GNSS遙感水汽的基本原理。

GNSS信號傳播方向的斜向水汽延時定義為測站與衛(wèi)星間信號傳播方向上的水汽折射率的積分。計算斜向水汽濕分量延時的理論表達形式如下［7]

式中：Nw為濕氣體折射率；Rd為干氣體常數(shù)；q為相對濕度；ρ為大氣的密度；k1＝0.00776（K／Pa）；k2＝0.704（K／Pa）；k2＝3740（K2／Pa）；ε為干氣體常數(shù)與濕氣體常數(shù)的比值。

2 空間折射率場的建立

2.1 渤海灣MM5模式三維氣象場數(shù)據(jù)的生成

試驗區(qū)域為渤海灣海域，數(shù)據(jù)的產(chǎn)生采用目前世界上先進的中尺度氣象模式MM5。

MM5模式是美國賓洲大學(xué)國家大氣研究中心（PSU／NCAR）在原有的流體靜力模式MM4基礎(chǔ)上發(fā)展的新一代中尺度非流體靜力模式，當(dāng)前采用版本為V3.7版本，它是目前國際上公認的、先進的、應(yīng)用成熟的有限區(qū)域模式。相比WRF等其他模式，MM5的特點是適合中尺度的模擬，計算效率高，適合長時間中尺度的模擬統(tǒng)計分析。模式的網(wǎng)格設(shè)置為三重嵌套網(wǎng)格，粗細三個網(wǎng)格的水平分辨率分別為27km、9km、3km.垂直方向上共22層，D01區(qū)域時間步長是60s，D03區(qū)域時間步長為6.667秒，每次積分1440min（一天）。物理過程參數(shù)化使用了無積云對流參數(shù)化方案、高分辨率行星邊界層（MRF）方案、簡單冰顯式水汽方案、云輻射方案、5層土壤模式。

我們采用NCEP／NCAR的FNL水平分辨率為1.0°×1.0°的全球分析資料，為同化模擬提供較高分辨率的背景場和側(cè)邊界條件，同時加入GTS常規(guī)觀測和非常規(guī)資料，利用MM5模式的4DDA格點nudging同化，對渤海灣進行了高分辨率的數(shù)值同化模擬，得到試驗所需要的連續(xù)逐小時22層

式中：h為任意高度；N（h）為相應(yīng)高度的折射率；Ti（h）為切比雪夫多項式的遞推公式；ai切比雪夫多項式的系數(shù)；hi為粗略折射率場中第i層的高度。

第二步：將變量位高h根據(jù)公式轉(zhuǎn)換為變量τ，并根據(jù)移動窗口的個數(shù)確定多項式次數(shù)n的值（n＝3）氣象場。在渤海灣區(qū)域利用雙線性插值程序，將3 km網(wǎng)格上的溫度、水汽、高度場插值到1／60°的經(jīng)緯度網(wǎng)格上。

2.2 折射率格網(wǎng)的切比雪夫多項式擬合與內(nèi)插

根據(jù)MM5模擬得到的三維氣象數(shù)據(jù)計算出氣象場網(wǎng)格點的折射率，建立空間折射率場。對求得的折射率進行積分，由于產(chǎn)生的折射率網(wǎng)格太粗，求得的對流層濕延時難以滿足精度要求，因此需要在任一衛(wèi)星高度角方向，根據(jù)網(wǎng)格點的折射率采用移動窗口方法擬合出多個多項式，根據(jù)多項式對所需折射率進行內(nèi)插，以提高內(nèi)插精度。在多種擬合方法中，以網(wǎng)格點折射率為節(jié)點的切比雪夫多項式擬合內(nèi)插效果最佳［8]。

實現(xiàn)以網(wǎng)格點折射率為節(jié)點的切比雪夫多項式擬合，需要建立折射率關(guān)于位勢高度的切比雪夫多項式，首先，確定切比雪夫函數(shù)關(guān)系式

第三步：切比雪夫多項式中，根據(jù)遞推公式確定T

第四步：對每一個窗口的折射率內(nèi)節(jié)點建立誤差方程

式中：vi代表折射率誤差；ai代表系數(shù)；右側(cè)代表誤差方程的矩陣形式

最后，根據(jù)最小二乘法則計算每一個移動窗口內(nèi)多項式的系數(shù)陣M

根據(jù)計算出的每個窗口內(nèi)的切比雪夫多項式，以100m的高度間隔內(nèi)插得到任意高度角方向上的折射率。建立起更高分辨率的空間折射率場。

3 斜向水汽延時的確定

區(qū)域性折射率場的高度，并不能覆蓋接收機到GNSS衛(wèi)星方向的大氣高度。MM5模式得到的折射率場最大高度為17km，而水汽對GNSS信號的延時影響主要集中在一萬米以下的大氣層。因此，計算GNSS信號方向的斜向水汽延時需要分為兩部分，折射率場延時和高空超出折射率場范圍的延時。

對于折射率場范圍內(nèi)的水汽延時，在內(nèi)插得到的濕折射率的基礎(chǔ)上采用式（6）計算任意衛(wèi)星高度角方向的水汽延時

式中：e為衛(wèi)星高度角；swdf（e）表示e高度角方向上氣象場范圍內(nèi)的水汽延時；Nwi為第i層由氣象場數(shù)據(jù)計算得到的濕折射率。

對于超出折射率場范圍的水汽延時，考慮到這一部分對信號延時的影響很小，基于氣象場頂層氣象數(shù)據(jù)采用了Saastamoinen模型（式7）計算得到天頂方向的水汽延時［9－11]

式中：φ為氣象場頂層與信號方向交點緯度（以弧度為單位）；H 是交點的正高（km）；P為頂層交點大氣壓（以mbar為單位）。

然后根據(jù)Niell映射函數(shù)將天頂濕延時映射到對應(yīng)高度角方向，從而得到最終斜向水汽延時Swd.

式中，Niell映射函數(shù)系數(shù)a、b、c與測站的位置、年積日和高度角e有關(guān)的量，可以通過相應(yīng)的經(jīng)驗公式和不同緯度的經(jīng)驗值內(nèi)插得到。

4 確定水汽延時的精度分析

4.1 空間折射率內(nèi)插精度分析

空間濕折射率的擬合與內(nèi)插精度決定著斜向水汽延時的精度，通過移動窗口的切比雪夫擬合與內(nèi)插得到更高分辨率的空間濕折射率場。表1給出了內(nèi)插值與MM5數(shù)據(jù)計算得到的濕折射率值的比較結(jié)果。從比較結(jié)果來看：在不同高度角下內(nèi)插得到的濕折射率與MM5數(shù)據(jù)計算得到的濕折射率最大差值為0.3146，最小差值近乎為零。差值的最大標準中誤差為0.2591，最小標準中誤差為0.0651，取得了較高的精度。因此，采用移動窗口的切比雪夫多項式方法對空間折射率進行擬合內(nèi)插是可行的。

表1 不同高度角方向濕折射率內(nèi)插值與MM5計算值比較統(tǒng)計結(jié)果

4.2 水汽延時的外符合精度檢核

為了評價斜向水汽延時的計算精度，對氣象場數(shù)據(jù)計算天頂水汽總延時和濕延時分別與GPS結(jié)果進行了比較，GPS數(shù)據(jù)處理采用了瑞士伯爾尼大學(xué)的Bernese5.0軟件。由圖1和圖2給出了環(huán)渤海灣區(qū)域的測站HB09、HC03的比較結(jié)果。GPS天頂總延時與氣象模式結(jié)果仍存在差異，二者之間最大差異2cm，差異絕大多數(shù)都小于1cm；對于HC03濕延時的比較結(jié)果，發(fā)現(xiàn)與測站HB09趨勢一致，但差異明顯減小，最大為10mm，基本都控制在5mm左右，由此可以得出：基于MM5模式成果積分計算的信號傳播方向的水汽延時結(jié)果是可信的，具有比較理想的精度。

5 結(jié) 論

切比雪夫多項式逼近具有收斂快、擬合程度高的優(yōu)點，采用切比雪夫的移動窗口擬合法能夠提高空間折射率的內(nèi)插精度，得到了更高分辨率的空間濕折射率場。結(jié)果表明：MM5模式數(shù)據(jù)計算的斜向水汽延時結(jié)果可信，精度較高，與GPS計算結(jié)果差異不大而且具有較強的一致性，由MM5模式計算的斜向水汽延時除以天頂方向的結(jié)果得到的映射函數(shù)值較為可靠，完全可以采用MM5模式計算的映射值作為評價當(dāng)前高精度GNSS數(shù)據(jù)處理軟件中各個映射函數(shù)模型的標準。

［1]Bevis M S，Chiswell B S.Mapping zenith wet delays onto perceptible water［J].J.Appl.Meteoro.，1994（33）：379－386.

［2]張雙成，劉經(jīng)南，葉世榕，等.顧及雙差殘差反演GPS信號方向的斜路徑水汽含量［J].武漢大學(xué)學(xué)報·信息科學(xué)版，2009，34（1）：100－104.

［3]Ware R，Alber C，Rocken C，et al.Sensing integrated water vapor along GPS ray paths［J].Geophys.Res.Lett.，1997，24（4）：417－420.

［4]Alber C，Ware R，Rocken C，et al.Obtaining single path phase delays from GPS double differences［J].Geophys.Res.Lett.，2000，27（17）：661－664.

［5]宋淑麗.地基GPS.網(wǎng)對水汽三維分布的監(jiān)測及其在氣象學(xué)中的應(yīng)用［D].上海：中國科學(xué)院上海天文臺，2004.

［6]劉焱雄，陳永奇.GPS氣象學(xué)中垂直干分量延時的精確確定［J].測繪學(xué)報，2000，29（2）：172－180.

［7]De Haan S，Van der Marel H.The influence on GPS estimates of NWP－derived mapping functions［J].Physicas and Chemistry of the Earth，2004，29（2／3）：159－166.

［8]奚梅成，劉儒勛.數(shù)值分析方法［M].合肥：中國科學(xué)技術(shù)大學(xué)出版社，1995.

［9]Liu Y X.Remote sensing of atmospheric water vapor using GPS data in the hong Kong region［D].Hong Kong Polytechnic University（People’s Republic Of China），2000年.

［10]李國平.地基GPS遙感大氣可降水量及其在氣象中的應(yīng)用研究［D].四川：西南交通大學(xué)，2007年.

［11]Boehm J，Werl B，Schuh H.Troposhere mapping function for GPS and very long baseline interferometry from European Centre for Medium－Range Weather Forecasts operational analysis data［J].Journal of Geophysical Research，2006，111（2）：406－415.