王 昊 丁 楠 張文淵 馮遵德

1 江蘇師范大學(xué)地理測繪與城鄉(xiāng)規(guī)劃學(xué)院，江蘇省徐州市上海路101號(hào)，221116 2 中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇省徐州市大學(xué)路1號(hào)，221116

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

1 GNSS水汽層析的基本原理

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

(1)

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

(2)

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

2 GNSS水汽層析的約束條件變權(quán)代數(shù)重構(gòu)算法

2.1 代數(shù)重構(gòu)算法

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

(3)

2.2 約束條件變權(quán)的代數(shù)重構(gòu)算法

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

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

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

(4)

(5)

3 GNSS水汽層析實(shí)驗(yàn)方案

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

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

表1 GAMIT10.7解算策略

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

(6)

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

迭代終止的條件設(shè)置為：

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

(7)

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

4 GNSS水汽層析實(shí)驗(yàn)結(jié)果分析

4.1 層析解算結(jié)果的綜合分析

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

表2 2種方案層析解算結(jié)果的RMSE、MAE和IQR對(duì)比

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

圖1為2種實(shí)驗(yàn)方案解算獲得的層析水汽密度與探空水汽密度對(duì)照散點(diǎn)圖。由圖可見，二者均呈橄欖形，其中方案B的散點(diǎn)相較于方案A表現(xiàn)出更強(qiáng)的聚集性。方案B散點(diǎn)的RMSE比方案A低19.558%，方案B的確定系數(shù)R2高于方案A，且方案B擬合曲線表達(dá)式的斜率也更接近于1。通過散點(diǎn)圖的對(duì)比進(jìn)一步表明，VWART解算結(jié)果的精度明顯優(yōu)于傳統(tǒng)ART。

圖1 2種方案層析水汽密度與探空水汽密度對(duì)照散點(diǎn)圖

圖2為2種方案層析解算結(jié)果殘差的正態(tài)分布曲線?？梢钥闯?，方案B殘差正態(tài)分布曲線的頂點(diǎn)位置高于方案A，說明方案B有更多的殘差值位于0附近；方案B殘差正態(tài)分布曲線具有更低的方差值和更接近于0的期望值，說明方案B的殘差分布更加集中。方案A有相當(dāng)一部分的殘差分布于[-1，1]區(qū)間外，甚至有部分殘差分布于[-4，4]區(qū)間外；而方案B的大部分殘差分布于[-1，1]內(nèi)。由此可知，方案B殘差分布情況的集中性和穩(wěn)定性整體優(yōu)于方案A，說明相比于傳統(tǒng)ART，VWART解算結(jié)果的整體穩(wěn)定性和可靠性有十分顯著的提高。

圖2 2種方案層析結(jié)果殘差的正態(tài)分布曲線

4.2 各歷元、各高度的層析解算結(jié)果分析

圖3為2種方案各層析歷元解算結(jié)果的RMSE對(duì)比。實(shí)驗(yàn)時(shí)段內(nèi)共進(jìn)行62組層析實(shí)驗(yàn)，其中方案A的RMSE最大為3.074 g/m3，最小為0.501 g/m3；方案B的RMSE最大為2.473 g/m3，最小為0.410 g/m3。VWART層析解算結(jié)果的最大RMSE和最小RMSE相較于傳統(tǒng)ART分別改進(jìn)19.551%和18.164%。

圖3 2種方案各層析歷元解算結(jié)果的RMSE對(duì)比

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

圖4 2種方案層析結(jié)果的分層RMSE和相對(duì)誤差

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

4.3 不同降水活動(dòng)的層析解算結(jié)果分析

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

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

5 結(jié) 語

1)VWART解算結(jié)果的精度和穩(wěn)定性均優(yōu)于傳統(tǒng)ART，在以探空數(shù)據(jù)為參考標(biāo)準(zhǔn)的情況下，其解算結(jié)果的RMSE、MAE和IQR相較于傳統(tǒng)ART分別降低20.334%、18.126%和11.995%；在以ECMWF數(shù)據(jù)為參考標(biāo)準(zhǔn)的情況下其RMSE、MAE和IQR分別降低36.625%、29.096%和25.299%。VWART的散點(diǎn)圖表明其具有更高的確定系數(shù)和更低的RMSE，解算結(jié)果的殘差分布更集中，殘差正態(tài)分布曲線擁有更低的方差和更接近于0的期望值。

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

3)與傳統(tǒng)ART反演的水汽廓線進(jìn)行對(duì)比發(fā)現(xiàn)，VWART反演的層析水汽廓線與探空水汽廓線具有更強(qiáng)的一致性，表現(xiàn)出優(yōu)于傳統(tǒng)ART的大氣水汽反演效果，且VWART在無雨、小雨和暴雨等不同天氣條件下均表現(xiàn)出優(yōu)秀的大氣水汽重構(gòu)能力。