基于GNSS的武漢區(qū)域電離層TEC建模

吳偉銓 周晨 禹智斌 左平兵 趙家奇 趙正予

（1. 武漢大學(xué)電子信息學(xué)院空間物理系，武漢 430072；2. 哈爾濱工業(yè)大學(xué)空間科學(xué)與應(yīng)用技術(shù)研究院，深圳 518055）

引 言

電離層天頂方向總電子含量(vertical total electron content, VTEC)是衛(wèi)星導(dǎo)航系統(tǒng)的重要產(chǎn)品. 近年來不少優(yōu)秀的學(xué)者利用經(jīng)典的電離層模型，對VTEC進(jìn)行了建模研究. Lanyi等[1]提出利用2階或3階多項(xiàng)式模型計算單站VTEC；Schaer等[2]利用球諧級數(shù)模型進(jìn)行全球電離層建模；Yuan等[3-4]提出廣義三角級數(shù)函數(shù)模型，該模型在小區(qū)域VTEC 建模中有很好的效果. 在各種電離層模型被提出之后，又有不少專家學(xué)者對這些模型進(jìn)行了比較分析. 文獻(xiàn)[5]對多項(xiàng)式模型、三角級數(shù)模型和低階球諧級數(shù)模型之間的一致性和參數(shù)設(shè)置局限性進(jìn)行了分析，指出了它們的優(yōu)缺點(diǎn). 文獻(xiàn)[6]根據(jù)應(yīng)用環(huán)境對電離層模型進(jìn)行了分類，并提出了模型的一些問題和改進(jìn)措施. 文獻(xiàn)[7]指出了多項(xiàng)式模型和三角級數(shù)模型不適用于大區(qū)域的電離層建模. 文獻(xiàn)[8-9]表明在模型參數(shù)選取適當(dāng)?shù)那闆r下，球諧級數(shù)模型適用于各種大小區(qū)域的電離層建模. 文獻(xiàn)[10]基于球諧級數(shù)模型對山東區(qū)域進(jìn)行電離層建模，并采用單頻精密單點(diǎn)定位驗(yàn)證了該區(qū)域電離層模型的有效性. 電離層在不同高度上的電子密度會有很大差異, 文獻(xiàn)[11]提出了一種三層的電離層模型，可以反映不同高度區(qū)域的VTEC情況. 在電離層建模中，通常要采用一定的估計策略來對模型系數(shù)進(jìn)行估計. 常采用的是分段常數(shù)(piece-wise constant, PWC)估計策略和分段線性(piece-wise linear, PWL)估計策略. 文獻(xiàn)[12]指出，只有在低緯度地區(qū)的電離層模型改正精度才會對估計策略比較敏感，而這種影響是隨緯度升高而減小的. 文獻(xiàn)[13]在總結(jié)全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System, GNSS)電離層延遲修正成果的基礎(chǔ)上，從建模及誤差分析等方面介紹了我國在這一領(lǐng)域的研究進(jìn)展. 但是，上述文獻(xiàn)都沒有對電離層建模方法和估計策略進(jìn)行綜合對比分析，也沒有研究地磁活動對電離層建模的影響. 因此，本文在上述研究的基礎(chǔ)上，對武漢地區(qū)進(jìn)行電離層建模，目的是比較不同建模方法和估計策略在該區(qū)域中的適用性，以及探究地磁活動和電離層建模的相關(guān)性.

國際GNSS服務(wù)(International GNSS Service,IGS)的電離層聯(lián)合分析中心(Ionosphere Associate Analysis Centers, IAACs)可以提供高精度的電離層產(chǎn)品，主要包括歐洲軌道確定中心(Center for Orbit Determination in Europe, CODE)[8]、美國噴氣推進(jìn)實(shí)驗(yàn)室(Jet Propulsion Laboratory, JPL)[14]、歐洲太空局(European Space Agency/European Space Operations Center, ESA/ESOC)[15]、西班牙加泰羅尼亞理工大學(xué)(Universitat Politècnica de Catalunya, UPC)[16]、加拿大自然資源部門(Natural Resources Canada, NRCan)、中國科學(xué)院測量與地球物理研究所和光電研究院(Chinese Academy of Science, CAS)[17]及武漢大學(xué)衛(wèi)星導(dǎo)航與定位技術(shù)研究中心(Wuhan University, WHU)[18].本文選擇CODE提供的電離層產(chǎn)品作為參考源，其采用15階球諧級數(shù)進(jìn)行全球的電離層建模，能提供時間分辨率為1 h、空間分辨率為2.5°×5°的全球電離層格網(wǎng)地圖(global ionospheric map, GIM)[19]. 選取以武漢為中心的十七個連續(xù)運(yùn)行參考站(continuous operation reference system, CORS)在不同地磁活動情況下的數(shù)據(jù)進(jìn)行試驗(yàn)，并與CODE公布的GIM作對比，對于分析建立武漢地區(qū)的電離層模型、削弱電離層延遲誤差具有重要意義[10].

傳統(tǒng)的二維電離層模型是區(qū)域電離層建模的經(jīng)典手段，但是，它們只能根據(jù)經(jīng)度和緯度擬合區(qū)域內(nèi)的VTEC，卻不能反映高度上的電子密度變化情況.而通過三層電離層模型可以從經(jīng)度、緯度、高度三個維度上對電離層進(jìn)行建模，并且可以有效提高電離層地圖的擬合精度[11]. 此外，多層電離層模型對空間目標(biāo)監(jiān)視和低軌衛(wèi)星航天測控具有重要應(yīng)用價值.本文通過研究三種二維電離層模型和估計策略對VTEC計算精度的影響，探索適用于武漢區(qū)域的電離層模型、模型參數(shù)和估計策略，以及分析模型精度和地磁活動的相關(guān)性. 在此研究的基礎(chǔ)上，進(jìn)一步通過三層電離層模型，來分析不同高度區(qū)域的VTEC情況.

1 電離層模型

二維電離層模型包括多項(xiàng)式模型、三角級數(shù)模型、球諧級數(shù)模型、樣條函數(shù)模型、廣義三角級數(shù)模型. 下面主要介紹多項(xiàng)式模型和球諧級數(shù)模型，以及在二維電離層模型基礎(chǔ)上提出的三層電離層模型.

1.1 多項(xiàng)式模型

多項(xiàng)式模型[1]將VTEC看作如下表達(dá)式：

式中：Tv(φ,S)是電離層穿刺點(diǎn)(ionosphere piece point,IPP)處的VTEC；m和n是 模型階數(shù)；Eij是模型系數(shù)；φ是地理緯度； φ0是IPP覆蓋中心的地理緯度；S是太陽時角；S0是IPP覆蓋中心的太陽時角.

太陽時角S具有如下形式：

式中：θ是地理經(jīng)度，單位是(°)；t是UTC時間，單位是h.

通常求解一個多項(xiàng)式模型用2～3 h，所以要每2～3 h計算一組模型系數(shù).

1.2 球諧級數(shù)模型

球諧級數(shù)模型[2]的表達(dá)式如下：

式中：n是模型階數(shù)；(sinφ) 是i度j階的締和勒讓德函數(shù)；φ是IPP的緯度；s是IPP的經(jīng)度；Aij和Bij是模型系數(shù).

求解球諧級數(shù)模型也是使用2～3 h. 區(qū)域電離建模時，可以采用低階的球諧級數(shù)模型.

1.3 三層電離層模型

三層電離層模型[11]把總電子含量(total electron content, TEC)看成與經(jīng)度θ、緯度φ和高度h都相關(guān)的函數(shù)，其表達(dá)式如下：

式中：T是發(fā)送機(jī)到接收機(jī)路徑上的TEC；Rtra和Rrec分別是發(fā)送機(jī)和接收機(jī)的位置；Ne(θ,φ,h)是電子密度函數(shù). 即TEC是電子密度沿發(fā)送機(jī)和接收機(jī)路徑的積分. 電子密度函數(shù)可看成水平基函數(shù)和垂直基函數(shù)的線性組合，表達(dá)式如下：

式中：cij是 參數(shù)；Hi是水平方向的電子密度基函數(shù)，可沿用在二維電離層模型中使用的函數(shù)；Vj是垂直方向的電子密度基函數(shù)，可采用國際參考電離層(international reference ionosphere, IRI)模型[20]輸出的電子密度作為該函數(shù)；NH和NV是兩個基函數(shù)的階數(shù).

三層電離層模型把VTEC看作如下表達(dá)式：

式中：Tv(φ,S,h1,h2)代 表緯度φ、 太陽時角S、高度h1到h2之 間的VTEC；Tv(φ,S,0,∞)是二維電離層模型計算的結(jié)果；cj是參數(shù). 即通過垂直方向電子密度函數(shù)的積分可以得到與高度、緯度、經(jīng)度都相關(guān)的VTEC.

2 武漢地區(qū)電離層建模

本文對武漢周邊17個CORS 2016-01-01—05的數(shù)據(jù)進(jìn)行建模，如圖1所示. 圖2是2015-12-31—2016-01-07的地磁活動指數(shù)(Dst和Kp). 可以看到，01-01處在磁暴恢復(fù)相期，而01-05處在地磁平靜時期，選取這兩天數(shù)據(jù)有利于研究不同地磁活動水平對區(qū)域電離層建模的影響.

圖1 武漢周邊17個CORS分布圖Fig. 1 Distribution map of 17 continuously operating reference stations around Wuhan

圖2 2015-12-31—2016-01-07 Dst指數(shù)和Kp指數(shù)Fig. 2 Dst index and Kp index from Dec. 31st,2015 to Jan. 7th, 2016

2.1 建模方法

本文基于GNSS觀測數(shù)據(jù)，使用載波相位平滑偽距法[14]提高偽距觀測精度，衛(wèi)星截止仰角設(shè)置為15°，投影函數(shù)采用標(biāo)準(zhǔn)投影函數(shù)，時間采樣率為30 s，電離層薄層高度設(shè)置為450 km.

本文利用參考文獻(xiàn)[1]中提出的2階多項(xiàng)式模型、文獻(xiàn)[5]中使用的4階球諧級數(shù)模型，以及文獻(xiàn)[21]中使用的4×3階多項(xiàng)式模型這三種模型進(jìn)行區(qū)域VTEC建模. 估計策略分別采用PWC估計策略和PWL估計策略[12]. PWC假定VTEC模型系數(shù)在適用時間內(nèi)不變；PWL假定模型系數(shù)在這段時間線性變化，通常需要兩組系數(shù)，分別代表該時間段開始和結(jié)束時的模型系數(shù)，則任意時刻的模型系數(shù)都能由這兩組系數(shù)線性表示. 將VTEC模型系數(shù)和硬件延遲偏差(differential code biases, DCB)聯(lián)合建立法方程[22]，采用最小二乘法，每2 h計算一組系數(shù). 得到聯(lián)合DCB后，通過假設(shè)所有衛(wèi)星的DCB之和為零來分離接收機(jī)和衛(wèi)星DCB. 按照0.1°的間隔劃分網(wǎng)格點(diǎn)，求得模型系數(shù)后再反代入網(wǎng)格點(diǎn)，計算出各網(wǎng)格點(diǎn)處VTEC.

2.2 結(jié)果和分析

對2016-01-01的數(shù)據(jù)進(jìn)行建模，采用PWC估計策略時，4×3階多項(xiàng)式模型、2階多項(xiàng)式模型、4階球諧級數(shù)模型的VTEC建模結(jié)果如圖3所示.

圖3 2016-01-01三種模型采用PWC估計策略的VTEC建模結(jié)果Fig. 3 Modeling diagram of three models and piece-wise constant strategy on Jan. 1st, 2016

對2016-01-05的數(shù)據(jù)進(jìn)行建模，采用PWC估計策略時，4×3階多項(xiàng)式模型、2階多項(xiàng)式模型、4階球諧級數(shù)模型的VTEC建模結(jié)果如圖4所示.

圖4 2016-01-05三種模型采用PWC估計策略的VTEC建模結(jié)果Fig. 4 Modeling diagram of three models and piece-wise constant strategy on Jan. 5th, 2016

從圖3～4的二維電離層模型結(jié)果可以看出：在0:00—6:00UT時，電離層電子濃度逐漸增大，并在6:00UT時達(dá)到峰值，而在16:00—24:00UT電子濃度都偏??；低緯度地區(qū)的電子濃度比高緯度地區(qū)高；發(fā)生磁暴時，電離層電子濃度顯著增大，說明地磁活動會影響電離層的電子濃度.

鑒于CODE在硬件延遲及電離層解算方面的優(yōu)越性，采用CODE公布的GIM作為外部對比源. CODE提供的該區(qū)域2016-01-01和01-05的VTEC結(jié)果如圖5所示.

圖5 2016-01-01和01-05 CODE公布的VTEC建模圖Fig. 5 Modeling diagram published by CODE on Jan. 1st and Jan. 5th, 2016

通過三種模型的VTEC建模圖和CODE公布的建模圖進(jìn)行對比可以發(fā)現(xiàn)，三種模型都能達(dá)到較好的擬合效果. 為了進(jìn)一步分析這些模型在地磁平靜時期和地磁擾動時期的建模精度變化，對武漢單站VTEC在PWC估計策略下的計算結(jié)果和CODE公布值進(jìn)行了比較，結(jié)果如圖6所示.

圖6 2016-01-01和01-05武漢單站VTEC建模結(jié)果和CODE公布結(jié)果對比Fig. 6 Comparison of VTEC modeling results and code results of Wuhan single station on Jan. 1st and Jan. 5th, 2016

從圖6可以看出，三種模型的建模結(jié)果幾乎重疊在一起，且和CODE公布的VTEC變化趨勢基本一致，但VTEC大小依然存在一定的誤差. 原因是CODE進(jìn)行的是全球電離層擬合，在中國區(qū)域用的站不到十個，在小區(qū)域范圍內(nèi)擬合精度不夠.

鑒于此原因，本文將VTEC通過電離層投影函數(shù)轉(zhuǎn)化成斜向總電子含量(slant total electron content,STEC)，與實(shí)測值比較計算30天的殘差值均方根(root mean square, RMS)，以此來分析模型和估計策略的精度. 其中殘差值的表達(dá)式如下：

式中：R是殘差值；TS是STEC的觀測值；fm是電離層投影函數(shù)，定義為STEC和VTEC的比；TV是VTEC；D是DCB.

統(tǒng)計2016-01-01—30的STEC殘差值RMS，在PWL和PWC估計策略下的結(jié)果如圖7所示.

圖7 兩種估計策略下的RMS分布Fig. 7 RMS distribution with two estimation strategies

從圖7的殘差RMS結(jié)果可以看出：無論采用哪種估計策略，效果最好的是4階球諧級數(shù)模型，其次是4×3階多項(xiàng)式模型，而2階多項(xiàng)式模型的效果最差，這與其模型參數(shù)較少有關(guān)；三種模型在PWL估計策略下的效果要優(yōu)于PWC估計策略，說明PWL估計策略更適用于區(qū)域的VTEC建模. 當(dāng)使用PWL估計策略和4階球諧級數(shù)模型時，RMS集中在0～1.5 TECU，區(qū)域VTEC建模達(dá)到較高的精度.

對比01-01和01-05的殘差，發(fā)現(xiàn)01-05三種模型的RMS都比01-01對應(yīng)的RMS要小. 另外，從圖6可以看出，同一時刻，模型的計算結(jié)果會比CODE公布值小，除了在01-01T6:00UT，此時VTEC也恰好達(dá)到峰值. 從以上的分析可以推斷出，地磁擾動會對VTEC計算結(jié)果產(chǎn)生干擾. 一方面，地磁活動強(qiáng)的時候GNSS觀測誤差會增大；另一方面，VTEC模型的計算精度也會降低.

谷子生育期為80～140 d，適合小麥、大麥等收后復(fù)種。筆者在內(nèi)蒙古農(nóng)牧業(yè)科學(xué)院試驗(yàn)地曾進(jìn)行飼用谷子不同生長時期全株?duì)I養(yǎng)物質(zhì)測定試驗(yàn)，結(jié)果表明：初花期預(yù)干全株的蛋白質(zhì)含量要明顯高于成熟期的蛋白質(zhì)含量，而初花期預(yù)干全株的纖維含量要明顯低于成熟期纖維含量，初花期預(yù)干全株的干物質(zhì)量與成熟期相當(dāng)，所以飼用谷子在初花期就可以進(jìn)行收獲，大大縮短了生長周期，是短日期作物收獲后復(fù)種的優(yōu)良作物之一。

以上都是基于經(jīng)度和緯度的二維電離層建模.然而，二維電離層模型不能反映VTEC與高度之間的關(guān)系. 為了進(jìn)一步研究底部電離層、峰高電離層和頂部電離層的VTEC情況，本文在上述研究的基礎(chǔ)上利用三層電離層模型對武漢地區(qū)進(jìn)行建模. 參考文獻(xiàn)[11]中對電離層高度的分層方法，將電離層分為三個區(qū)域：90～350 km，350～600 km，600～1 500 km. 利用IRI模型計算出每個IPP三個高度區(qū)域的VTEC比例，再代入4階球諧級數(shù)模型在PWC估計策略下的計算結(jié)果中，構(gòu)建出三層電離層模型.2016-01-01和01-05的建模結(jié)果如圖8所示.

圖8 2016-01-01和01-05三層電離層模型建模Fig. 8 Modeling diagram of three layer ionosphere model on Jan. 5th, 2016

從圖8可以看出：90～350 km高度區(qū)間的電離層TEC最大，其次是350～600 km，而在600 km以上，TEC會變得很稀少. 在得出該結(jié)果后，本文進(jìn)一步研究三層電離層VTEC隨時間的變化規(guī)律. 利用2015-12-31—2016-01-06七天的數(shù)據(jù)，對武漢上空的VTEC進(jìn)行了計算，得出三個高度區(qū)間以及總的VTEC結(jié)果如圖9所示.

從圖9可以看出：在每天VTEC達(dá)到峰值時刻前后，底部電離層的VTEC要明顯大于峰高電離層和頂部電離層的VTEC. 當(dāng)VTEC在谷值時，底部VTEC和峰高VTEC較為接近. 同時，不同高度的VTEC對地磁活動變化的響應(yīng)也有明顯區(qū)別，底部TEC和峰高TEC受地磁活動的影響較大.

圖9 2015-12-31—2016-01-06武漢上空三層電離層VTEC結(jié)果Fig. 9 VTEC results of three layers over Wuhan from Dec.31st, 2015 to Jan. 6th, 2016

3 結(jié) 論

本文針對區(qū)域電離層建模的三種模型和兩種估計策略，分別選取地磁平靜和地磁擾動的兩天數(shù)據(jù)進(jìn)行計算，并以CODE公布的GIM和STEC實(shí)測值作為對比源進(jìn)行分析，得出一些結(jié)論.

1)在區(qū)域VTEC建模中，三種二維電離層模型的建模結(jié)果和CODE公布的GIM基本一致，達(dá)到較好的擬合精度.

2)地磁擾動會對模型的擬合結(jié)果產(chǎn)生干擾，降低模型的擬合精度.

4)當(dāng)使用PWC估計策略時，所有模型的擬合效果都會變差. 因此，在該區(qū)域VTEC建模中建議使用PWL估計策略.

5)從三層電離層模型結(jié)果可以看出，電離層的電子含量集中在90～600 km高度區(qū)間，不同高度的VTEC對地磁活動變化的響應(yīng)有很大區(qū)別.

文中的試驗(yàn)是基于武漢周邊地區(qū)兩天的數(shù)據(jù)來進(jìn)行的，沒有考慮不同的地區(qū)和更長時間的情況. 而且，模型的精度受實(shí)際應(yīng)用場景和觀測數(shù)據(jù)密集程度等多種因素的影響，為了達(dá)到詳細(xì)的分析對比，需要進(jìn)行更加完善的試驗(yàn)設(shè)計.