地基GPS水汽層析的投影面算法

丁　楠，張書畢,2

1. 中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，江蘇 徐州 221008； 2. 江蘇省資源環(huán)境信息工程重點(diǎn)實(shí)驗(yàn)室，江蘇 徐州 221008

?

地基GPS水汽層析的投影面算法

丁楠1，張書畢1,2

1. 中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，江蘇 徐州 221008； 2. 江蘇省資源環(huán)境信息工程重點(diǎn)實(shí)驗(yàn)室，江蘇 徐州 221008

Foundationsupport：TheNationalNaturalScienceFoundationofChina(No. 41504032);TheNationalNaturalScienceFoundationofJiangsuProvinceofChina(No.BK20150175);TheSpecializedResearchFundfortheDoctoralProgramofHigherEducationofChina(No. 20130095110022)

水汽層析技術(shù)在研究氣候變化、極端天氣預(yù)警、輔助數(shù)值天氣預(yù)報(bào)等方面發(fā)揮著重要作用。常規(guī)的水汽層析技術(shù)在計(jì)算層析方程系數(shù)時(shí)，需要進(jìn)行大量的求交運(yùn)算。對(duì)此，本文提出了一種投影面算法，提高了運(yùn)算速度和反演的精度。本文采用的試驗(yàn)數(shù)據(jù)是香港衛(wèi)星定位參考站網(wǎng)(SatRef)提供的GPS信號(hào)數(shù)據(jù)，與傳統(tǒng)算法相比，投影面算法的計(jì)算速度更快、計(jì)算量更小，不受層析區(qū)域大小影響，水汽層析的結(jié)果與探空數(shù)據(jù)具有良好的一致性。

水汽層析; 信號(hào)路徑; 像素; 體素遍歷算法; 代數(shù)重構(gòu)算法

水汽是描述大氣狀況的基本特征參數(shù)，大氣中的水汽含量只占全球水循環(huán)系統(tǒng)總量的百分之一左右，卻是其中時(shí)空變化最快的活躍因子。大氣水汽分布及其變化特征的觀測(cè)，對(duì)于研究氣候變化有著重要意義。常規(guī)的大氣水汽探測(cè)手段如無線電探空，觀測(cè)間隔一般為12h，時(shí)間并不連續(xù)，且觀測(cè)站點(diǎn)相隔一般都在200～300km。地基微波輻射計(jì)在有濃云時(shí)穿透能力下降，特別是有降水發(fā)生時(shí)更會(huì)產(chǎn)生較大誤差。這些都限制了高時(shí)空分辨率、高精度的水汽信息獲取。地基GPS水汽層析技術(shù)補(bǔ)充了常規(guī)觀測(cè)方法的不足，可以獲取高時(shí)空分辨率，探測(cè)精度1～2mm[1]的水汽資料。

地基GPS水汽層析技術(shù)由文獻(xiàn)[2]首先提出，它給出了水汽層析方程的水平、垂直約束條件以及層析區(qū)域垂直方向分層的理論依據(jù)。而在層析區(qū)域網(wǎng)格劃分方面，文獻(xiàn)[3—4]對(duì)垂直方向分層及其分辨率進(jìn)行了專門的研究。通過層析結(jié)果與數(shù)值天氣預(yù)報(bào)模型[5]以及探空數(shù)據(jù)[6]的比對(duì)，證實(shí)層析技術(shù)的可靠性，而層析結(jié)果與這兩個(gè)系統(tǒng)之間存在的差異性[7]是由于觀測(cè)站網(wǎng)形的幾何結(jié)構(gòu)不同造成的。長(zhǎng)期研究還表明層析模型在垂直分辨率[8-10]以及探測(cè)高層水汽(4000m以上)[11]方面有著局限性。目前可以通過增加COSMIC(掩星)數(shù)據(jù)[12]、地面氣象觀測(cè)數(shù)據(jù)[13-14]等來改善層析在垂直分辨率上的不足。

水汽層析技術(shù)的原理是通過地基GPS網(wǎng)絡(luò)以及高精度GPS數(shù)據(jù)處理軟件獲取濕延遲觀測(cè)值，濕延遲和傾斜路徑上的水汽含量(slantwatervapor，SWV)之間可以進(jìn)行轉(zhuǎn)換，用于區(qū)域水汽場(chǎng)重建。將地基GPS網(wǎng)天頂方向上的空間，在水平和垂直方向上離散化為包含水汽信息的體素；通過SWV構(gòu)建層析方程組就可以解算出體素內(nèi)的水汽參數(shù)。而包含濕延遲觀測(cè)值的GPS信號(hào)穿過對(duì)流層時(shí)，受到對(duì)流層的折射影響，GPS信號(hào)要發(fā)生彎曲和延遲，其中信號(hào)的彎曲量很小[15]，由此引起的路徑總延遲(slanttotaldelay，STD)可以忽略[16]。所以在進(jìn)行GPS水汽層析時(shí)，一般將信號(hào)傳播路徑簡(jiǎn)化為某歷元時(shí)刻，衛(wèi)星位置到GPS接收機(jī)天線位置的空間直線路徑。

建立層析觀測(cè)方程組是層析技術(shù)的核心。構(gòu)建方程組就要確定每條GPS信號(hào)穿過對(duì)流層離散化后的體素位置，獲取信號(hào)穿過這些體素的截距長(zhǎng)度，由體素的位置坐標(biāo)和對(duì)應(yīng)體素內(nèi)信號(hào)長(zhǎng)度確定方程的系數(shù)矩陣。傳統(tǒng)算法[17]是求出信號(hào)路徑與層析區(qū)域內(nèi)所有面的交點(diǎn)，通過求出兩點(diǎn)之間的距離確定對(duì)應(yīng)體素內(nèi)信號(hào)路徑的長(zhǎng)度。傳統(tǒng)算法無法準(zhǔn)確地找出層析區(qū)域內(nèi)和信號(hào)相交的平面，而需要對(duì)所有平面求交后才能判定，增加了不必要的計(jì)算量，而當(dāng)層析區(qū)域增大時(shí)，每條信號(hào)的求交運(yùn)算也都會(huì)相應(yīng)增加。

針對(duì)傳統(tǒng)方法存在的上述問題，本文提出了一種投影面算法。該算法在體素大小確定后，即使擴(kuò)大層析區(qū)域也不會(huì)增加任何計(jì)算量，而且可以準(zhǔn)確找出層析區(qū)域內(nèi)與信號(hào)相交的平面進(jìn)行求交運(yùn)算，避免了計(jì)算的冗余。新算法同時(shí)采用代數(shù)重構(gòu)算法(algebraicreconstructiontechnique，ART)進(jìn)行層析反演，并根據(jù)水汽層析的特性對(duì)方程組進(jìn)行分組，通過按次序分組將迭代值投影到各組超平面，提高了反演的迭代速度和精度。本文采用香港衛(wèi)星定位參考站網(wǎng) (SatRef)提供的GPS信號(hào)數(shù)據(jù)進(jìn)行試驗(yàn)。采集9個(gè)香港參考站不同時(shí)段的觀測(cè)數(shù)據(jù)，結(jié)果表明投影面算法的計(jì)算速度更快、計(jì)算量更小并且不受層析區(qū)域大小影響。

1　投影面算法

投影面算法主要包括基于投影面輔助的直線遍歷體素和分組投影ART兩部分，直線遍歷體素算法減少了計(jì)算求交平面的數(shù)量，而分組投影ART則提高了層析反演的迭代速度以及結(jié)果精度。

1.1基于投影面輔助的直線遍歷體素

圖1　直線遍歷像素Fig.1　Line traverse pixels

首先定義某一點(diǎn)高度為該點(diǎn)與起點(diǎn)A所在像素底邊的垂線距離；Gi定義為當(dāng)前像素右上角點(diǎn)、ti是過Gi點(diǎn)的豎線與直線AB的交點(diǎn)。從起點(diǎn)A(x1,y1)所在的當(dāng)前像素(X,Y)開始判斷直線經(jīng)過的像素，當(dāng)直線AB位于當(dāng)前像素(X,Y)右上角G1點(diǎn)之上時(shí)，直線必然經(jīng)過與當(dāng)前像素相鄰的(X,Y+1)和(X+1,Y+1)。并且(X+1,Y+1)像素就是下一步的當(dāng)前體素；而若如圖1所示情況，當(dāng)直線AB在G1點(diǎn)之下時(shí)，直線AB經(jīng)過當(dāng)前像素右邊的(X+1,Y)，選擇該像素并作為下一步的當(dāng)前像素，至于此時(shí)直線是否還經(jīng)過(X+1,Y+1)則需進(jìn)一步判斷。因此算法在每一步遍歷時(shí)分兩種路徑：直線在參照點(diǎn)Gi之上：(X,Y)→(X,Y+1)→(X+1,Y+1)和直線在參照點(diǎn)之下：(X,Y)→(X+1,Y)。

當(dāng)判定直線穿過像素的路徑后，當(dāng)前像素變?yōu)槁窂缴献詈笠粋€(gè)到達(dá)的像素。再按照上述的方法在下一個(gè)單元確定直線穿過像素的路徑。判斷直線與參照點(diǎn)之間的上下關(guān)系，是二維直線像素遍歷的關(guān)鍵。從起點(diǎn)像素底邊算起，G1點(diǎn)的高度是Yv，而直線AB與過G1點(diǎn)豎線交點(diǎn)t1的初始高度為

(1)

HGi=n·Yv

(2)

(3)

上面是討論Yv/Xv大于dy/dx的情況，還有7種情況與其類似，這7種情況對(duì)應(yīng)的參照點(diǎn)Gi高度是Xv或Yv的整數(shù)倍，而對(duì)應(yīng)的ti點(diǎn)高度則各不相同。表1是直線在二維平面內(nèi)遍歷像素的8種情況分類及對(duì)應(yīng)ti點(diǎn)高度的通式。

表1　直線在二維平面內(nèi)遍歷像素的8種情況

當(dāng)Yv/Xv>dy/dx并且Zv/Xv>dz/dx時(shí)，X軸判定為主坐標(biāo)軸。下面空間直線的遍歷分兩大類情況討論，第1類，當(dāng)沿直線從A點(diǎn)所在當(dāng)前體素(X,Y,Z)穿過下一個(gè)體素時(shí)，X坐標(biāo)總是增加或減少一個(gè)體素單元；若此時(shí)算法判定Y和Z坐標(biāo)只有其中一個(gè)坐標(biāo)變化(增加或減少一個(gè)體素單位)或者不變化，則由二維直線遍歷像素算法的擴(kuò)展就可以解決。Y和Z坐標(biāo)不變化的情況很簡(jiǎn)單，這里詳細(xì)描述Y和Z坐標(biāo)只有其中一個(gè)坐標(biāo)變化的情況，算法思想是由XY平面和XZ平面的二維直線遍歷像素算法獲得空間直線遍歷體素的位置。如圖2(a)所示，X軸為主軸，直線起點(diǎn)A在體素(X,Y,Z)中。直線AB與體素的交點(diǎn)1、2、3，分別在角a、b、c對(duì)應(yīng)的平面上。將直線AB和體素分別投影到XY和XZ平面，在對(duì)應(yīng)的投影面用二維直線的像素遍歷算法分別計(jì)算投影直線經(jīng)過像素的路徑。圖2(b)中直線AB在XY平面投影為A1B1，由二維直線遍歷像素算法可知A1B1經(jīng)過像素的路徑為：(X,Y)→(X+1,Y)。圖2(c)中直線在XZ平面投影為A2B2，其經(jīng)過像素路徑為：(X,Z)→(X,Z+1)→(X+1,Z+1)。兩個(gè)投影平面路徑的當(dāng)前像素和終點(diǎn)像素分別對(duì)應(yīng)相同的體素，XY平面的路徑實(shí)際上應(yīng)該是(X,Y)→(X,Y)→(X+1,Y)，由于前兩步路徑所在像素相同，所以并未反映在二維投影中。如圖2(a)所示，參照兩個(gè)投影面的路徑可以得到三維空間直線AB經(jīng)過體素的路徑為：(X,Y,Z)→(X,Y,Z+1)→(X+1,Y,Z+1)。

第2類情況是當(dāng)Y和Z坐標(biāo)在一步之內(nèi)都變化(增加或減少)。如圖3(a)所示，直線穿過角a、b、c、d所在平面的交點(diǎn)為1、2、3、4?？芍?dāng)前體素(X,Y,Z)和完成一次遍歷后，到達(dá)的下個(gè)當(dāng)前體素是(X+1,Y+1,Z+1)。圖3(b)中二維直線遍歷像素算法可知A1B1經(jīng)過像素的路徑為(X,Y)→(X,Y+1)→(X+1,Y+1)。圖3(b)中直線經(jīng)過像素路徑為(X,Z)→(X,Z+1)→(X+1,Z+1)。第2類情況中，必然經(jīng)過的體素為：當(dāng)前體素(X,Y,Z)、(X,Y+1,Z+1)和下一個(gè)當(dāng)前體素(X+1,Y+1,Z+1)；而從(X,Y,Z)到達(dá)(X,Y+1,Z+1)存在兩種不同的選擇，因此存在兩種不同的路徑從當(dāng)前體素到達(dá)下一個(gè)當(dāng)前體素,路徑L1：(X,Y,Z)→(X,Y+1,Z)→(X,Y+1,Z+1)→(X+1,Y+1,Z+1)和路徑L2：(X,Y,Z)→(X,Y,Z+1)→(X,Y+1,Z+1)→(X+1,Y+1,Z+1)。兩條路徑的區(qū)別在第2個(gè)遍歷的體素上，直線是經(jīng)過體素(X,Y+1,Z)還是(X,Y,Z+1)，通過本文提出的投影面輔助方法可以找出正確的路徑。當(dāng)遇到第2類情況時(shí)，需判斷(X,Y+1,Z)和(X,Y,Z+1)哪一個(gè)是直線穿過的體素。此時(shí)將圖3(a)投影到兩個(gè)次軸Y、Z組成的YZ平面上。如圖3(d)所示，通過二維直線遍歷體素算法可知投影直線A3B3經(jīng)過的路徑為：(Y,Z)→(Y+1,Z)→(Y+1,Z+1)，而體素(X,Y+1,Z)和(X,Y,Z+1)在YZ平面上投影為(Y+1,Z)、(Y,Z+1)；與對(duì)應(yīng)的正確路徑對(duì)比可發(fā)現(xiàn)并沒有經(jīng)過像素(Y,Z+1)，因此路徑2是錯(cuò)誤的，去除體素(X,Y,Z+1)，正確的路徑為L(zhǎng)1。

圖2　第1類空間直線遍歷體素Fig.2　First type of line traverse voxel

圖3　第2類空間直線遍歷體素Fig.3　Second type of line traverse voxel

1.1.1求交平面計(jì)算優(yōu)化

采用投影面輔助直線遍歷體素算法,可以準(zhǔn)確獲取信號(hào)穿過的平面,雖然同傳統(tǒng)算法相比,計(jì)算量有所減小，但還是需要較多的計(jì)算量。因此對(duì)該算法進(jìn)行優(yōu)化，從而進(jìn)一步減少計(jì)算量。投影面算法在獲取觀測(cè)方程系數(shù)即體素內(nèi)信號(hào)截距長(zhǎng)度時(shí)，除起點(diǎn)所在體素內(nèi)的截距需要單獨(dú)計(jì)算外，其他分為兩大類情況。首先按照信號(hào)穿過的體素坐標(biāo)中的Z坐標(biāo)進(jìn)行分層，第1類情況如圖4(a)所示，體素的Z坐標(biāo)不重復(fù)，體素(X,Y,Z)與(X,Y,Z+1)內(nèi)的信號(hào)截距長(zhǎng)度一樣，都是信號(hào)穿過體素底面和頂面兩個(gè)交點(diǎn)之間的長(zhǎng)度，因此該情況無論直線穿過多少體素，只需計(jì)算一次截距長(zhǎng)度。第2類情況，Z坐標(biāo)值重復(fù)，將Z坐標(biāo)數(shù)值相同的，按照遍歷順序分為一組，每組除必定穿過的對(duì)應(yīng)體素的起始面(底面)和最終面(頂面)外，還需確定信號(hào)穿過的側(cè)面。如圖4(b)，直線遍歷體素順序?yàn)?X,Y,Z)→(X-1,Y,Z)→(X-1,Y,Z+1)→(X-1,Y+1,Z+1)→(X-2,Y+1,Z+1)，按照情況可分組為(X,Y,Z)→(X-1,Y,Z)以及(X-1,Y,Z+1)→(X-1,Y+1,Z+1)→(X-2,Y+1,Z+1)。第1組中，對(duì)應(yīng)底面與頂面分別為z1=(Z-1)·Zv和z3=Z·Zv，其下標(biāo)表示直線穿過該平面的次序，(X,Y,Z)與(X-1,Y,Z)不同的坐標(biāo)為X1=X與X2=X-1，則穿過的側(cè)面可通過式(4)計(jì)算

x2=min(X1,X2)·Xv

(4)

式中，min(X1,X2)表示取X1與X2中較小的數(shù)值，Xv為體素與X軸平行的邊長(zhǎng)。

圖4　計(jì)算截距長(zhǎng)度的方法Fig.4　Method for calculating the length of the ray across voxels

同理第2組有3個(gè)體素，共穿過2個(gè)側(cè)面。對(duì)照上述方法可知，第1個(gè)側(cè)面為y4=Y·Yv，第2個(gè)側(cè)面為x5=(X-2)·Xv。第2組底面為第1組的頂面，無需計(jì)算，其頂面為z6=(Z+1)·Zv。因此第2類情況，可以通過直線方程與平面方程z1、x2、z3、y4、x5、z6進(jìn)行求交運(yùn)算，獲取交點(diǎn)坐標(biāo)后，計(jì)算體素內(nèi)截距長(zhǎng)度，而無需對(duì)圖4(b)中相關(guān)的10個(gè)平面全部進(jìn)行求交運(yùn)算。由于香港參考站位置密集，層析區(qū)域較小，因此層析區(qū)域內(nèi)的信號(hào)高度角普遍較高，信號(hào)穿過的體素多屬于第1種情況，使用本文的方法可以顯著提高計(jì)算速度。

1.2分組投影ART

代數(shù)重構(gòu)算法的基本原理是從給定的初值x0出發(fā)(上標(biāo)表示迭代次數(shù))，通過式(5)將x0按順序依次投影到方程代表的超平面H上，通過連續(xù)的投影，到達(dá)所有超平面的交點(diǎn)位置，該位置表示方程組的解

(5)

由于衛(wèi)星星座測(cè)站分布等引起的層析方程組解算病態(tài)問題，需要增加約束條件(如水平、垂直約束方程組)減小方程條件數(shù)以及避免解的多值性。觀測(cè)方程組與約束方程組構(gòu)成了GNSS水汽層析方程組

(6)

SWV=∫sHabsds

(7)

式中，s表示信號(hào)路徑；Habs表示絕對(duì)濕度。

代數(shù)重構(gòu)算法進(jìn)行迭代計(jì)算時(shí)，方程的投影次序會(huì)對(duì)最終迭代解的精度和迭代收斂速度造成影響。本節(jié)按照4類方程組的可靠性進(jìn)行分組，將信息量大的觀測(cè)方程組排在第一組，使解的結(jié)果偏向觀測(cè)方程組提供的信息，有利于精確、高效地重建層析區(qū)域；輔助觀測(cè)方程組和垂直約束方程組雖然都具有較高的精度，但一般包含信息量較少，不利于獲取準(zhǔn)確的全局信息，這兩種方程組的排列順序根據(jù)它們獲取信息方式和數(shù)量不同可以調(diào)換分組順序；水平約束方程組的加入，雖然可以恢復(fù)沒有信號(hào)穿過體素塊內(nèi)的水汽參數(shù)，但在一定程度上降低了層析結(jié)果的精度，因此放在最后。分組投影ART算法將方程組按照OSVH順序進(jìn)行排序，并利用式(5)進(jìn)行迭代求解，在滿足停止條件后獲得最終的迭代解。

2　地基GPS信號(hào)處理實(shí)例與分析

2.1試驗(yàn)方案

本節(jié)利用投影面算法對(duì)GPS數(shù)據(jù)進(jìn)行處理,繪制信號(hào)穿過體素的射線圖。對(duì)應(yīng)2014年3月25日，UTC 6時(shí)的層析區(qū)域范圍為113.841°E—114.425°E，22.07°N—22.618°N，層析最高層為10 km。選取香港衛(wèi)星定位參考站網(wǎng)(SatRef)9個(gè)測(cè)站30 min的觀測(cè)數(shù)據(jù)，數(shù)據(jù)采樣間隔5 min。研究區(qū)域沿東西及天頂方向劃分為4×4×10共160個(gè)體素。香港9個(gè)參考站的層析網(wǎng)格劃分和GPS信號(hào)遍歷體素示意圖。

由于在不同時(shí)段,GPS信號(hào)接收數(shù)量不同,而GPS信號(hào)的密集程度以及分布情況對(duì)層析區(qū)域大小和位置有所影響。在信號(hào)密度較低時(shí)(信號(hào)和區(qū)域內(nèi)體素的比例)應(yīng)適當(dāng)縮小層析區(qū)域，而當(dāng)信號(hào)密度較高時(shí)，可擴(kuò)大層析區(qū)域來獲取更多的水汽信息。圖6對(duì)應(yīng)2014年3月25日，UTC 0時(shí)(圖6(a))、6時(shí)(圖6(b))、12時(shí)(圖6(c))GPS信號(hào)遍歷俯視圖。黑框?qū)嵕€內(nèi)的區(qū)域位置大小是相同的，但隨著時(shí)間的變化，理想的層析區(qū)域(黑色虛線方框內(nèi)區(qū)域)發(fā)生了變化。UTC 0時(shí)由于該時(shí)段接受衛(wèi)星信號(hào)較少，層析區(qū)域較小。UTC 6時(shí)和UTC 12時(shí)，由于這兩個(gè)時(shí)段接收的衛(wèi)星信號(hào)較多且分布均勻，所以相應(yīng)的層析區(qū)域擴(kuò)大。

2.2投影面算法對(duì)計(jì)算求交平面速度的分析

為驗(yàn)證投影面算法的速度，將其與傳統(tǒng)算法在圖6所示的3個(gè)時(shí)段的求交計(jì)算進(jìn)行比較。表2通過同時(shí)存在于3個(gè)時(shí)段層析區(qū)域內(nèi)的一條模擬信號(hào)，利用傳統(tǒng)算法和投影面算法計(jì)算其求交平面數(shù)量以及所需計(jì)算時(shí)間，證明了投影面算法隨層析區(qū)域擴(kuò)大并不影響其速度也不會(huì)增加計(jì)算量，而傳統(tǒng)算法會(huì)隨層析區(qū)域的擴(kuò)大而增加計(jì)算量和計(jì)算時(shí)間。

圖5　GPS信號(hào)處理Fig.5　GPS network for traversing voxel

圖6　GPS信號(hào)分布俯視圖Fig.6　Ground-based GPS signal distribution in the looking-down perspective

表3通過對(duì)3個(gè)時(shí)段所有層析區(qū)域內(nèi)GPS信號(hào)，在兩種算法下計(jì)算的求交平面?zhèn)€數(shù)和計(jì)算速度進(jìn)行比較，證明了投影面算法與傳統(tǒng)算法相比，具有更快的計(jì)算速度，并且通過1.1.1節(jié)提出的對(duì)求交平面分類的優(yōu)化方法，減少了求交平面的數(shù)量。

表2兩種算法一條信號(hào)計(jì)算量與計(jì)算速度對(duì)比

Tab.2ComparisonofcalculatedamountandspeedofoneGPSsignalforbothalgorithms

算法一條信號(hào)計(jì)算量(求交平面?zhèn)€數(shù))一條信號(hào)計(jì)算速度/sUTC0UTC6UTC12UTC0UTC6UTC12傳統(tǒng)算法1920220.7250.8010.861投影面算法4440.1860.1860.186

表3兩種算法的計(jì)算量與計(jì)算速度對(duì)比

Tab.3Comparisonofcalculatedamountandspeedforbothalgorithms

算法所有信號(hào)計(jì)算量(求交平面?zhèn)€數(shù))所有信號(hào)計(jì)算速度/sUTC0UTC6UTC12UTC0UTC6UTC12傳統(tǒng)算法92151642018216330.582618.506683.216投影面算法8922169249532.11977.81189.493

2.3投影面算法對(duì)水汽反演精度的分析

試驗(yàn)按照1.2節(jié)提出的分組表示方法，選擇OVSH、VOSH、SOVH、HSVO4種不同分組方案，利用2014年3月25日，UTC0時(shí)HKKP氣象站數(shù)據(jù)與4種方案計(jì)算的對(duì)應(yīng)體素內(nèi)水汽參數(shù)進(jìn)行比較分析，證明了分組投影ART算法的合理性。圖7(a)是4種方案與探空數(shù)據(jù)(RS)擬合的水汽廓線對(duì)比，從圖中可以看出投影面算法采用的OVSH與RS數(shù)據(jù)的水汽廓線基本一致，擬合水汽廓線最大誤差出現(xiàn)在2130m(0.853 2g/m3)，主要原因是水汽在1500～2500m之間出現(xiàn)違反指數(shù)遞減規(guī)律的波動(dòng)。圖7(b)是不同方案計(jì)算的各層平均水汽值同對(duì)應(yīng)RS數(shù)據(jù)做差獲得的誤差曲線圖，方案OVSH對(duì)應(yīng)的誤差曲線基本在0附近波動(dòng)，而其他幾個(gè)方案則誤差較大，HSVO方案由于迭代結(jié)果偏向前3組約束方程提供的先驗(yàn)信息，導(dǎo)致層析結(jié)果與實(shí)際不符，SOVH與VOSH層析誤差也有較大波動(dòng)，從層析結(jié)果來看，方程組的順序?qū)游龇囱菥染哂兄匾绊憽?/p>

圖8是各組方案的RMSE曲線圖，橫軸表示迭代次數(shù)，反映迭代收斂速度；縱軸表示RMSE，反映迭代結(jié)果的精度。結(jié)合表4各方案RMSE以及迭代次數(shù)，可以看出與水汽廓線對(duì)應(yīng)，OVSH的RMSE值最小(0.486 8g/m3)，HSVO的精度最低(3.799 5g/m3)。本次試驗(yàn)采用的垂直約束方程組V利用經(jīng)驗(yàn)?zāi)Ｐ瞳@取，精度較低。而輔助觀測(cè)方程組S是利用Rinexmet獲取的地表氣象數(shù)據(jù)，精度較高。因此采用SOVH方案使結(jié)果偏向S提供的信息時(shí)，要比偏向V的VOSH方案精度高。而HSVO解算結(jié)果精度最低，主要是由于H組成的約束條件自身含有較大誤差，但其迭代次數(shù)只有26次，因?yàn)橄冗M(jìn)行迭代的HSV方程組包含水汽信息更全面，使迭代結(jié)果更易快速收斂，但無法保證結(jié)果的精度。

圖7　水汽層析結(jié)果Fig.7　Results of Water Vapor Tomography

圖8　不同方案RMSE值Fig.8　RMSE for different schemes

分組方案OVSHSOVHVOSHHSVORMSE/(g/cm3)0.48680.70890.99683.7995迭代次數(shù)38414726

3　結(jié)　論

本文通過對(duì)直線遍歷體素時(shí)的輔助投影面以及迭代計(jì)算時(shí)投影到超平面的順序進(jìn)行研究，根據(jù)水汽層析的特性提出了投影面算法。與傳統(tǒng)算法相比，利用1.1.1節(jié)優(yōu)化方法，投影面算法在計(jì)算求交平面時(shí)，計(jì)算速度更快、計(jì)算量更小，并且不受層析區(qū)域大小影響，計(jì)算效率上的優(yōu)勢(shì)會(huì)隨著層析區(qū)域的擴(kuò)大而更加明顯，隨著各國(guó)GPS網(wǎng)的加密與建設(shè)，GPS水汽層析的區(qū)域勢(shì)必會(huì)進(jìn)一步擴(kuò)大，本文提出的直線遍歷體素算法可以高效地獲取層析的系數(shù)矩陣，為今后進(jìn)行實(shí)時(shí)大區(qū)域地基GPS水汽層析做好準(zhǔn)備。而在水汽層析反演計(jì)算方面，新算法通過分組排序的方式，使迭代解的結(jié)果偏向觀測(cè)方程組提供的可靠信息，提高了ART算法的精度，反演結(jié)果與HKKP提供的探空數(shù)據(jù)具有良好的一致性。試驗(yàn)分析發(fā)現(xiàn)優(yōu)先對(duì)包含完備信息的方程組(如約束方程組)進(jìn)行迭代計(jì)算，雖然可以達(dá)到快速收斂的目的，但若該組方程誤差較大，先驗(yàn)信息不準(zhǔn)確，其結(jié)果也會(huì)偏向誤差較大的方程提供的信息，最終降低解算的結(jié)果。由于ART算法本質(zhì)是通過逐個(gè)迭代方程獲得迭代解，因此除按照方程組性質(zhì)進(jìn)行分組外，對(duì)分組后每個(gè)層析方程的排序也是今后需要研究的問題之一。

[1]JIANGP,YESR,LIUYY,etal.NearReal-timeWaterVaporTomographyUsingGround-basedGPSandMeteorologicalData:Long-termExperimentinHongKong[J].AnnalesGeophysicae, 2014, 32(8): 911-923.

[2]FLORESA,RUFFINIG,RIUSA. 4DTroposphericTomographyUsingGPSEstimatedDelays[C]∥ProceedingsoftheIEEE1999InternationalGeoscienceandRemoteSensingSymposium, 1999.Hamburg:IEEE, 1999, 1: 675-677.

[3]SHRESTHASM.InvestigationsintotheEstimationofTroposphericDelayandWetRefractivityUsingGPSMeasurements[R].Calgary,Alberta:UniversityofCalgary, 2003.

[4]HOYLEVA.DataAssimilationfor4-DWetRefractivityModellinginaRegionalGPSNetwork[R].Calgary,Alberta:UniversityofCalgary, 2005.

[5]PANYT,PESECP,STANGLG.AtmosphericGPSSlantPathDelaysandRayTracingthroughNumericalWeatherModels:AComparison[J].PhysicsandChemistryoftheEarth,PartA:SolidEarthandGeodesy, 2001, 26(3): 183-188.

[6]BIYanmeng,MAOJietai,LIChengcai.PreliminaryResultsof4-DWaterVaporTomographyintheTroposphereUsingGPS[J].AdvancesinAtmosphericSciences, 2006, 23(4): 551-560.

[7]TROLLERM,LEUENBERGERD,BROCKMANNE,etal.GPS-tomography:ResultsandAnalysesoftheOperationalDeterminationofHumidityProfilesoverSwitzerland[J].GeophysicalResearchAbstracts,2007,9:11-18.

[8]BENDERM,RAABEA.PreconditionstoGroundBasedGPSWaterVapourTomography[J].AnnalesGeophysicae, 2007, 25(8): 1727-1734.

[9]BOCKO,DOERFLINGERE,MASSONF,etal.GPSWaterVaporProjectAssociatedtotheESCOMPTEProgramme:DescriptionandFirstResultsoftheFieldExperiment[J].PhysicsandChemistryoftheEarth, 2004, 29(2-3): 149-157.

[10]PERLERD,GEIGERA,HURTERF. 4DGPSWaterVaporTomography:NewParameterizedApproaches[J].JournalofGeodesy, 2011, 85(8): 539-550.

[11]CHAMPOLLIONC,MASSONF,BOUINMN,etal.GPSWaterVapourTomography:PreliminaryResultsfromtheESCOMPTEFieldExperiment[J].AtmosphericResearch, 2005, 74(1-4): 253-274

[12]XIAP,CAIC,LIUZ.GNSSTroposphereTomographyBasedonTwo-stepReconstructionsUsingGPSObservationsandCOSMICProfiles[J].AnnalesGeophysicae, 2013, 31(10): 1805-1815.

[13]ROHMW,BOSYJ.TheQualityofMeteorologicalObservationsandTroposphericDelayfromEPN/IGSPermanentStationsLocatedintheSudetyMountainsandintheAdjacentAreas[J].ActaGeodynamicaetGeomaterialia, 2007, 4(4): 145-152.

[14]VEDELH,HUANGXiangyu.ImpactofGroundBasedGPSDataonNumericalWeatherPrediction[J].JournaloftheMeteorologicalSocietyofJapan, 2004, 82(1B): 459-472.

[15]丁金才.GPS氣象學(xué)及其應(yīng)用[M]. 北京: 氣象出版社, 2009.

DINGJincai.GPSMeteorologyandItsApplications[M].Beijing:ChinaMeteorologicalPress, 2009.

[16]BEVISM,BUSINGERS,HERRINGTA,etal.GPSMeteorology:RemoteSensingofAtmosphericWaterVaporUsingtheGlobalPositioningSystem[J].JournalofGeophysicalResearch:Atmospheres, 1992, 97(D14): 15787-15801.

[17]劉春華. 三維投影矩陣的刻畫及迭代重建的加速研究[D]. 太原: 中北大學(xué), 2008.

LIUChunhua.TheStudyofBuildingThree-dimensionalProjectionMatrixandAccelerateIterativeImageReconstruction[D].Taiyuan:NorthUniversityofChina, 2008.

[18]BRESENHAMJE.AlgorithmforComputerControlofaDigitalPlotter[J].IBMSystemsJournal, 1965, 4(1): 25-30.

[19]AMANATIDESJ,WOOA.AFastVoxelTraversalAlgorithmforRayTracing[J].Eurographics, 1987, 87(3): 10.

(責(zé)任編輯：陳品馨)

修回日期： 2016-05-19

E-mail：metdingnan@163.com

Land-basedGPSWaterVaporTomographywithProjectionPlaneAlgorithm

DINGNan1,ZHANGShubi1, 2

1.SchoolofEnvironmentScienceandSpatialInformatics,ChinaUniversityofMiningandTechnology,Xuzhou221008,China; 2.KeyLaboratoryofResourcesandEnvironmentalInformationEngineeringofJiangsuProvince,ChinaUniversityofMiningandTechnology,Xuzhou221008,China

Watervaportomographyisofcrucialtheoreticalimportancetoandwidelyappliedinmeteorologyandclimatology,forinstance,extremeweatherforecastingandnowcasting,numericalweatherpredictionandclimatechangestudy.Conventionalwatervaportomographystrategyforcalculatingequationcoefficientneedsextensivecalculation.Tosolvethisproblem,projectionplanealgorithmisdevelopedforimprovingthespeedofcalculationandtheaccuracyofinversion.TheGPSsignaldatausedinthisarticlecomesfromHongKongSatRef.Theresult,whichiscomparedtothetraditionalalgorithm,indicatedthatnewalgorithmreducesthecomputingtimegreatlyandisnotaffectedbysizeoftomographyareaandthereisconsistencyintheresultsoftomographicinversionandradiosonde.

watervaportomography;signalpath;pixel;voxeltraversalalgorithm;algebraicreconstructiontechnique

DINGNan(1989—),male,PhDcandidate,majorsinGNSSwatervaportomography.

P228

A

1001-1595(2016)08-0895-09

國(guó)家自然科學(xué)基金(41504032); 江蘇省自然科學(xué)基金(BK20150175); 高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金(20130095110022)

2016-01-07

丁楠(1989—)，男，博士生，研究方向?yàn)镚NSS水汽層析。

引文格式：丁楠，張書畢.地基GPS水汽層析的投影面算法[J].測(cè)繪學(xué)報(bào)，2016,45(8)：895-903.DOI:10.11947/j.AGCS.2016.20160010.

DINGNan,ZHANGShubi.Land-basedGPSWaterVaporTomographywithProjectionPlaneAlgorithm[J].ActaGeodaeticaetCartographicaSinica,2016,45(8):895-903.DOI:10.11947/j.AGCS.2016.20160010.