快速網(wǎng)格路徑算法在對流層層析中的應(yīng)用

楊 飛，郭際明，史俊波

（1.武漢大學(xué) 測繪學(xué)院，湖北 武漢 430079）

快速網(wǎng)格路徑算法在對流層層析中的應(yīng)用

楊 飛1，郭際明1，史俊波1

（1.武漢大學(xué) 測繪學(xué)院，湖北 武漢 430079）

在GNSS對流層層析應(yīng)用中，常規(guī)的網(wǎng)格輻射路徑算法是利用空間平面與空間直線求截距。本文介紹了一種快速求解算法，用于計(jì)算網(wǎng)格輻射路徑，構(gòu)造層析系數(shù)矩陣。結(jié)果表明，快速求解方法與常規(guī)方法解算精度相當(dāng)，當(dāng)射線條數(shù)大于500時(shí)，計(jì)算效率明顯提高。

對流層層析；層析網(wǎng)格；系數(shù)矩陣；網(wǎng)格輻射路徑

GNSS對流層層析是求解對流層水汽密度，從而反映對流層水汽在三維空間分布情況的方法[1]。對流層層析首先要?jiǎng)澐謱游鼍W(wǎng)格，建立層析觀測方程[2]，求解層析方程系數(shù)矩陣，該系數(shù)是每條GNSS射線穿過的層析網(wǎng)格內(nèi)的截距，精確的對流層層析網(wǎng)格輻射路徑算法是層析計(jì)算的關(guān)鍵。陳宏斌[3]提出了利用空間直線與空間平面的關(guān)系求解算法；胡金玉[4]提出了利用衛(wèi)星信號路徑與緯度面、經(jīng)度面、高度面聯(lián)合求解算法。這些都屬于常規(guī)的逐個(gè)平面解算方法，沒有詳盡闡述解算過程。在放射性治療應(yīng)用中，Siddon[5]提出了一種三維CT數(shù)組輻射路徑的快速算法。本文將醫(yī)學(xué)CT中的算法用于GNSS對流層層析網(wǎng)格輻射路徑的解算，利用已知測站坐標(biāo)和衛(wèi)星坐標(biāo)求解各段路徑在總路徑中的比例關(guān)系，可以得到與常規(guī)方法相同精度的結(jié)果，提高解算效率。

1 三維層析網(wǎng)格系數(shù)矩陣常規(guī)方法

層析觀測方程[6]如式（1）所示：

式中，nl、nm、nh分別為經(jīng)度、緯度和高程方向劃分的網(wǎng)格數(shù)量；Sqi,j,k為第q條衛(wèi)星信號穿過（i，j，k）網(wǎng)格單元的路徑長度；ρi,j,k為相應(yīng)立體網(wǎng)格的水汽密度；SWDq為第q條衛(wèi)星信號傾斜路徑延遲[7]。將式（2）寫成矩陣形式：

式中，SWD是n×l階傾斜路徑濕延遲向量，x為m×l階未知參數(shù)向量。A為n×m階層析系數(shù)矩陣，其第j行第i列的元素表示第j條射線穿過網(wǎng)格i的長度[8]。所以構(gòu)建層析觀測方程解算層析結(jié)果，需要進(jìn)行層析網(wǎng)格輻射路徑長度的求解。

對流層層析常采用逐個(gè)平面求解的方法。如圖1所示，將衛(wèi)星信號射線AB視為空間直線，將經(jīng)度、緯度、高程面視為空間平面，依次求空間直線AB與各個(gè)空間平面交點(diǎn)P1、P2、P3，分別得到線段P1P2、P2P3的長度及其所在網(wǎng)格的位置。

圖1 衛(wèi)星信號射線穿過層析網(wǎng)格示意圖

2 層析網(wǎng)格輻射路徑快速求解算法

快速求解算法將第一個(gè)交點(diǎn)P1作為起點(diǎn)，通過遞歸循環(huán)得到其他所有交點(diǎn)，并求出每個(gè)格網(wǎng)路徑長度具體步驟如下：

1）點(diǎn)線面的初始化。將地面測站點(diǎn)A與衛(wèi)星位置點(diǎn)B的信號簡化為射線AB，利用式（3）參數(shù)化，參數(shù)α在A點(diǎn)時(shí)為0，在B點(diǎn)時(shí)為1；將層析區(qū)域劃分為（Nx-1，Ny-1，Nz-1）個(gè)網(wǎng)格單元，利用式（4）將經(jīng)度、緯度和高程面參數(shù)化：

式中，X1、Y1、Z1表示點(diǎn)A的三維坐標(biāo)；Xp（i）、Yp（i）、Zp（i）表示網(wǎng)格化后的各個(gè)經(jīng)度面、緯度面、高程面；dx、dy、dz分別表示經(jīng)度方向、緯度方向和高程方向相鄰兩個(gè)平面間的距離。

2）求取參數(shù)集合{α}，即各個(gè)穿刺點(diǎn)在整條衛(wèi)星信號中的比例因子。首先，利用上述參數(shù)化點(diǎn)面計(jì)算類似求取αy（i）、αz（i）；接著，利用式（5）由上述參數(shù)值求取αmin和αmax：

然后，利用式（6）求解參數(shù)集合{αx}、{αy}、{αz}：

其中需要利用式（7）得到（imin，imax）、（jmin，jmax）、（kmin，kmax）：

利用類似公式求取相應(yīng)y、z和j、k參數(shù)。最后，利用式（8）得到參數(shù)集合{α}：

3）利用上述參數(shù)集合求取射線與三維網(wǎng)格的交點(diǎn)距離Sqi,j,k。設(shè)m和m-1為射線與一個(gè)網(wǎng)格的兩個(gè)交點(diǎn)，則存在如下關(guān)系：

4）確定網(wǎng)格[i（m），j（m），k（m）]在全體三維網(wǎng)格中的位置，與網(wǎng)格交點(diǎn)m和m-1有下列關(guān)系式：

5）經(jīng)過4）的計(jì)算，得到衛(wèi)星信號穿過層析網(wǎng)格的位置（i，j，k）及長度構(gòu)成層析方程系數(shù)矩陣[10]。

3 實(shí)驗(yàn)驗(yàn)證

為驗(yàn)證快速算法的準(zhǔn)確性和效率，以武漢CORS網(wǎng)絡(luò)進(jìn)行層析建模。經(jīng)度范圍是113.8°～115.0°，緯度范圍是30.0°～31.0°，高度范圍是0～10 km。首先，選取測站W(wǎng)HHN和WHXZ于 2013-05-06的GPS觀測數(shù)據(jù)，計(jì)算衛(wèi)星相對于測站的高度角和方位角。然后，劃分層析網(wǎng)格，水平分辨率設(shè)為0．2°，垂直分辨率設(shè)為500 m，共6×5×20個(gè)網(wǎng)格，如圖2所示。最后，分別選取兩個(gè)測站的500條衛(wèi)星信號進(jìn)行層析網(wǎng)格輻射路徑求解，比較2種方法的解算精度。

圖2 地面層析網(wǎng)格劃分及測站位置圖

3.1 精度驗(yàn)證

本次層析實(shí)驗(yàn)中，測站W(wǎng)HHN和WHXZ的500 條衛(wèi)星信號分別穿過網(wǎng)格單元的次數(shù)為10 520次和10 451次。其中，一條衛(wèi)星信號最多穿過24個(gè)網(wǎng)格單元，最少穿過20個(gè)網(wǎng)格單元。以常規(guī)方法解算結(jié)果作為參考，計(jì)算快速解算方法中各測站不同數(shù)量衛(wèi)星信號射線穿過層析網(wǎng)格的平均誤差和RMS，如圖3所示。

圖3 快速解算法與常規(guī)解算方法的平均誤差和RMS

快速求解算法的平均誤差隨著衛(wèi)星信號射線數(shù)量的增加而趨于穩(wěn)定，測站W(wǎng)HXZ穩(wěn)定于0．248 m（圖3a），測站W(wǎng)HHN穩(wěn)定于0．232 m（圖3b）。另一方面，兩個(gè)測站的RMS也隨著衛(wèi)星信號射線數(shù)量的增加而趨于穩(wěn)定，其穩(wěn)定值分別為0．417 m（圖3c）和0．371 m（圖3d）。實(shí)驗(yàn)結(jié)果說明，快速解算方法與常規(guī)解算方法結(jié)果差值的平均誤差在25 cm左右，RMS在40 cm左右。相對于格網(wǎng)長度（≥500 m），平均誤差比例小于等于0．05%， RMS不超過0．08%，可以用于構(gòu)建層析系數(shù)矩陣。

3.2 快速求解算法的效率驗(yàn)證

為了驗(yàn)證快速算法的效率，選取不同數(shù)量的射線進(jìn)行層析網(wǎng)格求解，統(tǒng)計(jì)解算時(shí)間，每種情況進(jìn)行10 次計(jì)算取平均時(shí)間，結(jié)果如圖4所示。

圖4 兩種方法解算時(shí)間對比圖

從圖4可以看出，無論衛(wèi)星信號射線條數(shù)的多少（0～15 000），快速解算方法所用時(shí)間都小于常規(guī)解算方法。當(dāng)射線條數(shù)小于500時(shí)，時(shí)間效率的差別不明顯。以WHHN測站為例，當(dāng)射線條數(shù)為500時(shí)，兩種方法解算時(shí)間都小于2 s，分別為1．545 s和0．125 s，快速解算法比常規(guī)解算法快1．420 s。隨著射線條數(shù)的逐漸增加，兩種方法的效率差別更加明顯；當(dāng)射線為1 000條時(shí)，快速解算方法比常規(guī)解算法快2．753 s；當(dāng)射線為15 000條時(shí)，快速解算方法比常規(guī)解算法快39．518 s。

實(shí)驗(yàn)結(jié)果表明，在GNSS對流層層析網(wǎng)格解算過程中，無論測站位置和衛(wèi)星信號數(shù)量，尤其在衛(wèi)星信號數(shù)量大時(shí)，快速解算法能提高計(jì)算效率。

4 結(jié) 語

三維層析網(wǎng)格系數(shù)矩陣是GNSS對流層層析解算的關(guān)鍵。對比常規(guī)方法和快速求解方法，2種方法的精度相當(dāng)。當(dāng)射線條數(shù)增多時(shí)，快速求解方法的解算效率明顯提高。

[1] Flores A,Ruffini G,Rius A．4D Tropospheric Tomography Using GPS Slant Wet Delays[J]．Annales Geophysicae, 2000, 18(2):223-234

[2] 曹玉靜,劉晶淼,梁宏,等．基于地基GPS層析大氣水汽資源的方法研究[J]．自然資源學(xué)報(bào),2010,25(10):1 786-1 796

[3] 陳宏斌．地基GPS層析水汽三維分布研究[D]．成都:西南交通大學(xué),2014

[4] 胡金玉．基于CORS參考站對區(qū)域水汽的研究[D]．成都:西南交通大學(xué),2014

[5] Siddon R．Fast Calculation of the Exact Radiological Path for a Three-Dimensional CT Array[J]． Medical Physics,1985,12(2):252-255

[6] 范士杰． GPS海洋水汽信息反演及三維層析研究[D]．武漢:武漢大學(xué),2013

[7] 王曉英．地基GNSS層析對流層水汽若干關(guān)鍵技術(shù)研究[D]．南京:南京信息工程大學(xué),2013

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

[9] Censor Y． Finite Series-Expansion Reconstruction Methods[C]．Proceedings of the IEEE,1983, 71(3):409-419

[10] Christiaens M,Sutter B D,Bosschere K D,et al． A Fast, Cache-Aware Algorithm for the Calculation of Radiological Paths Exploiting Subword Parallelism[J]．Journal of Systems Architecture,1999, 45(10):781-790

P228.4

B

1672-4623（2016）06-0045-03

10.3969/j.issn.1672-4623.2016.06.015

楊飛，碩士，研究方向?yàn)閰^(qū)域?qū)α鲗咏！?/p>

2015-06-30。

項(xiàng)目來源：國家自然科學(xué)基金資助項(xiàng)目（41474004）。