劉大鵬，岳建平

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

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

劉大鵬，岳建平

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

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

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

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

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

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

1 電離層層析原理

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

(1)

(2)

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

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

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

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

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

1)確定待反演區(qū)域范圍。確定經(jīng)度、緯度、高度范圍及在這3個(gè)方向上的間隔。經(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ù)精密星歷及測(cè)站坐標(biāo)，生成站星連線序列，并為每一條站星連線定義一個(gè)有效點(diǎn)序列，用于存儲(chǔ)射線與曲面的有效交點(diǎn)(即剛好落在測(cè)站和衛(wèi)星之間，且位于反演區(qū)域內(nèi)的交點(diǎn))。

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

(3)

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

(4)

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

①經(jīng)度方向：沿經(jīng)度方向構(gòu)成的曲面是一個(gè)過(guò)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)成的曲面是一個(gè)頂點(diǎn)在原點(diǎn)，旋轉(zhuǎn)軸為z軸的圓錐面，如圖2(b)所示，它的錐半角與緯度B互余，則該曲面方程為

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

(6)

事先認(rèn)定射線GS為有效射線，若GS所在直線與該圓錐面至少有一個(gè)交點(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ù)實(shí)際情況可以判定，GS所在直線與此球面必有2個(gè)交點(diǎn)，選定在點(diǎn)G和S之間的一個(gè)為有效點(diǎn)，并將其加入射線GS的有效點(diǎn)序列，剔除另一點(diǎn)[11]。

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

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

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

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

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

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

(8)

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

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

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

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

3 數(shù)值模擬實(shí)驗(yàn)

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

圖6 實(shí)驗(yàn)區(qū)域測(cè)站分布示意圖

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

表1中絕對(duì)誤差E的定義為

(9)

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

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

4 結(jié) 論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[14] 孫曉安，陳淑珍，徐繼生.電離層斷層成像數(shù)值模擬[J].武漢大學(xué)學(xué)報(bào)(理學(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

國(guó)家自然科學(xué)基金資助項(xiàng)目(41174002)

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

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

10.19349/j.cnki.issn1006-7949.2017.09.011

P228

A

1006-7949(2017)09-0051-05