GNSS水汽層析的約束條件方程變權代數重構算法

王 昊 丁 楠 張文淵 馮遵德

1 江蘇師范大學地理測繪與城鄉規劃學院，江蘇省徐州市上海路101號，221116 2 中國礦業大學環境與測繪學院，江蘇省徐州市大學路1號，221116

ART通過投影的方式使迭代解逐步逼近層析方程組的解算域。層析方程組除了包含GNSS信號的觀測方程外，通常還額外附加約束條件方程作為偽觀測方程，以解決層析方程組的不適定性問題[5]。約束條件方程在ART迭代初期具有恢復無信號體素塊內水汽參數的重要作用，但實際精度相對較低。隨著迭代的進行，所有無信號體素塊內的水汽參數逐漸恢復，此時迭代解若繼續以相同的投影修正權重向低精度約束條件超平面投影，容易造成迭代解偏離GNSS信號觀測方程組解算域的現象出現。為使ART的迭代解向精度更高的層析方程超平面投影，可引入高度角定權法控制GNSS信號觀測方程組的投影修正量[4]，而目前尚未提出控制附加約束條件方程組投影修正權重的方法。針對上述問題，本文提出一種約束條件變權的代數重構算法VWART，該算法以相鄰2次迭代過程中約束條件方程觀測值與重構值的殘差為依據，對該約束條件方程的投影權重進行實時優化，以降低因迭代解偏離GNSS真觀測方程組解算域而造成的誤差，提高ART迭代結果的精度和可靠性。同時，利用香港地區CORS網GNSS觀測數據進行水汽層析實驗，對約束條件方程權重不變的傳統ART和VWART進行對比驗證，并對2種算法反演的層析結果進行精度評定。

1 GNSS水汽層析的基本原理

GNSS水汽層析技術是根據層析區域內衛星信號傳播路徑上的斜路徑水汽含量(SWV)來反演該區域內的三維水汽場[6]，將對流層層析區域離散化，以每條信號的SWV為觀測值建立層析區域內信號傳播路徑長度、水汽密度參數與SWV的函數關系，即

(1)

式中，SWVi為第i條衛星信號傳播路徑上的斜路徑水汽含量，lij為第i條信號穿過第j個體素塊的直線長度，Dj為第j個體素塊內的待求水汽密度，n為總體素塊個數。式(1)模型構建的層析方程組的系數矩陣往往嚴重秩虧，因此一般通過附加約束條件的方法來解決層析方程組的不適定性問題[7]。常見的約束條件包括垂直約束和水平約束，其中垂直約束采用指數函數[8]，水平約束采用高斯加權函數[9]。結合觀測方程和約束條件方程的層析方程組可表示為：

(2)

式中，SWVObs為斜路徑水汽含量的觀測值向量；AObs、AV、AH分別為觀測值系數矩陣、垂直約束系數矩陣和水平約束系數矩陣；D為水汽密度參數向量；n為水汽密度參數個數；Δ為觀測噪聲；m為方程總個數，其數值等于各部分方程個數o、v、h之和。

2 GNSS水汽層析的約束條件變權代數重構算法

2.1 代數重構算法

ART運用迭代的方式獲得方程組的迭代解，以迭代初值為起點，依次向各方程所代表的多維超平面投影，實現迭代解向超平面相交解算域的逐步逼近，在滿足迭代終止條件時停止，便可輸出方程組的迭代解。該算法在解算水汽層析方程組時避免系數矩陣的求逆運算，迭代速度快、解算效率高，具有較高的穩定性和較強的抗噪性，即使在觀測方程較少、信號噪聲較大的不利條件下，也能獲得穩定且高精度的迭代解，因此被成功應用于電離層和對流層層析解算中[10]。GNSS水汽層析中傳統ART迭代公式為：

(3)

2.2 約束條件變權的代數重構算法

GNSS水汽層析的層析方程組一般由觀測方程組和約束條件方程組2部分構成。由于受到衛星和測站的幾何分布以及層析區域建模時格網劃分等因素的影響，層析區域中存在許多無衛星信號穿過的體素塊，約束條件方程組在ART迭代解算初期有助于恢復這些體素塊內的水汽參數。隨著迭代的進行，所有體素塊逐漸恢復水汽參數，此時約束條件方程組的作用將大幅降低，過強的約束條件反而會使迭代解偏離真觀測方程交會的解算域，進而造成較大誤差。

考慮到上述問題，本文在傳統ART的基礎上對約束條件方程組引入一種隨迭代實時變化的權重模型，提出一種新的約束條件變權代數重構算法VWART。VWART的迭代公式于觀測方程組部分(1≤i≤o，o為觀測方程組數量)見式(3)，于約束條件方程組部分(o+1≤i≤m，m為層析方程組總個數)為：

開店設鋪，做漁藥漁需物資買賣經營，常常需要那些看似不經意的妙招。在市場經營的日常細碎活動中，那些聰慧機智的漁需營銷店家往往新招迭出，靠著智慧做出自個兒生意亮點，收獲事半功倍的利好。

(4)

(5)

3 GNSS水汽層析實驗方案

本文利用香港CORS網中19個測站的GNSS觀測數據進行水汽層析實驗，層析區域水平范圍為113°50′38.4″～114°23′2.4″E、22°12′14.4″～22°33′50.4″N，水平分辨率為5′24″ ×5′24″[11]；垂直范圍以京士柏探空氣象站(HKKP)的探空水汽密度為依據，將水汽密度為0.1 g/m3的高度作為對流層層頂。2009～2018年共計10 a的8月探空水汽密度統計結果表明，對流層層頂平均高度為10 560 m。采用上疏下密的非均勻分層方式對垂直方向建模[6,12]，建模后將層析區域在垂直方向上劃分為13層，其中第1～6層的垂直分辨率為350 m，第7～13層的垂直分辨率分別為383 m、452 m、551 m、708 m、990 m、1 659 m、3 717 m。

研究時段為2019-08(共31 d)，屬于香港地區的夏季，氣候炎熱，受熱帶氣旋“韋帕”及“白鹿”的影響，共發生3次大規模降雨。實驗利用GAMIT10.7軟件對19個測站的GNSS原始觀測數據進行處理，解算策略如表1所示。

表1 GAMIT10.7解算策略

本文分別基于約束條件方程權重不變的傳統ART和本文提出的VWART方法設計2種方案進行層析實驗，解算水汽密度參數。ART中初值的設置是影響迭代解算結果的重要因素，高精度的迭代初值可同時提高迭代速度和層析結果質量。本文以HKKP探空站2009～2018年每年8月高精度探空數據獲得的水汽廓線信息作為迭代初值，利用原始探空數據精確計算出位于探空氣球高度處的水汽密度，對分布于各分層高度區間內的探空水汽密度取平均值，即可獲得各層水汽密度的迭代初值。利用HKKP探空站的原始探空數據計算HKKP站位置天頂方向的探空水汽密度ρRS，具體公式為：

(6)

式中，Rv=461.5 J·kg-1·K-1；T為溫度；e為水汽壓，其數值可由世界氣象組織推薦的Goff-Gratch公式計算得到[13]。

迭代終止的條件設置為：

‖xk+1-xk‖2<10-5

(7)

式中，xk為第k次的迭代解，當相鄰2次迭代解滿足式(7)時，迭代終止。

4 GNSS水汽層析實驗結果分析

4.1 層析解算結果的綜合分析

本文以原始探空數據計算獲得的探空水汽密度為參考值，對層析反演結果進行驗證。探空數據具有精度高、垂直分辨率高等特性，因此以探空數據作為參考標準對層析反演結果進行驗證切實可行[14-15]。由于探空數據具有低時間分辨率的局限性，且HKKP探空站的探空數據僅UTC 00:00和UTC 12:00兩個時刻可被獲取，故本文僅利用探空水汽密度對實驗時段內每日UTC 00:00～00:30及UTC 12:00～12:30兩個時段的層析結果進行分析。除探空數據資料外，本文還以ECMWF數據資料為標準對層析反演結果進行驗證。ECMWF的第4代產品ERA-Interim可提供全球再分析資料，其中包含溫度、相對濕度等用于計算水汽密度的格網數據資料，最小空間分辨率可達7′30″ ×7′30″[7,16]。表2(單位g/m3)為分別以探空數據和ECMWF數據為參考，傳統ART(方案A)和VWART(方案B)層析解算結果的RMSE(均方根誤差)、MAE(平均絕對誤差)和IQR(四分位區間)對比。

表2 2種方案層析解算結果的RMSE、MAE和IQR對比

由表2可知，無論以哪種數據為參考標準，方案B的RMSE、MAE和IQR均小于方案A。以探空數據為參考時，VWART對RMSE、MAE和IQR的改進率分別為20.334%、18.126%和11.995%；以ECMWF數據為參考時，VWART的改進率分別為36.625%、29.096%和25.299%，說明相較于傳統ART方法，VWART具有更高的層析解算精度。

圖1為2種實驗方案解算獲得的層析水汽密度與探空水汽密度對照散點圖。由圖可見，二者均呈橄欖形，其中方案B的散點相較于方案A表現出更強的聚集性。方案B散點的RMSE比方案A低19.558%，方案B的確定系數R2高于方案A，且方案B擬合曲線表達式的斜率也更接近于1。通過散點圖的對比進一步表明，VWART解算結果的精度明顯優于傳統ART。

圖1 2種方案層析水汽密度與探空水汽密度對照散點圖

圖2為2種方案層析解算結果殘差的正態分布曲線。可以看出，方案B殘差正態分布曲線的頂點位置高于方案A，說明方案B有更多的殘差值位于0附近；方案B殘差正態分布曲線具有更低的方差值和更接近于0的期望值，說明方案B的殘差分布更加集中。方案A有相當一部分的殘差分布于[-1，1]區間外，甚至有部分殘差分布于[-4，4]區間外；而方案B的大部分殘差分布于[-1，1]內。由此可知，方案B殘差分布情況的集中性和穩定性整體優于方案A，說明相比于傳統ART，VWART解算結果的整體穩定性和可靠性有十分顯著的提高。

圖2 2種方案層析結果殘差的正態分布曲線

4.2 各歷元、各高度的層析解算結果分析

圖3為2種方案各層析歷元解算結果的RMSE對比。實驗時段內共進行62組層析實驗，其中方案A的RMSE最大為3.074 g/m3，最小為0.501 g/m3；方案B的RMSE最大為2.473 g/m3，最小為0.410 g/m3。VWART層析解算結果的最大RMSE和最小RMSE相較于傳統ART分別改進19.551%和18.164%。

圖3 2種方案各層析歷元解算結果的RMSE對比

圖4和圖5進一步比較不同方案解算大氣水汽密度的垂直精度。由圖4可見，以探空數據計算的垂直水汽廓線為參考，方案B反演的垂直水汽廓線的RMSE和相對誤差在各高度層均低于方案A。在4 000 m以下的高度區間內，方案B對RMSE數值的降低尤為顯著；而在地表附近的高度區間內，方案B對RMSE的優化接近50%，進一步證實VWART的優越性。由圖5可見，方案B各高度層殘差分布的集中性優于方案A，僅第4、5、6三層出現較為明顯的異常值；而方案A在第1～6層均出現明顯的殘差異常值，且異常值數量遠多于方案B，其殘差分布的總體情況相較于方案B也更加疏松。由此可見，VWART大幅提高了各垂直高度處層析解算結果的穩定性和精度。

圖4 2種方案層析結果的分層RMSE和相對誤差

圖5 2種方案層析結果殘差的分層箱形圖

4.3 不同降水活動的層析解算結果分析

對流層大氣水汽含量與降水活動存在密切聯系，為深入分析VWART方法在不同降水活動條件下的適用性，本文將3種天氣條件下2種方案反演的層析水汽廓線與探空數據和ECMWF數據計算的水汽廓線進行對比分析。研究時段及天氣條件分別為：2019-08-05(無雨)、2019-08-14(小雨)、2019-08-25(暴雨)，反演水汽廓線所使用的數據由HKKP探空站位置、ECMWF格網點位置和體素塊位置共同決定。由于HKKP探空站位于114°6′50.4″～114°12′14.4″E、22°17′38.4″～22°23′2.4″N范圍內，故使用位于該區域的第1～13層體素塊層析結果進行層析水汽廓線的反演。同時，該區域僅存在一處位于114°7′30″E、22°22′30″N的ECMWF格網點，因此使用該格網點處的ECMWF數據進行層析水汽廓線的對比分析。圖6為3種天氣條件下層析水汽廓線與探空數據和ECMWF數據的對比。可以看出，不同天氣條件下的大氣水汽含量各不相同，整體表現為雨天大氣水汽含量高、無雨天大氣水汽含量低。2種方案均能正確反映出大氣水汽的總體分布情況，反演的層析水汽廓線中方案B與探空數據和ECMWF數據表現出更強的一致性。在無雨的天氣條件下，大氣水汽含量較少，方案A的垂直水汽廓線在2 000～6 000 m處與探空水汽廓線的差異遠大于方案B；在小雨的天氣條件下，方案B反演的層析水汽廓線與探空水汽廓線在1 000 m以下的高度區間內具有更好的貼合度；在暴雨的天氣條件下，方案A反演的層析水汽廓線與探空水汽廓線在3 000～6 000 m高度區間內存在較大差異，而方案B與探空水汽廓線的差異較小。綜上所述，VWART在不同天氣條件下均具有優秀的適用性，且相比于傳統ART，其更能提高大氣水汽的反演質量。

圖6 不同天氣條件下2種方案反演的層析水汽廓線與探空數據和ECMWF數據對比

5 結 語

1)VWART解算結果的精度和穩定性均優于傳統ART，在以探空數據為參考標準的情況下，其解算結果的RMSE、MAE和IQR相較于傳統ART分別降低20.334%、18.126%和11.995%；在以ECMWF數據為參考標準的情況下其RMSE、MAE和IQR分別降低36.625%、29.096%和25.299%。VWART的散點圖表明其具有更高的確定系數和更低的RMSE，解算結果的殘差分布更集中，殘差正態分布曲線擁有更低的方差和更接近于0的期望值。

2)VWART對各高度區間的層析解算質量均有所改進，對地表附近高度區間層析解算結果精度的提高尤為顯著。在研究時段的所有層析歷元中，VWART的RMSE最大值和最小值相較于傳統ART分別改進19.551%和18.164%。

3)與傳統ART反演的水汽廓線進行對比發現，VWART反演的層析水汽廓線與探空水汽廓線具有更強的一致性，表現出優于傳統ART的大氣水汽反演效果，且VWART在無雨、小雨和暴雨等不同天氣條件下均表現出優秀的大氣水汽重構能力。