一種改進(jìn)的電離層層析投影矩陣生成算法

劉大鵬，岳建平

(河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098)

一種改進(jìn)的電離層層析投影矩陣生成算法

劉大鵬，岳建平

(河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098)

正確建立投影矩陣是電離層層析的必要條件。由于反演區(qū)域較大且觀測數(shù)據(jù)繁多，同時考慮到離散后格網(wǎng)的復(fù)雜性，傳統(tǒng)的投影矩陣的計算算法中存在大量的重復(fù)計算，嚴(yán)重影響計算效率。文中對傳統(tǒng)算法流程進(jìn)行調(diào)整，在不改變計算結(jié)果的情況下，消除冗余計算，大大提高了計算效率。

電離層層析；投影矩陣；重復(fù)計算；計算效率

電離層是日地空間中的重要組成部分，它是由于地球大氣受紫外線輻射及X射線等的影響發(fā)生電離形成的。電離層中分布的電子會對穿過其中的GPS信號傳播路徑產(chǎn)生影響，從而產(chǎn)生電離層延遲誤差。傳統(tǒng)的電離層改正模型認(rèn)為電離層是一個距地面一定高度的無限薄的連續(xù)閉合曲面。然而這種模型只能反映電離層的水平結(jié)構(gòu)，對于垂直結(jié)構(gòu)及GPS信號在電離層內(nèi)部的具體情況無法反映[1-2]。基于GPS的電離層層析技術(shù)可以有效重構(gòu)電離層的空間三維結(jié)構(gòu)，同時了解和掌握電離層電子密度的三維分布狀態(tài)，有助于提高電離層延遲改正的精度[3-5]。

電離層層析模型一般分為函數(shù)基與像素基兩類。函數(shù)基模型通常利用一組函數(shù)來描述電離層中電子密度的空間分布；而像素基則是將待反演區(qū)域劃分為一系列的網(wǎng)格(像素)，然后在所選擇的參考框架和反演時段內(nèi)，假定每個網(wǎng)格內(nèi)的電子是均勻分布的，從而進(jìn)行電子密度反演[6-8]。

像素基電離層層析時，須求算投影矩陣，本文針對傳統(tǒng)投影矩陣求解算法中計算冗余的問題，提出“先整體，后局部”的思想，對算法進(jìn)行改進(jìn)，有效減少了多余計算，提高解算效率。

1 電離層層析原理

電離層層析就是根據(jù)反演區(qū)域內(nèi)一系列站星連線上的TEC值反演該區(qū)域內(nèi)電子密度在空間中的分布。TEC值就等于電離層中電子密度沿GPS信號傳播路徑的線積分[9]，即

(1)

(2)

式中:Aij為第i條GPS射線在第j個網(wǎng)格內(nèi)的截距，若某射線未穿過該網(wǎng)格，則相應(yīng)的截距定為0；xj為第j個網(wǎng)格中的電子密度；εi為第i條GPS射線的觀測噪聲與離散誤差。

2 改進(jìn)的投影矩陣生成算法

式(2)中系數(shù)矩陣Aij又稱為投影矩陣。通常計算Aij的方法為：計算射線與每個網(wǎng)格的6個包圍面的6個交點(diǎn)，識別并保留網(wǎng)格內(nèi)的交點(diǎn)，從而計算該射線在相應(yīng)網(wǎng)格內(nèi)的截距，進(jìn)而擴(kuò)展到整個反演區(qū)域，最終構(gòu)成投影矩陣。然而這種“先局部，后整體”的做法有如下2種缺點(diǎn)：①就單個網(wǎng)格而言，一條射線與網(wǎng)格的交點(diǎn)數(shù)通常為0個或2個(不考慮射線恰好通過該網(wǎng)格頂點(diǎn)的情況)，這樣每條射線與每個網(wǎng)格就進(jìn)行了至少4次多余計算。②對于整個待反演區(qū)域的眾多網(wǎng)格而言，相鄰的網(wǎng)格有公共面，這個公共面就會多次與GPS射線求交點(diǎn)，從而產(chǎn)生大量的重復(fù)計算。針對這2個問題，本文對傳統(tǒng)算法進(jìn)行調(diào)整，提出一種“先整體，后局部”的方法，可以較好地避免計算冗余的問題。所謂“先整體，后局部”，就是將在待反演區(qū)域內(nèi)定義一系列的曲面，計算每一條射線與每個曲面的交點(diǎn)坐標(biāo)，在同一條射線上，相鄰兩個交點(diǎn)所夾的線段長度就是該射線在某個網(wǎng)格內(nèi)的截距，只須識別該截距的網(wǎng)格歸屬，就可以確定它在投影矩陣中的位置，從而計算投影矩陣。該方法的優(yōu)點(diǎn)在于它讓每條射線與每個曲面只計算一次交點(diǎn)，從而減少計算量，提高解算效率。雖然加入了“識別”的步驟，但該環(huán)節(jié)計算量很小，通過合理的算法設(shè)計，幾乎不影響計算效率。算法流程如圖1所示。

圖1 改進(jìn)算法流程

新算法計算投影矩陣需要經(jīng)過以下7個步驟：

1)確定待反演區(qū)域范圍。確定經(jīng)度、緯度、高度范圍及在這3個方向上的間隔。經(jīng)緯度范圍根據(jù)需求確定，跨度通常為10°～20° ，考慮到電離層電子密度沿經(jīng)度方向的變化小于緯度方向，通常將網(wǎng)格間隔在經(jīng)緯度方向上分別定為1°～1.5°和0.5°；高度范圍通常設(shè)定為100～1 000 km，高度間隔取50～100 km。

2)根據(jù)精密星歷及測站坐標(biāo)，生成站星連線序列，并為每一條站星連線定義一個有效點(diǎn)序列，用于存儲射線與曲面的有效交點(diǎn)(即剛好落在測站和衛(wèi)星之間，且位于反演區(qū)域內(nèi)的交點(diǎn))。

3)在反演區(qū)域內(nèi)，分別沿經(jīng)度、緯度和高度方向生成3個曲面序列。由于所有數(shù)據(jù)均是通過GPS觀測獲得，故而計算都是在WGS-84坐標(biāo)系下進(jìn)行，因此，有必要進(jìn)行空間直角與大地坐標(biāo)間的轉(zhuǎn)換，從而得到測站與衛(wèi)星的大地坐標(biāo)值[10]。具體轉(zhuǎn)換公式為

(3)

式中:r為待反演區(qū)域平均地球半徑；H為某高度面的高度。設(shè)地面測站G的坐標(biāo)為(XG,YG,ZG)，衛(wèi)星S的坐標(biāo)為(XS,YS,ZS)，則這條GPS射線的方程為

(4)

要計算射線與某網(wǎng)格的截距，須先計算射線與每個曲面的交點(diǎn)。沿經(jīng)度、緯度和高度3個方向構(gòu)成的曲面類型各異，從而須分別計算。

①經(jīng)度方向：沿經(jīng)度方向構(gòu)成的曲面是一個過Z軸的平面，如圖2(a)所示，它與起始子午面夾角為經(jīng)度L，則該曲面方程為

tanL·x-y=0.

(5)

事先認(rèn)定射線GS為有效射線(GS與反演區(qū)域內(nèi)所有曲面的交點(diǎn)都在反演區(qū)域內(nèi))，若穿刺點(diǎn)c在點(diǎn)G和S之間，則判定c是有效點(diǎn)，并將其加入射線GS的有效點(diǎn)序列，否則就將其剔除。

②緯度方向：沿緯度方向構(gòu)成的曲面是一個頂點(diǎn)在原點(diǎn)，旋轉(zhuǎn)軸為z軸的圓錐面，如圖2(b)所示，它的錐半角與緯度B互余，則該曲面方程為

tan2B·(x2+y2)-z2=0.

(6)

事先認(rèn)定射線GS為有效射線，若GS所在直線與該圓錐面至少有一個交點(diǎn)，則在點(diǎn)G和S之間的為有效點(diǎn)，將其加入射線GS的有效點(diǎn)序列，并剔除另外一點(diǎn)。

③高度方向：在待反演區(qū)域內(nèi)，將距地面同一高度的點(diǎn)構(gòu)成的曲面看作球心在原點(diǎn)的球面，如圖2(c)所示，則該球面方程為

x2+y2+z2=(r+H)2.

(7)

式中:r為反演區(qū)域平均地球半徑；H為該高度面的高度。事先認(rèn)定射線GS為有效射線，根據(jù)實際情況可以判定，GS所在直線與此球面必有2個交點(diǎn)，選定在點(diǎn)G和S之間的一個為有效點(diǎn)，并將其加入射線GS的有效點(diǎn)序列，剔除另一點(diǎn)[11]。

圖2 穿刺點(diǎn)坐標(biāo)計算示意圖

4)排序。每條射線上的有效點(diǎn)序列只是按照求穿刺點(diǎn)的順序依次存入有效點(diǎn)，在空間分布上并沒有任何規(guī)則，這樣的順序不便于投影矩陣元素(截距)的計算，因此應(yīng)將有效點(diǎn)序列按一定的原則進(jìn)行重排。如圖3所示，選取1條射線與8個網(wǎng)格相交作為研究對象，共產(chǎn)生4個有效點(diǎn)A,B,C,D，可以發(fā)現(xiàn)，這4個點(diǎn)與測站G的距離是依次遞增的。因此，只要將一條射線的有效點(diǎn)序列中的點(diǎn)按照與測站G的距離由近及遠(yuǎn)依次排列，即可使序列中相鄰索引值的兩點(diǎn)之間的距離就是該射線在某網(wǎng)格內(nèi)的截距。

圖3 有效交點(diǎn)排序示意圖

5)網(wǎng)格編號。每個網(wǎng)格在三維空間中都有一個具體的位置，而在投影方程中的每一行中，網(wǎng)格是一維排布的。這就需要通過一定的數(shù)學(xué)方法，實現(xiàn)網(wǎng)格在空間中的三維位置與其在投影方程中的一維位置的互相轉(zhuǎn)換。本文首先鎖定一個高度面，同時鎖定一個緯度鏈，沿經(jīng)度開始排布，完成一條緯度鏈后，轉(zhuǎn)入下一個緯度，直到完成一個高度面的排序，然后轉(zhuǎn)入下一個高度面，直到處理完所有網(wǎng)格[12]，編號結(jié)果如圖4所示。

圖4 網(wǎng)格編號結(jié)果示意圖

6)識別。經(jīng)過3)、4)兩個步驟，每條射線在格網(wǎng)中的截距已經(jīng)確定，但還無法知曉每段截距對應(yīng)的網(wǎng)格索引值。設(shè)某網(wǎng)格的經(jīng)度范圍為(Lmin，Lmax)，緯度范圍為(Bmin，Bmax)，高度范圍為(Hmin，Hmax)，則判定一點(diǎn)是在該網(wǎng)格中的依據(jù)為

(8)

根據(jù)立體幾何知識可以知道，若一個線段在某個網(wǎng)格內(nèi)，那么它的中點(diǎn)也一定在該網(wǎng)格內(nèi)[13-14]。對于一條線段，計算出其中點(diǎn)的三維直角坐標(biāo)后，將其轉(zhuǎn)化為大地坐標(biāo)，即可判斷該線段所屬網(wǎng)格索引值。

圖5為網(wǎng)格三維位置索引與一維編號互相轉(zhuǎn)化的流程圖，LID、BID及HID分別表示某網(wǎng)格在經(jīng)度、緯度和高度方向上的索引值；Lnum、Bnum及Hnum表示反演區(qū)域在經(jīng)度、緯度及高度方向上劃分網(wǎng)格的個數(shù)。運(yùn)算符int表示取整。根據(jù)圖5，只要輸入相應(yīng)的參數(shù)，即可判斷線段所屬網(wǎng)格的三維空間位置及一維索引值，并實現(xiàn)互相轉(zhuǎn)化[15]。

圖5 網(wǎng)格三維位置索引與一維編號互相轉(zhuǎn)化流程圖

7)建立投影矩陣。根據(jù)投影矩陣的定義，將截距按照投影矩陣的定義進(jìn)行排列，構(gòu)成投影矩陣。

3 數(shù)值模擬實驗

為驗證改進(jìn)算法在計算投影矩陣中的可行性及算法的高效性，本文設(shè)計了數(shù)值模擬實驗。選定實驗區(qū)域為東經(jīng)100°～110°，間隔為1°；緯度范圍為北緯20°～30°，緯度間隔取0.5°；高度范圍取100～1 000 km，高度間隔取50 km。從而劃定了3 600個網(wǎng)格。假設(shè)在實驗區(qū)內(nèi)架設(shè)了23臺接收機(jī)，其分布如圖6所示。在反演區(qū)域上空20 000±200 km范圍隨機(jī)設(shè)置8個點(diǎn)來模擬衛(wèi)星的位置。假定所有測站都能接收到8個衛(wèi)星信號，則在任意時刻，總有184條射線穿過反演區(qū)域。分別用傳統(tǒng)算法和改進(jìn)算法建立投影矩陣。計算結(jié)果及統(tǒng)計量如表1所示。

圖6 實驗區(qū)域測站分布示意圖

表1 傳統(tǒng)算法與改進(jìn)算法計算結(jié)果及計算量統(tǒng)計

表1中絕對誤差E的定義為

(9)

式中，a0,a1分別表示兩種算法生成的投影矩陣中的元素。

由于所有的測站與衛(wèi)星都在反演區(qū)域內(nèi)，無需考慮兩種算法在計算時點(diǎn)位判斷與刪除的差異。由表1可知，通過改進(jìn)算法計算投影矩陣的結(jié)果與改進(jìn)前完全相同，而計算量卻減少了80%以上，計算效率顯著提高。

4 結(jié) 論

投影矩陣計算中的冗余計算嚴(yán)重影響了計算效率，本文通過調(diào)整傳統(tǒng)算法中的計算順序，提出的“先整體，后局部”的計算思想，使每條射線與每個經(jīng)度面、緯度面及高度面只進(jìn)行一次交點(diǎn)計算，消除了多余計算。并通過模擬實驗驗證了改進(jìn)后算法的有效性、可靠性和優(yōu)越性。從計算結(jié)果來看，計算結(jié)果與改進(jìn)前的算法完全一致；就計算效率而言，由于計算量的降低，計算效率得到一倍以上的提升。

[1] 吳寒，姚宜斌，陳鵬，等.GNSS電離層層析成像算法[J].武漢大學(xué)學(xué)報(信息科學(xué)版)，2013，28(12)：1405-1408.

[2] 肖宏波，門守強(qiáng).電離層層析圖像重建的改進(jìn)算法[J].西安工業(yè)大學(xué)學(xué)報，2014，34(5)：345-348.

[3] 席超，蔡成林，李思敏，等.一種基于三角分區(qū)的廣域電離層改正新方法[J].武漢大學(xué)學(xué)報(信息科學(xué)版)，2015，40(3)：390-394.

[4] 章洪平.基于地基GPS的中國區(qū)域電離層監(jiān)測與延遲改正研究[D].上海：中國科學(xué)院上海天文臺，2006:1-100.

[5] 周忠謨，易杰軍，周琪.GPS衛(wèi)星測量原理與應(yīng)用[M].北京：測繪出版社，2004:1-200.

[6] 聞德保.基于GNSS的電離層層析算法及其應(yīng)用[M].北京：測繪出版社，2013：7-8.

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

[8] 金雙根，朱文耀.用地基GPS觀測資料映射三維電離層剖面[J].大地測量與地球動力學(xué)，2007,27(5)：49-53.

[9] 蘭孝奇，李森，解坤.利用區(qū)域地基GPS的電離層層析成像研究[J].測繪學(xué)報，2012，37(4)：17-18.

[10] 楊帆.基于GPS的區(qū)域電離層層析算法及其應(yīng)用研究[D].沈陽：東北大學(xué)，2013：25-41.

[11] 陳必焰.電離層層析技術(shù)及其應(yīng)用[D].長沙：中南大學(xué)，2012：18-20.

[12] 聞德保，張嘯，張光勝，等.基于選權(quán)擬合法的電離層電子密度層析重構(gòu)[J].地球物理學(xué)報，2014,57(8)：2395-2403.

[13] 吳熊斌，徐繼生，馬淑英.一種改進(jìn)的電離層層析成像算法[J].地球物理學(xué)報，2000，43(1)：19-28.

[14] 孫曉安，陳淑珍，徐繼生.電離層斷層成像數(shù)值模擬[J].武漢大學(xué)學(xué)報(理學(xué)版)，1994(6)：57-64.

[15] 肖宏波，史小紅，薛昆.電離層層析成像技術(shù)與電離層反演[J].現(xiàn)代電子技術(shù)，2009(18)：131-133.

[責(zé)任編輯：劉文霞]

An improved generation algorithm for projection matrix in ionospheric tomography

LIU Dapeng,YUE Jianping

(School of Earth Science and Engineering,Hohai University,Nanjing 210098,China)

The establishment of projection matrix is an essential requirement in the ionospheric tomography.Because of the large inversion region and huge amount of observation data,and of the complex grid,the traditional algorithm leads to a lot of repetitive calculations,which seriously affects the computational efficiency.In this paper,the traditional algorithm flow is adjusted.The redundant computation can be eliminated and the computational efficiency be greatly improved without changing the calculation result.

ionospheric tomography;projection matrix;repetitive calculation;computational efficiency

2016-07-11

國家自然科學(xué)基金資助項目(41174002)

劉大鵬(1993-)，男，碩士研究生.

著錄：劉大鵬，岳建平.一種改進(jìn)的電離層層析投影矩陣生成算法[J].測繪工程，2017，26(9)：51-55.

10.19349/j.cnki.issn1006-7949.2017.09.011

P228

A

1006-7949(2017)09-0051-05