快速網格路徑算法在對流層層析中的應用

楊 飛，郭際明，史俊波

（1.武漢大學 測繪學院，湖北 武漢 430079）

快速網格路徑算法在對流層層析中的應用

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

（1.武漢大學 測繪學院，湖北 武漢 430079）

在GNSS對流層層析應用中，常規的網格輻射路徑算法是利用空間平面與空間直線求截距。本文介紹了一種快速求解算法，用于計算網格輻射路徑，構造層析系數矩陣。結果表明，快速求解方法與常規方法解算精度相當，當射線條數大于500時，計算效率明顯提高。

對流層層析；層析網格；系數矩陣；網格輻射路徑

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

1 三維層析網格系數矩陣常規方法

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

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

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

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

圖1 衛星信號射線穿過層析網格示意圖

2 層析網格輻射路徑快速求解算法

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

1）點線面的初始化。將地面測站點A與衛星位置點B的信號簡化為射線AB，利用式（3）參數化，參數α在A點時為0，在B點時為1；將層析區域劃分為（Nx-1，Ny-1，Nz-1）個網格單元，利用式（4）將經度、緯度和高程面參數化：

式中，X1、Y1、Z1表示點A的三維坐標；Xp（i）、Yp（i）、Zp（i）表示網格化后的各個經度面、緯度面、高程面；dx、dy、dz分別表示經度方向、緯度方向和高程方向相鄰兩個平面間的距離。

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

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

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

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

3）利用上述參數集合求取射線與三維網格的交點距離Sqi,j,k。設m和m-1為射線與一個網格的兩個交點，則存在如下關系：

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

5）經過4）的計算，得到衛星信號穿過層析網格的位置（i，j，k）及長度構成層析方程系數矩陣[10]。

3 實驗驗證

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

圖2 地面層析網格劃分及測站位置圖

3.1 精度驗證

本次層析實驗中，測站WHHN和WHXZ的500 條衛星信號分別穿過網格單元的次數為10 520次和10 451次。其中，一條衛星信號最多穿過24個網格單元，最少穿過20個網格單元。以常規方法解算結果作為參考，計算快速解算方法中各測站不同數量衛星信號射線穿過層析網格的平均誤差和RMS，如圖3所示。

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

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

3.2 快速求解算法的效率驗證

為了驗證快速算法的效率，選取不同數量的射線進行層析網格求解，統計解算時間，每種情況進行10 次計算取平均時間，結果如圖4所示。

圖4 兩種方法解算時間對比圖

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

實驗結果表明，在GNSS對流層層析網格解算過程中，無論測站位置和衛星信號數量，尤其在衛星信號數量大時，快速解算法能提高計算效率。

4 結 語

三維層析網格系數矩陣是GNSS對流層層析解算的關鍵。對比常規方法和快速求解方法，2種方法的精度相當。當射線條數增多時，快速求解方法的解算效率明顯提高。

[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]．自然資源學報,2010,25(10):1 786-1 796

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

[4] 胡金玉．基于CORS參考站對區域水汽的研究[D]．成都:西南交通大學,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]．武漢:武漢大學,2013

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

[8] 宋淑麗．地基GPS網對水汽三維分布的監測及其在氣象學中的應用[D]．上海:中國科學院上海天文臺,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

楊飛，碩士，研究方向為區域對流層建模。

2015-06-30。

項目來源：國家自然科學基金資助項目（41474004）。