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

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

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

0 引 言

目前，基于GNSS衛星信號反演水汽信息遙感方法形成的GNSS氣象學（GNSS MET）已經得到了廣泛關注。我國香港、上海、北京等地區都陸續建立了GNSS水汽監測網，國家海洋局也正在我國沿海地區建設GNSS業務化觀測系統。隨著國內外地基GNSS水汽監測系統的建立，GNSS水汽遙感技術逐步走向成熟，并實現業務化運行。

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

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

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

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

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

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

2 空間折射率場的建立

2.1 渤海灣MM5模式三維氣象場數據的生成

試驗區域為渤海灣海域，數據的產生采用目前世界上先進的中尺度氣象模式MM5。

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

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

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

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

2.2 折射率格網的切比雪夫多項式擬合與內插

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

實現以網格點折射率為節點的切比雪夫多項式擬合，需要建立折射率關于位勢高度的切比雪夫多項式，首先，確定切比雪夫函數關系式

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

第四步：對每一個窗口的折射率內節點建立誤差方程

式中：vi代表折射率誤差；ai代表系數；右側代表誤差方程的矩陣形式

最后，根據最小二乘法則計算每一個移動窗口內多項式的系數陣M

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

3 斜向水汽延時的確定

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

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

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

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

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

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

式中，Niell映射函數系數a、b、c與測站的位置、年積日和高度角e有關的量，可以通過相應的經驗公式和不同緯度的經驗值內插得到。

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

4.1 空間折射率內插精度分析

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

表1 不同高度角方向濕折射率內插值與MM5計算值比較統計結果

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

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

5 結 論

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

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

［2]張雙成，劉經南，葉世榕，等.顧及雙差殘差反演GPS信號方向的斜路徑水汽含量［J].武漢大學學報·信息科學版，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.網對水汽三維分布的監測及其在氣象學中的應用［D].上海：中國科學院上海天文臺，2004.

［6]劉焱雄，陳永奇.GPS氣象學中垂直干分量延時的精確確定［J].測繪學報，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]奚梅成，劉儒勛.數值分析方法［M].合肥：中國科學技術大學出版社，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遙感大氣可降水量及其在氣象中的應用研究［D].四川：西南交通大學，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.