一種聯合地基GNSS 和測高儀數據的電離層層析成像新算法

歐明 吳家燕 陳龍江 甄衛民 呂夢海

（1.中國電波傳播研究所，青島 266107；2.中電科(青島)電波技術有限公司，青島 266107）

引言

電離層是指地球上空60～1 000 km 高度大氣電離區域，穿越該區域的無線電波會產生折射、反射、散射和吸收等效應，從而影響導航、通信、測控、遙感等諸多無線電信息系統的性能[1].利用各類技術手段來獲取電離層的特征參量，研究電離層中的各種現象，從而揭示其背后的物理機制，具有重要的科學乃至政治、軍事和經濟意義[2].電子密度是表征電離層狀態最為直接的特征參量之一，獲取精確的電子密度變化即可量化電離層對衛星導航、通信、雷達的影響[3-4].

電離層層析成像(computerized ionospheric tomography,CIT)是隨著計算機硬件性能的提升和衛星電離層探測的進步而逐漸發展起來的一種非常重要的電離層遙感工具.CIT 是利用衛星測量的電離層積分總電子含量(total electron content,TEC)反演得到電子密度的過程[5].目前CIT 主要數據來源包括低軌衛星信標、中高軌GNSS 數據，而早期CIT 數據來源主要是低軌信標[6].近年來基于地基GNSS 臺站進行CIT 是主要的發展方向，國內外研究人員利用GNSS 三維CIT 技術，獲取了大量電離層水平、垂直方向的結構信息[7].

受限于掃描視角和臺站數目非均勻分布的影響，基于地基GNSS 數據的CIT 算法不同程度存在著電子密度成像垂直分辨率不高的問題[8-9].近年來，多手段聯合CIT 技術開始成為世界各國的研究熱點.Hajj 等提出了聯合地基和天基觀測反演電離層電子密度(ionospheric electron density,IED)的設想[10]；Rius等嘗試了地基GPS 數據和MET 衛星掩星數據的CIT 試驗[11]；Dear 等提出采用地面測高儀優化GPS 數據實現IED CIT[12]；英國Bath 大學開發了多手段數據分析軟件(multi-instrument data analysis system,MIDAS)，利用模式化基函數(球諧函數與經驗正交函數)進行電離層表征，同時利用奇異值分解(singular value decomposition,SVD)的方法對表征函數的權重進行求解，以降低CIT 反演過程中由于數據稀疏性造成的矩陣求逆困難[13-14].為了提高三維CIT 的穩定性，Ma 等提出利用測高儀聯合GPS 觀測，結合神經網絡的方法獲取三維電子密度[15]；Yao 等提出了組合型CIT 算法，以解決單獨使用函數基CIT 模型和像素基CIT 模型存在的階數選擇困難、反演參數多、計算效率低的問題[16].肖銳選擇經驗正交函數(empirical orthogonal function,EOF)作為三維時變CIT 算法的基函數，從而大大降低了三維CIT 需要求解的未知數數量，提升了算法的穩定性和時效性[17]；聞德保針對GPS 基CIT 中的不適定問題，先后提出了多種解決方案，如選權擬合法等[18]；Li 等提出了一種新CIT 方法，旨在通過使用受約束的最小二乘來緩解CIT 矩陣的不適定問題，以提高IED 的反演精度[19]；Ssessanga N 等通過將測高儀數據結合同化技術中，以提高CIT 的穩定性并避免電子密度反演過程中的負值問題[20].由此可以看出，聯合多種手段觀測數據已成為提升CIT 分辨率的主要方法.

本文針對目前地基GNSS CIT 存在垂直分辨率較低和算法穩定性不足的問題，提出了一種聯合地基GNSS 和測高儀數據的三維CIT 新算法，該算法綜合了GNSS 電離層探測水平分辨率較高和測高儀電離層垂直探測分辨率較高的優點，以全球電離層無線電觀測網(global ionospheric radio observatory,GIRO)測高儀數據驅動更新IRI 模型，將驅動后的IRI 模型作為背景電離層模型，再利用改進的代數重建迭代(algebraic reconstruction technique,ART)算法進行CIT，有效提升了電離層三維電子密度重構的精度，相比僅利用地基GNSS 數據，該算法電離層hmF2的反演精度明顯提升，通過實測數據對本文算法的精度進行了驗證.

1 CIT 算法

1.1 CIT 原理

三維CIT 主要是指通過一系列衛星和地面接收機間無線電信號傳播路徑上的積分TEC 測量來重構區域內三維的IED 分布.地面GNSS 接收機所獲得的TEC 可以表示為沿信號傳播路徑上電子密度的積分，即：

式中：di為總電子含量；i為射線編號；Ne(r)為沿傳播路徑隨空間r變化的電子密度；S為地面接收機至衛星的視線路徑.

將反演區域按照經度、緯度和高度方向劃分為N個網格，采用級數展開型算法，利用一組空間基函數hk(r)將電子密度分布離散化，有：

式中：K表示展開所取的級數；ak為權重系數.由于ak不隨時間變化，有：

式中：Δsin表示第i個射線在第n個網格內的截距；gik由 權重系數ak與基函數h k(r)相乘得到.

在給定基函數的具體形式后，根據相應射線的幾何位置可以確定gik.基函數包括像素基函數和函數基函數兩大類.本文采用像素基函數，其將空間離散為有限數目的網格，且將網格中的電子密度當作常數分布.若網格內有無線電信號傳播路徑通過，則基函數設定為1，否則為0，即：

因此，電離層問題可以轉換為求解下列線性方程組的問題：

式中：向量d由TEC 觀測數據di組成；采用像素基函數后，待求量X即為電離層電子密度Ne；投影矩陣A由基函數加權的第i條信號傳播路徑在第j個網格內投影相對于參考路徑的增量gik組成.采用CIT 算法求得基函數權重，便可得到電子密度的分布[21].

1.2 改進的ART 算法

CIT 方程組，即式(5)的求解，目前常用的算法為ART 算法.ART 算法對初始估計反復進行迭代修正，直到某種設定的條件得到滿足為止，每一次修正針對一個方程進行.這類算法比較節省計算機內存資源且易于操作，應用較為廣泛[6].算法的計算流程如下：

1)首先設定電離層模型背景場初值X(0)，且

2)基于實測的電離層數據，背景場初值按照以下迭代方式求解 ：

式中：j為網格編號；分別為迭代k+1 次和k次的第j個網格的電子密度；λ為松弛因子；aij為第i條路徑在第j個網格內的投影；‖ai‖為第i條路徑總長.

由式(7)可以看出，ART 算法迭代過程中，方程的迭代順序、迭代輪次、松弛因子的具體取值均會對CIT 的最終結果造成一定的影響.松弛因子與ART 算法的收斂速度直接相關，其取值影響收斂時的震蕩幅度和收斂所需的迭代次數.

傳統ART 算法松弛因子取一個固定常數，迭代順序如何給定等方面存在較大的隨意性，使得CIT 結果存在穩定性較差的問題.為消除算法的不穩定性，本文對ART 算法進行相應的改進，具體實現方法如下：

1) λk采用與迭代輪次k相關的函數，使得迭代松弛因子能夠隨著迭代輪次的增加而逐步降低，從而保證迭代收斂的穩定性，設定方法如下：

2)設置迭代截止門限 tol，當 tol＜10-8時停止迭代，門限值計算方法為

3)改進ART 算法的迭代順序，算法不再按照i=mod(k,m)+1的順序進行迭代，而是采用隨機數的方法確定迭代順序，以消除傳統ART 算法采用固定迭代次序對電子密度反演精度的影響.

4)網格內同一高度層的電子密度加入水平約束，以實現對無射線穿越網格內電子密度的調整.水平約束采用雙線性插值算法實現：

式中：w表示雙線性插值所用的加權系數，其實現方法可參考文獻[22].

1.3 測高儀數據驅動的電離層背景場初值

背景場初值對三維CIT 的結果影響非常大，給出一個更為合理的背景初值將在一定程度上提高CIT 的精度.傳統的CIT 方法一般選擇經驗電離層模型輸出的電子密度值作為背景場，其中最為常用的模型即國際參考電離層模型IRI[23].由于地基GNSS用于三維CIT 的垂直分辨率較低，利用測高儀數據數據驅動IRI 可在高度方向上提升地基GNSS 三維CIT 的分辨率[24-25].

IRI 模型基于CCIR/URSI 系數計算全球電離層F2層特征參數foF2和hmF2.本文利用GIRO 發布的測高儀自動判讀數據，通過分步線性最優估計方法更新CCIR/URSI 系數[26]，使得IRI 模型計算的電離層foF2和hmF2與測高儀觀測數據之間的方差最小，進而獲得更準確的全球電子密度分布作為CIT 的背景初值.分步線性最優估計具體的實現方法可參考文獻[26].

2 數據來源

2.1 GNSS 數據

選擇中國及周邊區域約290 個地基GNSS 臺站的數據參與CIT，數據由IGS 和中國陸態網獲得，臺站的地理位置分布如圖1 所示，其中紅色實心方塊為地基GNSS 臺站.

圖1 CIT 使用的地基GNSS 臺站(紅色方塊)和精度評估用的測高儀臺站(藍色三角)分布Fig.1 Distributions of ground-based GNSS stations (red□) for CIT and ionosonde stations (blueΔ) for accuracy validation

在進行CIT 反演之前，首要的任務是利用地基GNSS 觀測文件獲取電離層TEC 信息.在計算TEC前，首先利用Bernese 軟件對RINEX 文件進行數據預處理.Bernese 軟件是瑞士伯爾尼大學天文研究所針對GNSS 高精度測繪研制的一款經典軟件，該軟件廣泛應用于GPS 數據后處理領域中[27].本文首先利用Bernese 軟件對RINEX 格式文件中GNSS的觀測數據進行周跳和異常值的探測，并采用雙頻的載波相位觀測數據平滑碼偽距，最后平滑處理后的觀測數據仍以RINEX 格式的文件輸出.這樣處理可較好剔除GNSS 數據中的部分異常觀測值，從而提升CIT 的質量.

利用意大利Gran Sasso 太空研究所研制的GNSS TEC 數據處理軟件對各臺站的RINEX 觀測數據進行處理，從而計算得到電離層TEC 信息[28].GNSS_2016 軟件通過讀取GNSS 觀測數據及相應GNSS 衛星的星歷，可輸出衛星編號、觀測時刻、衛星仰角、方位角、穿刺點經緯度及傾斜路徑上的TEC 值.

2.2 Madrigal TEC 數據

美國麻省理工大學Madrigal 能夠管理和提供各種形式的存檔的和實時的地球空間物理在線數據，是一個功能強大的數據庫.Madrigal 采用統計框架對全球數千個地基GNSS 數據進行處理，用于估計GNSS 非電離層差分延遲偏差，實現了1°×1°×5 min 分辨率電離層TEC 的高精度估計，在空間物理研究領域應用非常廣泛[29].本文選擇Madrigal TEC 數據產品對CIT 重構的TEC 地圖進行精度評估.

2.3 電離層測高儀數據

測高儀數據下載自全球電離層無線電觀測站GIRO.GIRO 是美國馬薩諸塞州Lowell 大學發起和建立的全球電離層數字測高儀站網.截止目前，GIRO的臺站已經覆蓋了全球27 個國家并在美國Lowell建立了一個數據處理中心LGDC，由LGDC 負責全球超過50 個測高儀臺站的數據處理(電離圖自動判讀與數據質量控制)、建模(IRI 實時同化模型IRTAM 和全球底部電離層時效同化模型GAMBIT)，以及DIDBase 等數據產品的發布[30].電離圖自動判讀數據即為DIDBase 數據之一.在CIT 研究中利用GIRO 自動判讀的測高儀數據對IRI 模型進行驅動，GIRO 包括了全球50 個臺站的數據，其中有3 個國內臺站和47 個國外臺站，具體臺站列表請見http://giro.uml.edu/.

同時，利用獨立的中國國防電波觀測站網19 個測高儀臺站數據對CIT 的精度進行評估.圖1 中藍色三角形即為精度評估所用的測高儀臺站的位置分布.需要說明的是，在本文電子密度精度評估時，采用全部19 個臺站的數據進行foF2評估.但是考慮到目前電波觀測網的測高儀的判讀數據不包含hmF2，需要利用POLAN 程序通過人工判讀電離圖的方式處理得到，因此本文僅選擇北京、蘇州和海口3 個不同緯度臺站的hmF2數據進行精度評估.

3 CIT 結果分析

3.1 電離層TEC 地圖精度分析

CIT 算法網格劃分，設定為水平方向經度和緯度上的網格大小為2.5°×2.5°，垂直方向450 km 以下高度網格間距為10 km，450～1 000 km 高度區間網格間距為25 km，時間分辨率設定為15 min.首先分析CIT 算法給出的電離層TEC 地圖的重構精度.為顯示方便，圖2～5 分別按照1 h 時間間隔，給出了2017-09-06T00:00—23:00UT 中國及周邊區域Madrigal、歐洲定軌中心CODE-GIM、IRI 模型及本文CIT 算法重構TEC 的時空分布特征.

從圖2～5 的對比可以看出：相比于Madrigal 的TEC 產品，CODE 給出的TEC 分布與觀測值基本形態較為符合，但磁赤道駝峰區附近的TEC 要偏低，且TEC 整體分布過于平滑，無法反映局部的TEC 中小尺度的變化特征；IRI 模型計算的TEC 則相比觀測結果存在系統性偏低的情況；而本文CIT 算法重構TEC 則與Madrigal 數據非常一致，且在大部分細節結構上與Madrigal 結果類似.

圖2 2017-09-06 中國及周邊區域Madrigal 給出的TEC 分布Fig.2 Distribution of Madrigal vertical TEC over China and surrounding areas on September 6,2017

圖3 2017-09-06 中國及周邊區域CODE-GIM 數據給出的TEC 分布Fig.3 Ionospheric TEC distribution of CODE-GIM over China and its surrounding regions on September 6,2017

圖4 2017-09-06 中國及周邊區域IRI 模型計算給出的TEC 分布Fig.4 Ionospheric TEC distribution of IRI model over China and its surrounding regions on September 6,2017

圖5 2017-09-06 中國及周邊區域CIT 重構的TEC 分布Fig.5 TEC distribution of CIT algorithm over China and its surrounding regions on September 6,2017

進一步地，對CODE-GIM、IRI 模型及CIT 三種方法的TEC 精度進行統計比較，評估結果如圖6 所示.從分析結果來看：CODE 電離層TEC 數據相對Madrigal 數據的平均誤差為4.0 TECU，標準差為3.2 TECU；IRI 模型的平均誤差為2.9 TECU，標準差為4.4 TECU；CIT 算法的平均誤差和標準差則分別下降為2.5 TECU 和3.0 TECU，精度提升效果非常明顯.

圖6 電離層TEC 地圖精度評估結果Fig.6 Accuracy evaluation results of different ionospheric TEC map

3.2 IED 精度分析

同樣地，根據地基GNSS 和測高儀數據進行聯合CIT，可以獲取區域三維時變的電子密度分布.圖7給出了2017-09-06T00:00UT 時刻中國及周邊區域上空三維電子密度剖面.其中圖7(a)所示為電子密度沿著經度面的切片；圖7(b)所示為峰值電子密度；圖7(c)所示為電子密度沿著高度面的切片；圖7(d)所示為電子密度沿著緯度面的切片.從電子密度空間形態來看，CIT 結果較好地再現了該區域電離層隨緯度、經度和高度方向的基本變化特征.同樣地，圖8 給出了經度100°上IED 隨時間-緯度-高度上的變化.可以看出，該區域電子密度的緯度和高度變化特征基本符合電離層的基本變化規律.

圖7 2017-09-06T0:00UT 中國及周邊區域上空CIT 的三維電子密度Fig.7 3D electron density distribution of CIT over China and its surrounding areas at 0:00UT on September 6,2017

圖8 2017-09-06T0:00UT 經度100°上空IED 隨時間-緯度-高度的變化Fig.8 The variation of ionospheric electron density over 100°E at 0:00UT on September 6,2017

進一步地，利用獨立的國防電波觀測站網測高儀的實測數據對CIT 前后的foF2和hmF2精度進行分析.圖9 給出了CIT 前后電離層foF2精度的對比結果.CIT 前，背景模型IRI 的foF2平均誤差為0.89 MHz，標準差為0.96 MHz；僅利用地基GNSS 進行CIT，其foF2平均誤差和標準差分別下降為0.75 MHz 和0.7 MHz；聯合GNSS 和測高儀進行CIT(GNSS＋IonCIT)，其foF2平均誤差和標準差則進一步下降為0.63 MHz 和0.61 MHz.從對比結果來看，僅利用GNSS 進行層析，也能在一定程度上提高電離層foF2的精度，而聯合測高儀數據后，則精度又有進一步的提升.

圖9 CIT 前后電離層foF2 誤差比較Fig.9 Comparison of ionospheric foF2 deviation before and after CIT

同樣的，圖10 給出了CIT 前后電離層hmF2精度的對比結果.從統計結果來看，僅利用地基GNSS 進行CIT，其hmF2平均誤差和標準差相比背景模型幾乎無明顯變化，平均誤差和標準差均在20.6 km 和16.5 km 上下，由此可以看出，如果僅僅考慮利用地基GNSS 進行CIT，將難以提升電子密度的垂直分辨率；而聯合GNSS 和測高儀進行CIT 后，其hmF2平均誤差和標準差則下降為14.8 km 和11.7 km，驗證了聯合測高儀數據后對CIT 精度的提升效果.

圖10 CIT 前后電離層hmF2 誤差比較Fig.10 Comparison of ionospheric hmF2 deviation before and after CIT

4 小 結

針對單獨使用地基GNSS 進行三維CIT 的不足，提出了一種聯合地基GNSS 和測高儀數據的三維CIT 方法.算法綜合了測高儀探測電離層垂直分辨率較高和地基GNSS CIT 水平分辨率較高的優點，以測高儀輔助驅動更新IRI 模型，將驅動后的IRI 模型作為背景電離層模型，再利用改進的ART 算法進行CIT；基于IGS 與中國陸態網地基GNSS 臺站及GIRO 測高儀數據實現了中國及周邊區域三維CIT.對CIT 結果獲取的TEC 和電子密度的形態學分析驗證了方法的可靠性；與獨立觀測的測高儀數據的對比結果表明：僅依賴地基GNSS 進行CIT 無法有效提升hmF2的重構精度，聯合測高儀數據進行CIT 后，電離層foF2和hmF2的重構精度均有明顯提升.

需要指出的是，本文用于CIT 成像精度評估的臺站主要分布于GNSS 臺站較為密集的區域.對于GNSS 觀測站非常稀少的區域，如海洋或者沙漠上空區域，本文CIT 成像算法的適用性還需要進一步的分析驗證，這也是本文下一步的研究方向.

致謝：本文GNSS 觀測數據從IGS 網站和中國大陸構造環境監測網絡(crustal movement observation network of China,CMNOC)獲取.全球電離層地圖CODE 數據從http://ftp.aiub.unibe.ch/CODE/網站獲取.麻省理工學院高精度Madrigal 數據由http:// cedar.openmadrigal.org/下載，作者在此表示感謝.