基于CHAMP、GRACE和COSMIC掩星數(shù)據(jù)的全球電離層hmF2建模研究

劉楨迪， 方涵先*， 翁利斌， 馬強(qiáng)， 張建彬

1 解放軍理工大學(xué)氣象海洋學(xué)院， 南京　211101 2 中國(guó)人民解放軍63655部隊(duì)， 烏魯木齊　841700

?

基于CHAMP、GRACE和COSMIC掩星數(shù)據(jù)的全球電離層hmF2建模研究

劉楨迪1， 方涵先1*， 翁利斌1， 馬強(qiáng)1， 張建彬2

1 解放軍理工大學(xué)氣象海洋學(xué)院， 南京211101 2 中國(guó)人民解放軍63655部隊(duì)， 烏魯木齊841700

基于CHAMP、GRACE和COSMIC的電離層hmF2掩星數(shù)據(jù)，采用最小二乘法建立了一個(gè)包含地磁和太陽(yáng)活動(dòng)、經(jīng)度、地方時(shí)、年積日和緯度信息的非線性多項(xiàng)式hmF2模型(Nonlinear Polynomial Peak Height Model—NPPHM).利用GRACE掩星數(shù)據(jù)對(duì)NPPHM與IRI2012進(jìn)行了獨(dú)立檢驗(yàn)，結(jié)果顯示這兩個(gè)模型在2008年與GRACE數(shù)據(jù)的相關(guān)系數(shù)分別為0.798和0.532，均方根誤差分別為25.97 km和44.56 km；在2012年，相關(guān)系數(shù)分別為0.732和0.488，均方根誤差分別為31.39 km和42.83 km.選取全球不同地區(qū)14個(gè)測(cè)高儀站點(diǎn)數(shù)據(jù)，并引入相似離度對(duì)這兩個(gè)模型的精度進(jìn)行了評(píng)估，結(jié)果表明，NPPHM的相似離度遠(yuǎn)小于IRI2012，更加接近測(cè)高儀觀測(cè)值.使用Athens站2003—2013年數(shù)據(jù)對(duì)模型進(jìn)行了檢驗(yàn)，結(jié)果顯示IRI2012與Athens站測(cè)高儀數(shù)據(jù)的平均偏差達(dá)到8.11%，NPPHM則只有3.53%，并且在太陽(yáng)活動(dòng)低年及每年10月NPPHM的精度要明顯高于IRI2012.此外，NPPHM也能夠較好模擬出hmF2的日變化、季節(jié)變化及赤道異常等特性.

電離層hmF2建模；非線性多項(xiàng)式；掩星數(shù)據(jù)

1　引言

電離層F2層的電子密度峰值高度(hmF2)是電離層物理中一個(gè)十分重要的參量.它對(duì)高頻電波在電離層內(nèi)的傳播和反射具有重要影響.hmF2的抬升或下降將直接影響到信號(hào)在電離層內(nèi)的傳播，可能改變其傳輸路徑，造成信號(hào)中斷.因此，對(duì)hmF2的研究將有助于無(wú)線電波頻譜管理.在構(gòu)建電離層電子密度剖面時(shí)，hmF2是其中一個(gè)至關(guān)重要的因素.所以，對(duì)峰值高度hmF2的建模具有重要的科學(xué)意義.

電離層hmF2存在較為復(fù)雜的變化，它在赤道地區(qū)的分布范圍一般是350～500 km，在中緯地區(qū)一般是250～350 km.hmF2也存在依賴于季節(jié)、地方時(shí)、太陽(yáng)活動(dòng)和地磁活動(dòng)的變化特性.此外，電場(chǎng)漂移及中性風(fēng)等對(duì)hmF2也具有一定的調(diào)制作用(Hoque and Jakowski，2012).

目前，國(guó)際參考電離層模型IRI是應(yīng)用廣泛的hmF2經(jīng)驗(yàn)?zāi)Ｐ?IRI模型是根據(jù)大量電離層探測(cè)數(shù)據(jù)得到的經(jīng)驗(yàn)預(yù)報(bào)模型，它綜合了全球180多個(gè)數(shù)字測(cè)高儀站點(diǎn)的資料，主要描述了海拔高度50～1500 km范圍內(nèi)平靜地磁場(chǎng)條件下非極區(qū)電離層的峰值高度、電子密度、電子溫度、離子成分、離子溫度和離子漂移的月平均值等重要參數(shù)(Bilitza et al., 1979).繼Rawer(1978)發(fā)布了第一個(gè)版本的IRI之后，學(xué)者們又提出了多個(gè)IRI的更新版本(Bilitza et al., 1979; Bilitza, 1990, 2015; Bilitza and Reinisch, 2008)，目前國(guó)際參考電離層的最新版本為IRI2012.IRI利用F2層傳輸因子M(3000)F2參量來(lái)得到hmF2，二者的轉(zhuǎn)換方程只有在近似假設(shè)電子密度與高度呈拋物型相關(guān)的條件下才能成立(Shubin et al., 2013)，但這一假設(shè)與實(shí)際情況并不完全相符.Magdaleno 等(2011)的研究表明，在高緯地區(qū)IRI模型值超出測(cè)高儀實(shí)測(cè)值的30%左右，在赤道附近則低估達(dá)40%左右.

近來(lái)，很多學(xué)者利用不同的觀測(cè)數(shù)據(jù)和方法對(duì)hmF2的建模進(jìn)行了研究.例如，Zhang等(2009，2010)使用經(jīng)驗(yàn)正交函數(shù)(EOF)、Altadill等(2013)使用球諧函數(shù)構(gòu)造了全球hmF2經(jīng)驗(yàn)?zāi)Ｐ停珺runini等(2013)則分別使用不同的數(shù)學(xué)表達(dá)來(lái)描述hmF2的全球變化和日變化.但是，這些模型都是基于地基測(cè)高儀數(shù)據(jù)建立的.雖然地面測(cè)站已經(jīng)積累了較長(zhǎng)時(shí)間的電離層觀測(cè)數(shù)據(jù)，但全球站點(diǎn)的數(shù)量極為有限，且全球分布不均(在海洋上幾乎沒(méi)有分布)，這是以上建模的不足之處.隨著掩星技術(shù)的發(fā)展，研究人員基于掩星數(shù)據(jù)也構(gòu)建了一些hmF2模型，例如Shubin等(2013)使用勒讓德函數(shù)和傅里葉展開(kāi)分別描述hmF2的空間變化和時(shí)間變化，得到了SMF2(the Satellite Model of the F2layer)模型，但是，SMF2模型缺少太陽(yáng)活動(dòng)變化項(xiàng)，只能局限于太陽(yáng)活動(dòng)低年使用.Hoque和Jakowski(2012)利用非線性多項(xiàng)式方程，建立了hmF2的NPHM(Neustrelitz Peak Height Model)模型.然而，NPHM模型只包含13個(gè)系數(shù)，并且模型缺少經(jīng)度變化項(xiàng)和地磁變化項(xiàng)，它與NeQuick模型(Radicella and Leitinger，2001)的比較表明兩者的相對(duì)誤差是接近的，可見(jiàn)NPHM在精度上提高不大.

本文利用CHAMP(Challenging Minisatellite Payload), GRACE(Gravity Recovery And Climate Experiment)和COSMIC(Constellation Observing System for Meteorology, Ionosphere and Climate)的掩星探測(cè)數(shù)據(jù)，基于目前對(duì)電離層hmF2的研究基礎(chǔ)，使用非線性多項(xiàng)式方程對(duì)全球hmF2建模.利用非線性多項(xiàng)式方程能夠構(gòu)建出形式相對(duì)簡(jiǎn)單、計(jì)算效率相對(duì)較快，同時(shí)精度較高的模型(Hoque and Jakowski, 2012).同NPHM模型相比，本文采用的非線性多項(xiàng)式包含315個(gè)系數(shù)，且除緯度、太陽(yáng)活動(dòng)、年積日和地方時(shí)外，更為全面的考慮了經(jīng)度和地磁活動(dòng)兩項(xiàng)影響因素.此外，本文將全球按照緯度劃分為71個(gè)緯度帶，在各緯度帶內(nèi)分別通過(guò)最小二乘擬合確定系數(shù)，從而回避了在模型中加入先驗(yàn)緯度變化信息，進(jìn)一步提高模型的計(jì)算效率和準(zhǔn)確性.使用GRACE數(shù)據(jù)、地面測(cè)高儀數(shù)據(jù)對(duì)模型進(jìn)行了獨(dú)立檢驗(yàn)，并利用NPPHM模擬了hmF2變化特性，評(píng)估NPPHM的精度和模擬能力.

2　數(shù)據(jù)

本文使用了掩星數(shù)據(jù)進(jìn)行建模和檢驗(yàn)，這部分?jǐn)?shù)據(jù)源由CHAMP衛(wèi)星2001年6月到2008年9月、GRACE衛(wèi)星2007年2月到2013年9月和COSMIC衛(wèi)星2006年4月到2014年12月的掩星剖面(EDP)組成.Hajj和Romans(1998)、Garcia-Fernandez 等(2003)和Lei等(2007)對(duì)掩星反演技術(shù)做出了詳細(xì)說(shuō)明.通過(guò)Abel變換可以從掩星觀測(cè)中獲得十分準(zhǔn)確的F2層峰值高度(hmF2)和峰值密度(NmF2).Yue等(2010)指出反演得到的hmF2與真實(shí)值之間的絕對(duì)標(biāo)準(zhǔn)差為8.9 km(夜間)和7.4 km(白天)，相對(duì)標(biāo)準(zhǔn)差為2%(夜間)和2%(白天).Lei 等 (2007)將掩星反演得到的hmF2與Millstone Hill站的非相干散射雷達(dá)觀測(cè)數(shù)據(jù)及Jicamarca站的垂測(cè)儀觀測(cè)數(shù)據(jù)進(jìn)行比較后發(fā)現(xiàn)，掩星反演的偏差范圍為10～30 km.Krankowski等(2011)利用部分歐洲垂測(cè)儀數(shù)據(jù)對(duì)掩星反演的hmF2做了進(jìn)一步研究發(fā)現(xiàn)，兩種數(shù)據(jù)的平均偏差為2.8 km.CHAMP和GRACE衛(wèi)星的掩星剖面由German Aerospace Center(DLR)Neustrelitz(http:∥swaciweb.dlr.de/)提供.COSMIC掩星剖面由COSMIC Data Analysis and Archive Center(CDAAC,http:∥cosmic-io.cosmic. ucar.edu/cdaac/index. html)提供.從掩星剖面中，可以獲取各個(gè)高度處的電子密度及精度較高的地理坐標(biāo).對(duì)于得到的剖面，根據(jù)Lei等(2007)的研究，本文篩選出，150 km≤hmF2≤600 km，103·cm-3≤NmF2≤8×106·cm-3，并且衛(wèi)星軌道高度(orbalt)滿足orbalt≥hmF2+50的剖面，從中提取出對(duì)應(yīng)的hmF2.最終由COSMIC得到有效數(shù)據(jù)約3.8×106個(gè)(約占總數(shù)據(jù)的92%)，由CHAMP得到有效數(shù)據(jù)約2.2×105個(gè)(約占總數(shù)據(jù)的5%)，由GRACE得到有效數(shù)據(jù)約1.3×105個(gè)(約占總數(shù)據(jù)的3%).圖1表示選用數(shù)據(jù)隨太陽(yáng)活動(dòng)水平、地理緯度、地理經(jīng)度、年積日、地方時(shí)和地磁指數(shù)的分布情況.從圖1可知，建模所用數(shù)據(jù)對(duì)經(jīng)度、年積日和地方時(shí)的覆蓋較為全面和均勻；中低緯數(shù)據(jù)較多，高緯數(shù)據(jù)較少；低太陽(yáng)活動(dòng)水平分布較多，高太陽(yáng)活動(dòng)水平分布較少，F(xiàn)10.7>150時(shí)分布極少；地磁活動(dòng)弱時(shí)的數(shù)據(jù)較多，地磁活動(dòng)強(qiáng)時(shí)的數(shù)據(jù)較少，Ap>20的數(shù)據(jù)極少.

圖1　數(shù)據(jù)隨太陽(yáng)活動(dòng)、年積日、地理緯度、地方時(shí)、地理經(jīng)度、地磁指數(shù)的分布情況Fig.1　Data distributions as F10.7, day of year, Geographical Latitude, Local Time，Geographical Longitude and geomagnetic index Ap

另外，本文也使用了地基測(cè)高儀hmF2數(shù)據(jù)進(jìn)行檢驗(yàn).這部分?jǐn)?shù)據(jù)由Space Physics Interactive Data Resource (SPIDR,http:∥spidr.ngdc.noaa.gov/spidr/)提供.我們選取了14個(gè)能夠代表全球不同區(qū)域的測(cè)站.由于測(cè)高儀數(shù)據(jù)只用于對(duì)模型的檢驗(yàn)，所以具體的測(cè)站坐標(biāo)及數(shù)據(jù)內(nèi)容將在第4節(jié)予以說(shuō)明.

3　方法和結(jié)果

3.1建模方法

采用非線性多項(xiàng)式方程對(duì)全球hmF2建模.非線性多項(xiàng)式方程已經(jīng)被成功應(yīng)用到頂部電離層標(biāo)高模型、電離層電子總含量模型及熱層大氣密度模型等多個(gè)經(jīng)驗(yàn)?zāi)Ｐ椭?Kutiev et al., 2006;Kakinami et al., 2008,2009; Liu et al., 2013).

為提高計(jì)算效率，同時(shí)避免在模型中加入先驗(yàn)緯度變化信息，本文采用劃分緯度帶分別建模的方法.將全球hmF2數(shù)據(jù)從90°S開(kāi)始向北按照地理緯度間隔5°進(jìn)行劃分.同時(shí)，為避免模型出現(xiàn)緯度變化不光滑的情況，使相鄰緯度帶之間重合2.5°的區(qū)域，最終將全球劃分為71個(gè)緯度帶.在各緯度帶內(nèi)分別擬合求相應(yīng)的系數(shù).每個(gè)緯度帶的中間緯度作為該緯度帶擬合所得系數(shù)對(duì)應(yīng)的緯度，其他緯度處的系數(shù)可以通過(guò)插值得到.表1給出了前5個(gè)緯度帶的具體情況.

表1　前5個(gè)緯度帶范圍及所得系數(shù)對(duì)應(yīng)的緯度Table 1　Ranges of the first 5 latitude zones and the corresponding latitude for each group of coefficients

除緯度外，影響hmF2的因素主要包括地理經(jīng)度(Lon)、太陽(yáng)活動(dòng)(F10.7)、地磁活動(dòng)(Ap)、地方時(shí)(LT)和年積日(Doy).非線性多項(xiàng)式模型是基于這些影響因素都是相互獨(dú)立的假定建立的.但是在研究中我們發(fā)現(xiàn)，將太陽(yáng)活動(dòng)、地磁活動(dòng)和地理經(jīng)度作為獨(dú)立項(xiàng)時(shí)，模型精度并沒(méi)有得到明顯提高，然而計(jì)算效率卻大大下降，所以本文將這三項(xiàng)因素作為一個(gè)獨(dú)立多項(xiàng)式加以考慮，這與Liu等(2013)的處理方法相似.具體模型形式如下：

hmF2=F1(F10.7,Ap,Lon)F2(LT)F3(Doy),

(1)

方程(1)中F1-3是分別表示各項(xiàng)變化且含有多個(gè)系數(shù)的多項(xiàng)式.F1表征hmF2隨太陽(yáng)活動(dòng)(F10.7)、地磁指數(shù)(Ap)和地理經(jīng)度(Lon)的變化.根據(jù)Liu 等 (2005) 和Müller 等 (2009) 的研究，地磁活動(dòng)項(xiàng)取為線性關(guān)系.Shubin (2015)指出太陽(yáng)活動(dòng)與hmF2之間近似為對(duì)數(shù)關(guān)系.所以F1具體形式如下：

(2)

F1含有(2M+3)個(gè)系數(shù).

F2表征hmF2隨地方時(shí)(LT)的日變化：

(3)

F2含有(2N+1)個(gè)系數(shù).

F3表征hmF2隨年積日(Doy)的季節(jié)變化：

(4)

F3含有(2P+1)個(gè)系數(shù).

需要注意的是，利用最小二乘原理進(jìn)行擬合確定系數(shù)時(shí)，先根據(jù)方程(2)—(4)將(1)式展開(kāi).方程(2)—(4)相乘后總共得到[(2M+3)×(2N+1)×(2P+1)]個(gè)獨(dú)立系數(shù).與Hoque和Jakowski(2012)的方法相似，方程(3)中，取N=3，分別表示hmF2的日變化、半日變化和1/3日變化；根據(jù)Lin等(2014)的研究，方程(4)中，取P=2，分別表示hmF2的年變化和半年變化.為獲取較高的精度，同時(shí)提高模型計(jì)算效率，本文對(duì)方程(1)中M的不同取值進(jìn)行了研究對(duì)比，最終選取M=3.模型在全球的71個(gè)緯度帶內(nèi)共有9×7×5×71=22365個(gè)系數(shù)，它們由掩星數(shù)據(jù)通過(guò)最小二乘法擬合得到.

3.2建模結(jié)果

為了檢驗(yàn)擬合效果，需要對(duì)NPPHM的殘差進(jìn)行分析.根據(jù)建模使用的掩星數(shù)據(jù)源對(duì)應(yīng)的緯度、經(jīng)度、年積日、地方時(shí)、太陽(yáng)活動(dòng)和地磁活動(dòng)，得到NPPHM及IRI2012對(duì)應(yīng)的模型值(mod)，計(jì)算模型值與原始數(shù)據(jù)(obs)間的相對(duì)偏差(RD)、平均相對(duì)偏差(Mean)和均方差(STD)：

(5)

(6)

(7)

(8)

圖2給出了NPPHM和IRI2012相對(duì)偏差的統(tǒng)計(jì)結(jié)果，實(shí)線是NPPHM的統(tǒng)計(jì)結(jié)果，虛線是IRI2012的統(tǒng)計(jì)結(jié)果.IRI2012相對(duì)偏差偏向正值區(qū)分布且在偏差較大區(qū)域分布較多，平均相對(duì)偏差達(dá)7.3%，均方差為14.99%.而NPPHM偏差呈現(xiàn)正態(tài)分布特征，大部分集中在0附近，向正負(fù)兩側(cè)迅速減少，平均相對(duì)偏差為0.81%，均方差為9.57%.由此可見(jiàn)，建模的結(jié)果是比較理想的.

4　模型檢驗(yàn)

利用建模數(shù)據(jù)源以外的掩星數(shù)據(jù)和電離層測(cè)高儀數(shù)據(jù)，對(duì)模型進(jìn)行獨(dú)立檢驗(yàn)，并同時(shí)與IRI2012的精度進(jìn)行對(duì)比.此外，利用NPPHM模擬hmF2的形態(tài)特征以檢驗(yàn)其表征hmF2變化特性的性能.

圖2　相對(duì)偏差分布直方圖Fig.2　Histogram of relative deviation

4.1GRACE掩星數(shù)據(jù)獨(dú)立檢驗(yàn)

分別選取GRACE 2008年和2012年的掩星數(shù)據(jù)，與NPPHM和IRI2012進(jìn)行比對(duì).需要指出的是，這一步進(jìn)行檢驗(yàn)的模型，是將2008年和2012年GRACE掩星數(shù)據(jù)分別從建模原始數(shù)據(jù)源中剔除后重新建模所得到的結(jié)果，所以驗(yàn)證數(shù)據(jù)是獨(dú)立的.2008年和2012年分別是太陽(yáng)活動(dòng)低年和較高年.此外，這些數(shù)據(jù)涵蓋了大量不同的經(jīng)緯度、地方時(shí)和年積日.圖3和圖4分別是2008年和2012年的比對(duì)結(jié)果，同時(shí)計(jì)算給出了相關(guān)系數(shù)(R)、平均相對(duì)偏差(Mean)和均方根誤差(RMS).均方根誤差(RMS)計(jì)算方法如下：

(9)

由圖3可以看出，2008年太陽(yáng)活動(dòng)低年期間，NPPHM與GRACE掩星觀測(cè)值的相關(guān)系數(shù)為0.798，高于IRI2012的0.532；兩個(gè)模型的平均相對(duì)偏差分別為1.93%和11.48%，表明NPPHM的偏差遠(yuǎn)小于IRI2012；均方根誤差分別為25.97 km和44.56 km，NPPHM比IRI2012降低了約42%.

圖3　2008年NPPHM與IRI2012 hmF2同GRACE hmF2的散點(diǎn)分布圖Fig.3　Scatter plots of NPPHM results and IRI2012 results versus GRACE hmF2 in 2008

圖4所示的2012年是太陽(yáng)活動(dòng)水平較高年份，NPPHM和IRI2012與觀測(cè)值的相關(guān)系數(shù)分別為0.732和0.488，前者較后者提高了約50%；平均相對(duì)偏差分別為0.19%和5.97%；均方根誤差分別為31.39 km和42.83 km，前者較后者降低了約27%.

從GRACE掩星數(shù)據(jù)的檢驗(yàn)結(jié)果來(lái)看，IRI2012經(jīng)驗(yàn)?zāi)Ｊ脚c掩星觀測(cè)數(shù)據(jù)相比存在較大的偏差，而NPPHM模型的偏差較小，尤其在太陽(yáng)活動(dòng)低年，NPPHM的偏差進(jìn)一步減小，這可能是因?yàn)榻Ｊ褂玫脑紨?shù)據(jù)大部分集中在太陽(yáng)活動(dòng)低年，高年數(shù)據(jù)較少.需要指出的是，當(dāng)電離層受到磁暴等異?；顒?dòng)影響時(shí)，峰值高度hmF2可能會(huì)出現(xiàn)極高或極低值，然而這部分的數(shù)據(jù)量很少且電離層暴變化非常復(fù)雜，所以NPPHM與IRI2012模型均在這些情況下出現(xiàn)了較大的偏差，這是經(jīng)驗(yàn)?zāi)Ｐ偷慕Ｔ硭鶝Q定的.

4.2電離層測(cè)高儀數(shù)據(jù)檢驗(yàn)

4.2.1全球hmF2分布檢驗(yàn)

選取代表全球不同區(qū)域的電離層測(cè)高儀數(shù)據(jù)，與NPPHM和IRI2012進(jìn)行比對(duì).分別選取了南、北半球各4個(gè)測(cè)站的3、6、9和12月的數(shù)據(jù)，分別代表了春、夏、秋、冬四個(gè)季節(jié).這些測(cè)站代表著南北半球的低、中、高緯地區(qū)及東西半球的美洲、歐洲和大洋洲地區(qū)，表2列出了選取測(cè)站的地理坐標(biāo).為了保證測(cè)站當(dāng)月的hmF2數(shù)據(jù)量足夠多，選取數(shù)據(jù)的時(shí)間段范圍較寬，包括2007年、2008年、2013年和2014年.其中2007年和2008年是太陽(yáng)活動(dòng)低年，2013年和2014年是太陽(yáng)活動(dòng)高年.圖5表示上述4年F10.7指數(shù)隨年積日的變化.

將測(cè)高儀hmF2的月均值(Ionosonde)與NPPHM及IRI2012的模型月均值進(jìn)行比較.為了全面評(píng)估模型值與觀測(cè)值在變化趨勢(shì)及數(shù)值大小上的一致程度，本文引入相似離度進(jìn)行分析.相似離度不僅能夠體現(xiàn)樣本間數(shù)值差異的大小，而且能夠體現(xiàn)樣本間形態(tài)上的相似，是一種比較全面的相似衡量標(biāo)準(zhǔn)(李開(kāi)樂(lè)，1986；翁利斌等，2011)，它在氣象相關(guān)領(lǐng)域得到了廣泛應(yīng)用.計(jì)算公式為

圖4　2012年NPPHM與IRI2012hmF2同GRACE hmF2的散點(diǎn)分布圖Fig.4　Scatter plots of NPPHM results and IRI2012 results versus GRACE hmF2 in 2012

圖5　2007年、2008年、2013年和2014年太陽(yáng)活動(dòng)變化Fig.5　The variations of F10.7 in 2007,2008,2013 and 2014

站點(diǎn)緯度/°N經(jīng)度/°E站點(diǎn)緯度/°N經(jīng)度/°EKwajalein9167.2Jicamarca-12.1-77PuertoRico18.5-67.2Madimbo-22.430.9Athens3823.5Louisvale-28.521.2Boulder40-105.3Hermanus-34.419.2MillstoneHill42.6-71.5Bundoora-37.72145.05Moscow55.537.3PortStanley-51.7-57.8College64.9-147.8Thule77.5-69.2

(10)

其中，

(11)

(12)

(13)

(14)

圖6給出了太陽(yáng)活動(dòng)低年的分析結(jié)果，圖7給出了太陽(yáng)活動(dòng)高年的分析結(jié)果.從圖6可以看出，在太陽(yáng)活動(dòng)低年，IRI2012普遍偏高，而NPPHM在南北半球的春、夏、秋、冬四個(gè)季節(jié)都與測(cè)高儀觀測(cè)值更加接近.NPPHM在反映hmF2的變化趨勢(shì)上也優(yōu)于IRI2012.例如Boulder站LT=23時(shí)左右測(cè)高儀觀測(cè)數(shù)據(jù)顯示出上升趨勢(shì)，NPPHM呈現(xiàn)出相同的趨勢(shì)，而IRI2012則出現(xiàn)了下降趨勢(shì)；Madimbo站LT=12時(shí)左右觀測(cè)數(shù)據(jù)出現(xiàn)下降，但這一變化十分微小且持續(xù)時(shí)間很短，可以忽略，NPPHM保持水平，而IRI2012則出現(xiàn)較為明顯的上升趨勢(shì).從相似離度的分析可以看出，NPPHM與觀測(cè)值的相似離度均小于IRI2012，且有較大提升.例如在Kwajalein站，IRI2012的相似離度是0.285，NPPHM的相似離度是0.156，提升45%以上.

圖6　太陽(yáng)活動(dòng)低年測(cè)高儀(點(diǎn))、NPPHM(實(shí)線)和IRI2012(虛線)hmF2月均值隨地方時(shí)的變化情況上四幅表示北半球，下四幅表示南半球；從左至右分別表示3月、6月、9月和12月；c1、c2分別表示IRI2012和NPPHM同測(cè)高儀的相似離度.Fig.6　Variations of monthly mean of NPPHM (solid lines) and IRI2012 (dotted lines) as a function of LT in comparisons with ionosonde observation (dots) in low solar activity yearsThe upper panel is the Northern Hemisphere, the lower panel is the Southern Hemisphere; fromleft to right is March, June, September and December, respectively; c1 is the analog deviation between IRI2012 and ionosonde observation, c2 is the analog deviation between NPPHM and ionosonde observation.

圖7　同圖6，太陽(yáng)活動(dòng)高年情況Fig.7　Same as Fig.6,but for high solar activity years

由圖7可知，在太陽(yáng)活動(dòng)高年，NPPHM相比于IRI2012更接近實(shí)測(cè)值.從相似離度的分析結(jié)果可以看出，在其中7個(gè)測(cè)站NPPHM的相似離度小于IRI2012.例如，在Puerto Rico站，IRI2012相似離度為0.182，而NPPHM的相似離度為0.097，提升47%左右.但是，在Athens站，NPPHM的相似離度略高于IRI2012，分別為0.069和0.053，這可能是太陽(yáng)活動(dòng)高年的建模數(shù)據(jù)較少造成的.

總體來(lái)看，不論在太陽(yáng)活動(dòng)低年還是太陽(yáng)活動(dòng)高年，NPPHM在絕大部分測(cè)站處的精度優(yōu)于IRI2012，在數(shù)值和形態(tài)上都與測(cè)高儀觀測(cè)值更加相似.同時(shí)我們也注意到，從全球來(lái)看，NPPHM在低緯赤道地區(qū)的相似離度較其他地區(qū)偏高，這可能是由于低緯赤道地區(qū)不滿足掩星數(shù)據(jù)反演過(guò)程中的球?qū)ΨQ假設(shè)，因此產(chǎn)生了一定的反演誤差.下一步我們將引入包括低緯赤道地區(qū)在內(nèi)的垂測(cè)儀hmF2數(shù)據(jù)，以進(jìn)一步改進(jìn)NPPHM的區(qū)域適用性.

4.2.2長(zhǎng)期hmF2變化檢驗(yàn)

選取Athens站(38°N，23.5°E)一個(gè)太陽(yáng)活動(dòng)周(2003—2013年)的觀測(cè)數(shù)據(jù)，對(duì)11年數(shù)據(jù)逐月(Month)逐世界時(shí)(UT)求出hmF2的月均值，以檢驗(yàn)NPPHM在一個(gè)完整的太陽(yáng)活動(dòng)周內(nèi)的表現(xiàn).圖8給出了11年間月均值的變化情況.圖8(中)出現(xiàn)的空缺部分是因?yàn)锳thens站缺少2008年8月、9月、10月和2012年6月所有時(shí)刻的數(shù)據(jù)以及2012年7月部分時(shí)刻的數(shù)據(jù).

由圖8可知，NPPHM能夠較為準(zhǔn)確地反映出11年間Athens站的月均值變化，高低值分布及變化趨勢(shì)都與測(cè)高儀數(shù)據(jù)比較一致.經(jīng)計(jì)算，2003—2013年IRI2012與Athens站測(cè)高儀數(shù)據(jù)的平均偏差(Mean)達(dá)到8.11%，NPPHM則只有3.53%，后者比前者降低了約56.5%.注意到，IRI2012在每年10月的UT=9—12時(shí)之間出現(xiàn)明顯的高值區(qū)，而測(cè)高儀數(shù)據(jù)則相應(yīng)的表現(xiàn)為低值區(qū)，NPPHM與測(cè)高儀數(shù)據(jù)表現(xiàn)一致.表3分別列出了各年10月UT=9、10、11、12時(shí)的月均值偏差.由表3可知，NPPHM的偏差在0附近波動(dòng)，有正有負(fù)，其絕對(duì)值絕大部分在20 km以下，最高為28.79 km；而IRI2012的偏差全部為正值，大部分在50 km以上，最高達(dá)83.37 km.2006—2010年間，測(cè)高儀數(shù)據(jù)在UT=6—18時(shí)較低，大部分在200 km以下，最低達(dá)180 km；在UT=0—5時(shí)及UT=19—23時(shí)之間較高，但大部分也在280 km以下，最高達(dá)300 km.NPPHM與測(cè)高儀觀測(cè)值相似，而IRI2012則明顯偏高，在UT=6—18時(shí)之間普遍偏高20 km以上，且出現(xiàn)了300 km左右的高值，在UT=0—5時(shí)及UT=19—23時(shí)之間，hmF2大多數(shù)分布在300～330 km之間.表4列出了2006—2010年間NPPHM和IRI2012在UT=19—23時(shí)之間的平均偏差.從表4可以看出，NPPHM的平均偏差絕對(duì)值均在2 km以下，而IRI2012均在25 km以上.

表3　UT=9、10、11、12時(shí)的月均值偏差(單位：km)Table 3　The deviation of monthly mean at UT=9,10,11,12 (unit: km)

注：每個(gè)月的上一行表示NPPHM偏差，下一行表示IRI2012偏差.

表4　UT=19—23之間的平均偏差Table 4　Mean deviation from UT=19 to UT=23

4.3hmF2變化特性檢驗(yàn)

為了檢驗(yàn)NPPHM表征hmF2隨緯度(Lat)、地方時(shí)(LT)和年積日(Doy)變化的性能，本文研究了NPPHM在確定經(jīng)度(Lon)、太陽(yáng)活動(dòng)(F10.7)和地磁指數(shù)(Ap)條件下的模型結(jié)果.當(dāng)Lon=70°W，Ap=5時(shí)，圖9(a1—a4)分別給出了F10.7=80條件下春分日(Doy=80)、夏至日(Doy=173)、秋分日(Doy=266)和冬至日(Doy=356)的hmF2分布情況;圖9(b1—b4)給出了F10.7=120條件下相同年積日的分布情況.從圖9可知，白天hmF2在磁赤道附近顯著增強(qiáng)，這與我們熟知的電離層赤道異?，F(xiàn)象相符，并且由圖9可知這種異常在LT=12時(shí)和LT=20時(shí)左右出現(xiàn)兩個(gè)峰值，這與Yue 等(2015)的研究結(jié)果一致.LT=12時(shí)的增強(qiáng)是由于此時(shí)太陽(yáng)天頂角最小，太陽(yáng)輻射最強(qiáng)，熱層大氣電離產(chǎn)生的等離子體在電場(chǎng)作用下向上漂移；LT=20時(shí)的增強(qiáng)是日落期間向上垂直漂移的反轉(zhuǎn)增強(qiáng)造成的.由圖9分析hmF2的日變化可知，磁赤道附近白天hmF2高于夜間，中緯及高緯地區(qū)夜間hmF2高于白天，這與我們已知的結(jié)論是相符的(熊年祿等，1999).造成這種現(xiàn)象的原因是在中高緯地區(qū)，子午風(fēng)在夜間產(chǎn)生很強(qiáng)的向上漂移，在白天產(chǎn)生向下漂移，造成了hmF2的晝夜差，而在低緯地區(qū)這種晝夜差逐漸減弱至反號(hào).由圖9a2、a4、b2、b4可知，hmF2在夏季半球高于冬季半球，這與已知的hmF2季節(jié)變化特征相符(熊年祿等，1999)，這可能是因?yàn)榘滋斓氖⑿酗L(fēng)場(chǎng)是夏季半球吹向冬季半球.

5　結(jié)論

本文利用非線性多項(xiàng)式構(gòu)建了全球hmF2經(jīng)驗(yàn)?zāi)Ｐ蚇PPHM，模型中包含315個(gè)系數(shù)，自變量包括太陽(yáng)活動(dòng)指數(shù)、地磁指數(shù)、地理經(jīng)度、地方時(shí)和年積日.將來(lái)自CHAMP、GRACE和COSMIC的掩星數(shù)據(jù)按照緯度劃分為71個(gè)不同緯度帶，分別在這些緯度帶內(nèi)通過(guò)最小二乘擬合最終得到71組系數(shù).經(jīng)過(guò)對(duì)NPPHM和IRI2012的檢驗(yàn)，本文得到以下結(jié)論：

(1) 分別利用NPPHM和IRI2012對(duì)建模使用的掩星數(shù)據(jù)源進(jìn)行模擬， NPPHM平均偏差遠(yuǎn)小于IRI2012，說(shuō)明建模的擬合效果較好.

(2) 利用建模數(shù)據(jù)源以外的GRACE掩星數(shù)據(jù)獨(dú)立檢驗(yàn).NPPHM和IRI2012在太陽(yáng)活動(dòng)低年(2008年)與GRACE數(shù)據(jù)的相關(guān)系數(shù)分別為0.798和0.532，均方根誤差分別為25.97 km和44.56 km；在太陽(yáng)活動(dòng)高年(2012年)的相關(guān)系數(shù)分別為0.732和0.488，均方根誤差分別為31.39 km和42.83 km.

(3) 選取代表全球不同區(qū)域的地面測(cè)高儀數(shù)據(jù)并引入相似離度對(duì)NPPHM和IRI2012進(jìn)行對(duì)比分析表明，NPPHM在全球絕大部分測(cè)站優(yōu)于IRI2012，它在數(shù)值和變化形態(tài)上都與測(cè)高儀觀測(cè)更加相似.

(4) 使用測(cè)高儀Athens站2003—2013年hmF2的月均值與NPPHM及IRI2012比較.長(zhǎng)期來(lái)看，IRI2012與測(cè)高儀數(shù)據(jù)偏差較大，達(dá)8.11%，而NPPHM高低值分布及變化趨勢(shì)都與測(cè)高儀數(shù)據(jù)比較一致，且偏差僅為3.53%.

圖8　2003—2013年Athens站IRI2012(上)、測(cè)高儀(中)和NPPHM(下)月均值隨世界時(shí)的變化Fig.8　Diurnal variations of monthly mean of IRI2012 (top)，ionosonde observation (middle) and NPPHM (bottom) as a function of UT from 2003 to 2013 in Athens

圖9　NPPHM隨緯度和地方時(shí)的變化(Lon=70°W，Ap=5)上一行表示太陽(yáng)活動(dòng)低年(F10.7=80)，下一行表示太陽(yáng)活動(dòng)高年(F10.7=120)；從左至右分別表示春分日(Doy=80)、夏至日(Doy=173)、秋分日(Doy=266)和冬至日(Doy=356).Fig.9　Variations of NPPHM as functions of latitude and local time(Lon=70°W，Ap=5)The upper panel represents low solar activity year (F10.7=80), the lower panel represents high solar activity year(F10.7=120); from left to right is spring equinox (Doy=80),summer solstice(Doy=173),autumn equinox(Doy=266) and wintersolstice(Doy=356).

(5) NPPHM能夠較好的模擬hmF2的日變化、季節(jié)變化及赤道異常等特性.

本文采用的非線性多項(xiàng)式模型也可以用于對(duì)電離層電子總含量(TEC)及F2層峰值密度(NmF2)的建模中.今后可以進(jìn)一步補(bǔ)充建模原始數(shù)據(jù)，尤其是太陽(yáng)活動(dòng)高年、高緯地區(qū)和來(lái)自赤道地區(qū)測(cè)高儀的數(shù)據(jù)，以使模型更加完善.

致謝感謝German Aerospace Center(DLR) Neustrelitz提供CHAMP和GRACE衛(wèi)星電離層探測(cè)數(shù)據(jù)，感謝CDAAC(COSMIC Data Analysis and Archive Center)提供COSMIC衛(wèi)星電離層探測(cè)數(shù)據(jù)，感謝SPIDR(Space Physics Interactive Data Resource)提供地面測(cè)高儀探測(cè)數(shù)據(jù).感謝國(guó)家自然科學(xué)基金(40505005)的支持.

Altadill D, Magdaleno S, Torta J M, et al. 2013. Global empirical models of the density peak height and of the equivalent scale height for quiet conditions.Adv.SpaceRes., 52(10): 1756-1769. Bilitza D, Eyfrig R, Sheikh N M. 1979. A global model for the height of the F2-peak using M3000 values from the CCIR numerical map.ITUTelecommun.J., 46: 549-553.

Bilitza D. 1990. International reference ionosphere 1990, NSSDC 22-90. Greenbelt, Maryland, USA: National Space Science Data Center.

Bilitza D, Reinisch B W. 2008. International reference ionosphere 2007: Improvements and new parameters.Adv.SpaceRes., 42(4): 599-609.

Bilitza D. 2015. The international reference ionosphere-Status 2013.Adv.SpaceRes., 55(8): 1914-1927.

Brunini C, Conte J F, Azpilicueta F, et al. 2013. A different method to update monthly medianhmF2values.Adv.SpaceRes., 51(12): 2322-2332.

Garcia-Fernandez M, Hernandez-Pajares M, Juan M, et al. 2003. Improvement of ionospheric electron density estimation with GPSMET occultations using Abel inversion and VTEC information.J.Geophys.Res., 108(A9): 1338. Hajj G A, Romans L J. 1998. Ionospheric electron density profiles obtained with the global positioning system: Results from the GPS/MET experiment.RadioSci., 33(1): 175-190.

Hoque M M, Jakowski N. 2012. A new global model for the ionospheric F2 peak height for radio wave propagation.Ann.Geophys., 30(5): 797-809.

Kakinami Y, Watanabe S, Oyama K I. 2008. An empirical model of electron density in low latitude at 600 km obtained by Hinotori satellite.Adv.SpaceRes., 41(9): 1495-1499.

Kakinami Y, Chen C H, Liu J Y, et al. 2009. Empirical models of total electron content based on functional fitting over Taiwan during geomagnetic quiet condition.Ann.Geophys., 27(8): 3321-3333.

Krankowski A, Zakharenkova I, Krypiak-Gregorczyk A, et al. 2011. Ionospheric electron density observed by FORMOSAT-3/COSMIC over the European region and validated by ionosonde data.J.Geod., 85(12): 949-964.

Kutiev I S, Marinov P G, Watanabe S. 2006. Model of topside ionosphere scale height based on topside sounder data.Adv.SpaceRes., 37(5): 943-950.

Lei J H, Syndergaard S, Burns A G, et al. 2007. Comparison of COSMIC ionospheric measurements with ground-based observations and model predictions: Preliminary results.J.Geophys.Res., 112(A7): A07308.Li K Y. 1986. A new similarity parameter and its application.ActaMeteor.Sinica(in Chinese), 44(2): 174-183.

Lin J, Yue X A, Zeng Z, et al. 2014. Empirical orthogonal function analysis and modeling of the ionospheric peak height during the years 2002-2011.J.Geophys.Res., 119(5): 3915-3929. Liu H, Lühr H. 2005. Strong disturbance of the upper thermospheric density due to magnetic storms: CHAMP observations.J.Geophys.Res., 110(A9): A09S29. Liu H X, Hirano T, Watanabe S. 2013. Empirical model of the thermospheric mass density based on CHAMP satellite observations.J.Geophys.Res., 118(2): 843-848.

Magdaleno S, Altadill D, Herraiz M, et al. 2011. Ionospheric peak height behavior for low, middle and high latitudes: A potential empirical model for quiet conditions-comparison with the IRI-2007 model.J.Atoms.SolarTerr.Phys., 73(13): 1810-1817. Müller S, Lühr H, Rentz S. 2009. Solar and magnetospheric forcing of the low latitude thermospheric mass density as observed by CHAMP.Ann.Geophys., 27(5): 2087-2099.

Radicella S M, Leitinger R. 2001. The evolution of the DGR approach to model electron density profiles.Adv.SpaceRes., 27(1): 35-40.

Rawer K, Ramakrishnan S, Bilitza D. 1978. International Reference Ionosphere, International Union of Radio Science. Special Report. Brussels, Belgium. Shubin V N, Karpachev A T, Tsybulya K G. 2013. Global model of the F2layer peak height for low solar activity based on GPS radio-occultation data.J.Atoms.SolarTerr.Phys., 104: 106-115. Shubin V N. 2015. Global median model of the F2-layer peak height based on ionospheric radio-occultation and ground-based Digisonde observations.Adv.SpaceRes., 56(5): 916-928.Weng L B, Fang H X, Xie Y Q, et al. 2011. Application of analog prediction method on ionospheric TEC short-term forecast.ChineseJ.SpaceSci. (in Chinese), 31(6): 747-753.

Xiong N L, Tang C C, Li X J. 1999. Introduction of the Ionospheric Physics (in Chinese). Wuhan: Wuhan University Press.

Yue X, Schreiner W S, Lei J, et al. 2010. Error analysis of Abel retrieved electron density profiles from radio occultation measurements.Ann.Geophys., 28(1): 217-222.

Yue X A, Schreiner W S, Kuo Y H, et al. 2015. Ionosphere equatorial ionization anomaly observed by GPS radio occultations during 2006-2014.J.Atoms.SolarTerr.Phys., 129: 30-40.

Zhang M L, Liu C, Wan W, et al. 2009. A global model of the ionospheric F2peak height based on EOF analysis.Ann.Geophys., 27(8): 3203-3212.

Zhang M L, Liu C X, Wan W X, et al. 2010. Evaluation of global modeling ofM(3000) F2andhmF2based on alternative empirical orthogonal function expansions.Adv.SpaceRes., 46(8): 1024-1031.

附中文參考文獻(xiàn)李開(kāi)樂(lè). 1986. 相似離度及其使用技術(shù). 氣象學(xué)報(bào), 44(2): 174-183. 翁利斌, 方涵先, 解妍瓊等. 2011. 相似預(yù)報(bào)法在電離層TEC短期預(yù)報(bào)中的應(yīng)用. 空間科學(xué)學(xué)報(bào), 31(6): 747-753.

熊年祿, 唐存琛, 李行健. 1999. 電離層物理概論. 武漢: 武漢大學(xué)出版社.

(本文編輯汪海英)

Global model of ionospheric hmF2based on CHAMPE, GRACE and COSMIC radio occultation

LIU Zhen-Di1， FANG Han-Xian1*， WENG Li-Bin1， MA Qiang1， ZHANG Jian-Bin2

1InstituteofMeteorologyandOceanography,PLAUniversityofScienceandTechnology,Nanjing211101,China2UnitNo. 63655ofPLA,Urumqi841700,China

For global ionospherichmF2modeling, a nonlinear polynomial model approach based on globalhmF2observational data form ionospheric radio occultation (IRO) measurements onboard CHAMP, GRACE, and COSMIC satellites is presented in this paper. The Nonlinear Polynomial Peak Height Model (NPPHM) is constructed by a nonlinear fit withhmF2measurements in least squares sense and describes the dependencies ofhmF2on geomagnetic activity, solar activity, geographical longitude, local time, day of year and geographical latitude. Using independent data from CHAMP satellite, quantitative analysis is made. The correlation coefficients for proposed model NPPHM and CHAMP data are 0.798 in 2008 and 0.732 in 2012, respectively. The corresponding coefficients for IRI2012 are 0.532 and 0.488. NPPHM shows root mean squared errors (RMS) of 25.97 km and 31.39 km in 2008 and 2012, respectively. The corresponding values for IRI2012 are 44.56 km and 42.83 km. Analog deviations are calculated to compare the NPPHM with IRI2012, using data from 14 different world wide ionosonde stations. The deviation of NPPHM is much less than that of IRI2012. Using observational data from the Athens station from 2003 to 2013, the mean deviation of IRI2012 is 8.11%, more than 3.53% of NPPHM. In low solar activity years and every October during the 11 years, the accuracy of NPPHM is much higher than that of IRI2012. What’s more, NPPHM can well present the daily variation, seasonal variation and equatorial ionization anomaly phenomenon ofhmF2.

IonospherichmF2modeling； Nonlinear polynomial； Radio occultation data

10.6038/cjg20161003.

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

劉楨迪,男,1991年生,碩士研究生,研究方向?yàn)殡婋x層物理.E-mail:lzdion@163.com

方涵先,男,1974年生,教授,博士生導(dǎo)師,研究方向?yàn)殡婋x層物理.E-mail:fanghxp@163.com

10.6038/cjg20161003

P352

2015-10-09，2016-01-25收修定稿

劉楨迪， 方涵先， 翁利斌等. 2016. 基于CHAMP、GRACE和COSMIC掩星數(shù)據(jù)的全球電離層hmF2建模研究. 地球物理學(xué)報(bào)，59(10):3555-3565,

Liu Z D， Fang H X，Weng L B,et al. 2016. Global model of ionospherichmF2based on CHAMPE, GRACE and COSMIC radio occultation.ChineseJ.Geophys. (in Chinese),59(10):3555-3565,doi:10.6038/cjg20161003.