王 靜，唐 義，張止戈，鄭旭麗，倪國(guó)強(qiáng)

北京理工大學(xué)光電學(xué)院 光電成像技術(shù)與系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室，北京 100081

1 引 言

電離層是近地衛(wèi)星運(yùn)行的主要場(chǎng)所，也是空間天氣對(duì)人類活動(dòng)影響的重要區(qū)域．由太陽(yáng)風(fēng)暴引起的電離層擾動(dòng)會(huì)影響GPS定位系統(tǒng)的精度、無(wú)線電通信的質(zhì)量和電力傳輸?shù)陌踩婋x層也是導(dǎo)彈、低軌衛(wèi)星和空間站的主要運(yùn)行場(chǎng)所，其環(huán)境狀況將直接影響運(yùn)行其中的飛行器壽命、功能實(shí)現(xiàn)以及宇航員的健康安全等［1］．由于電離層擾動(dòng)時(shí)常發(fā)生，變化快，動(dòng)態(tài)范圍大，因此有效地監(jiān)測(cè)電離層狀態(tài)，尤其是電子密度剖面（Electron Density Profile，EDP）和總電子含量（Total Electron Content，TEC）成為國(guó)內(nèi)外研究的重點(diǎn)．

傳統(tǒng)的探測(cè)方法，如電離層測(cè)高儀只能探測(cè)底部電離層，非相干散射雷達(dá)能探測(cè)800km 以下電子密度剖面，提供電離層多個(gè)參量信息，但其造價(jià)和運(yùn)行成本非常昂貴，且二者都受到站點(diǎn)分布、探測(cè)時(shí)間和頻率的限制，不能滿足全球高分辨力探測(cè)的需求．自20世紀(jì)70年代以來(lái)，美國(guó)NASA 開(kāi)始發(fā)展用于電離層極光和氣輝觀測(cè)的遠(yuǎn)紫外光譜成像遙感技術(shù)，因其良好的時(shí)空分辨力而備受該領(lǐng)域的研究者的青睞．過(guò)去，針對(duì)如何由遠(yuǎn)紫外成像光譜儀觀測(cè)數(shù)據(jù)反演得到電離層參量，國(guó)內(nèi)外學(xué)者進(jìn)行了大量研究［2－7］．在遠(yuǎn)紫外氣輝遙感EDP的反演問(wèn)題上，國(guó)外的研究大多是在電子密度剖面滿足Chapman函數(shù)的假設(shè)下進(jìn)行的，這使得到的峰值電子密度與峰值高度誤差較大；而國(guó)內(nèi)對(duì)這一問(wèn)題的研究尚未涉足．研究遠(yuǎn)紫外光譜遙感反演EDP 的方法，不但能夠豐富電離層監(jiān)測(cè)的數(shù)據(jù)來(lái)源，也將為未來(lái)我國(guó)同類光譜數(shù)據(jù)快速轉(zhuǎn)化為數(shù)據(jù)產(chǎn)品服務(wù)．

本研究基于搭載TIMED 飛船的全球紫外成像儀（Global Ultraviolet Imager，GUVI）探測(cè)得到的氧原子135．6nm 夜氣輝光譜數(shù)據(jù)，利用正則化和牛頓迭代法相結(jié)合的方法計(jì)算得到了自由形式的電離層電子密度剖面，之后將其與電離層電子密度的Chapman表達(dá)式相結(jié)合，得到Chapman擬合結(jié)果．

2 數(shù)據(jù)來(lái)源及觀測(cè)模型

2．1 數(shù)據(jù)來(lái)源

GUVI［8］于2001年12月搭載TIMED衛(wèi)星進(jìn)入近極地太陽(yáng)同步軌道，軌道周期為97．8 min，軌道高度625km，傾角為74．1°．本研究利用GUVI臨邊觀測(cè)得到的氧原子135．6nm 夜氣輝光譜數(shù)據(jù)反演得到電子密度剖面，其體輻射速率可表示為

式中，εmn表征了中和作用在氧原子135．6nm 夜氣輝中的貢獻(xiàn)，R1＝7．3×10－13、R2＝10－7、R3＝1．3×10－15、R4＝1．4×10－10分別是相關(guān)反應(yīng)的反應(yīng)速率，單位為cm3/s．

如圖1所示，是利用大氣模式NRL MSISE－00和國(guó)際參考電離層模型IRI－2011計(jì)算得到的1992年3月21日（30°N，0°E）處中和作用的貢獻(xiàn)εmn．氧原子135．6nm 夜氣輝中，約75%來(lái)源于輻射復(fù)合反 應(yīng)，25%來(lái)源于中和作用［9－11］，且由圖可知中和作用的貢獻(xiàn)僅在250km 以下才變得顯著．因此，在計(jì)算中可以忽略中和作用的影響，認(rèn)為氧原子135．6nm夜氣輝光譜數(shù)據(jù)僅由輻射復(fù)合反應(yīng)產(chǎn)生．另外，由于O＋是F層等離子體中最主要的成分（含量在99%以上），且F層呈整體電中性，在本研究中認(rèn)為

圖1 εmn 隨海拔高度和當(dāng)?shù)貢r(shí)間的變化Fig．1 εmnvariation with altitude and local time

因此，體輻射速率η可近似表示成

其中，［e］為電子密度．

2．2 觀測(cè)模型

在TIMED 飛船沿軌飛行的同時(shí)，GUVI掃描成像系統(tǒng)在穿軌方向以0．4°的步長(zhǎng)掃描與星下點(diǎn)成67．2°到80°的臨邊區(qū)域，在一次掃描中獲得32×14個(gè)像元值．如圖2所示，從儀器的成像方式來(lái)看，每個(gè)像元值是觀測(cè)方向上從電離層切點(diǎn)到GUVI所在位置的輻射之和．為了去除隨機(jī)噪聲帶來(lái)的影響，對(duì)沿軌方向的14個(gè)像元求平均賦值于中心像元．于是，一次掃描的32個(gè)像元與32個(gè)切點(diǎn)高度相關(guān)．若將電離層看作是球?qū)ΨQ分層結(jié)構(gòu)，則海拔高度為90～530km 的電離層均勻分成23層，可建立離散形式的觀測(cè)模型對(duì)一次掃描來(lái)說(shuō)，觀測(cè)向量B為32×1的列向量；權(quán)重矩陣W為32×23，其每行中各元素為各層輻射在觀測(cè)值中所占比例；體輻射速率η為23×1的列向量；N為附加的噪聲項(xiàng)．

圖2 GUVI臨邊觀測(cè)結(jié)構(gòu)示意圖Fig．2 Geometrical chart of GUVI limb observing

2．2．1 噪聲分析

光子散粒噪聲是GUVI探測(cè)數(shù)據(jù)噪聲的主要來(lái)源，本研究?jī)H考慮光子散粒噪聲的影響．光子散粒噪聲是與入射光子數(shù)的隨機(jī)變化有關(guān)的噪聲，在給定時(shí)間內(nèi)入射光子數(shù)服從泊松分布，光子散粒噪聲σshot可由信號(hào)電子數(shù)Se（以電子數(shù)計(jì)）的均方根值表示為

GUVI靈敏度為0．5counts/s/Rayleigh/pixel，像元積分時(shí)間為0．064s，考慮本模型采取了14像元合并的方式，則模型中附加噪聲項(xiàng)可表示為

其中，B取式（5）觀測(cè)向量中的觀測(cè)值．

2．2．2 權(quán)重矩陣

前面，在電離層為球?qū)ΨQ分層結(jié)構(gòu)這一假設(shè)下建立了觀測(cè)模型．權(quán)重矩陣每行中各個(gè)元素僅與各層到探測(cè)器的距離有關(guān)［12］．于是，在觀測(cè)模型中，可認(rèn)為所有位置的權(quán)重矩陣是相同的．權(quán)重矩陣的特性直接影響反演結(jié)果的準(zhǔn)確性，本研究計(jì)算得到法矩陣WTW的條件數(shù)為3469．一般認(rèn)為條件數(shù)小于100 時(shí)為良態(tài)；條件數(shù)在100到1000之間時(shí)為中等程度的病態(tài)；條件數(shù)超過(guò)1000 時(shí)存在嚴(yán)重的病態(tài)［13］．可見(jiàn)本研究中電子密度剖面反演問(wèn)題存在嚴(yán)重病態(tài)．

綜上，由建立的臨邊觀測(cè)模型得出電離層電子密度剖面反演問(wèn)題轉(zhuǎn)化為如下的最小二乘最優(yōu)問(wèn)題：

其中，b為觀測(cè)值與噪聲值之差．然而，由于觀測(cè)模型中權(quán)重矩陣W的病態(tài)性，往往使得由普通的最小二乘法求得的估計(jì)值與真值之間存在較大差別，因此，本研究采用正則化方法增加約束項(xiàng)來(lái)求解上述問(wèn)題．

3 反演方法

針對(duì)上述病態(tài)問(wèn)題，Tikhonov提出了正則化方法．正則化方法的關(guān)鍵是增加全部或部分參數(shù)（或參數(shù)改正數(shù)）加權(quán)平方和極小的條件來(lái)增加約束，克服不適定性，使解唯一且穩(wěn)定，即用相鄰的適定問(wèn)題的解去逼近原問(wèn)題的解．于是，為使式（8）有唯一穩(wěn)定的解，構(gòu)造準(zhǔn)則函數(shù)

使式（9）最小化的參數(shù)η即為所建立的觀測(cè)模型的正則化解．其中，Ω（η）稱為穩(wěn)定泛函，α是正則化參數(shù)（或稱為平滑因子、嶺參數(shù)）．在實(shí)際計(jì)算中，穩(wěn)定泛函Ω（η）可取不同的形式，這里取Ω（η）＝ηTHη，其中正則化矩陣H＝DT2D2，D2通?？扇挝痪仃?、拉普拉斯算子的一階和二階差分矩陣，這里取二階差分形式的拉普拉斯算子，即

目前，確定正則化參數(shù)的常用方法有L－曲線法和GCV 法（Generalized Cross－Validation）［14］，本研 究采用L－曲線法確定正則化參數(shù)．由于‖Wη－b‖2和‖D2η‖2都是正則化參數(shù)α的函數(shù)，選擇不同的α值，以‖Wη－b‖2為橫坐標(biāo)，‖D2η‖2為縱坐標(biāo)，經(jīng)過(guò)曲線擬合得到L－曲線，如圖3所示．

圖3 L－曲線Fig．3 L－curve

由于牛頓迭代法具有快速收斂的特點(diǎn)，本研究采用其尋找未知參數(shù)η的最優(yōu)解．迭代公式［15］如下

其中，J為Wη在ηk處的雅克比矩陣取值，C為觀測(cè)值噪聲的協(xié)方差矩陣．迭代終止條件為

其中，B0為輻射估計(jì)值．

4 反演結(jié)果及分析

4．1 與地基探測(cè)數(shù)據(jù)的對(duì)比

利用正則化和牛頓迭代相結(jié)合的方法，選取存在近同步傳統(tǒng)觀測(cè)的GUVI臨邊觀測(cè)氧原子135．6nm夜氣輝輻射數(shù)據(jù)，反演得到2002年10月4日和2004年7月2日自由形式的電離層電子密度剖面，詳細(xì)參數(shù)列于表1中，將地基近同步觀測(cè)數(shù)據(jù)及全解析性電離層模型FAIM［16］得到的電子密度剖面與本研究結(jié)果繪于圖4中．其中，2002年10月4日地基觀測(cè)數(shù)據(jù)來(lái)源于Millstone Hill非相干散射雷達(dá)（Incoherent Scatter Radar，ISR），Athens數(shù)字測(cè)高儀臺(tái)站只提供2004年7月2日的峰值電子密度和峰值高度，其地基觀測(cè)結(jié)果是由Chapman函數(shù)擬合得到的．

表1 位置參數(shù)列表Table 1 Parameters of selected locations

如圖4a 和4c所示，與全解析性電離層模型FAIM 相比，本研究得到的結(jié)果在海拔200km 以上與地基方式探測(cè)結(jié)果更為接近，然而，在海拔高度200km 以下，卻與其結(jié)果相差甚遠(yuǎn)．造成此現(xiàn)象的原因可能有兩個(gè)：一是海拔高度200km以下，未考慮中和作用在氧原子135．6nm 夜氣輝輻射中的影響；二是GUVI觀測(cè)資料和地基探測(cè)結(jié)果在空間（經(jīng)度上相差7°左右）和時(shí)間上并非完全一致．盡管海拔200km 以下本研究反演得到的電子密度偏差很大，但是其峰值電子密度和峰值高度吻合得很好．

一般情況下，可用Chapman函數(shù)來(lái)描述電離層中的電子密度［17］

其中，z為海拔高度，Nm和hm分別為峰值電子密度和峰值高度，標(biāo)高H取值為54km．將本研究計(jì)算得到的峰值電子密度和峰值高度代入式（13）中，得到Chapman型電子密度剖面，如圖4b 所示．由圖4b和4d可見(jiàn)，利用本研究提出的正則化解法得到峰值電子密度和峰值高度，結(jié)合Chapman型表達(dá)式得到的電子密度剖面與地基方式探測(cè)結(jié)果符合得很好．

綜上所述，利用正則化與牛頓迭代相結(jié)合的方法由遠(yuǎn)紫外夜氣輝OI 135．6nm 光譜數(shù)據(jù)反演得到的自由形式EDP，與全解析性電離層模型FAIM 結(jié)果相比，具有明顯的優(yōu)越性，但是其與地基觀測(cè)結(jié)果，尤其是在海拔200km 以下仍存在一定差異．盡管如此，本研究得到的峰值電子密度和峰值高度與地基觀測(cè)結(jié)果吻合得很好．因此，利用本研究方法得到的電子密度剖面代替?zhèn)鹘y(tǒng)電離層探測(cè)方法，彌補(bǔ)了非相干散射雷達(dá)和數(shù)字測(cè)高儀臺(tái)站稀少，難以實(shí)時(shí)獲得高空間分辨力全球電離層電子密度結(jié)構(gòu)和特性的不足．

4．2NmF2隨磁暴變化分析

利用本研究提出的算法反演了2002年9月29日到10月3日GUVI探測(cè)得到的5°N 處的夜氣輝數(shù)據(jù)，初步分析電離層電子密度隨磁暴的變化情況．國(guó)際上，一般采用Dst指數(shù)來(lái)描述磁暴：－50＜Dst≤－30 為小磁暴，－100＜Dst≤－50 為中等磁暴，－200＜Dst≤－100為大磁暴，Dst≤－200為特大磁暴．圖5中繪出了2002年9月29日到10月3日的Dst（來(lái)源于美國(guó)國(guó)家海洋大氣局NOAA）變化曲線．從圖5可以看出，9月29日處于磁寧?kù)o期，9月30日發(fā)生了小磁暴，10月1日發(fā)生了大磁暴，10月2日到10 月3 日為磁暴的恢復(fù)相，Dst指數(shù)開(kāi)始回升．

圖5 2002年9月29日—10月3日Dst指數(shù)變化Fig．5 Dst index for Sep 29—Oct 3，2002

如圖6所示為反演得到的2002年9月29日到10月3日EDPs，在圖像橫軸下方標(biāo)注數(shù)據(jù)的確切經(jīng)度和獲取時(shí)間，當(dāng)存在軌道缺失時(shí)采用鄰近插值得到所需數(shù)據(jù)且不標(biāo)記其經(jīng)度和獲取時(shí)間．在此，以磁寧?kù)o期9月29日（圖6a）為參考，初步分析磁暴對(duì)電離層電子密度的影響情況．

9月30日凌晨4時(shí)起Dst指數(shù)顯著下降，持續(xù)到10時(shí)在－20nT 附近上下波動(dòng)．如圖6b所示，經(jīng)度208°～257°附近出現(xiàn)電離層赤道異常（Equatorial Anomaly）峰值高度處的電子密度明顯大于9月29日相應(yīng)位置值；13∶00到20∶00發(fā)生了小磁暴，相應(yīng)位置處的電子密度較9月29日也有所提升．

10月1日3時(shí)至10時(shí)Dst指數(shù)在－20～－40nT范圍內(nèi)波動(dòng)，西半球經(jīng)度201°～299°附近電離層赤道異常區(qū)范圍顯著變大（圖6c）；之后Dst指數(shù)迅速下降，并在17時(shí)左右達(dá)到最小值，Dst指數(shù)在－160nT上下起伏并一直持續(xù)到午夜，電離層F 層的電子密度普遍變大，且向上抬升．

凌晨到下午16時(shí)10月2 日的Dst指數(shù)比10月1日的更小，因此在2002年10月2日（圖6d）西半球經(jīng)度242°～291°內(nèi)的電離層赤道異常范圍較10月1日進(jìn)一步擴(kuò)大，然而在193°～242°內(nèi)的電子密度卻有所下降．16時(shí)以后Dst指數(shù)開(kāi)始回升，峰值高度或峰值電子密度出現(xiàn)大范圍下降．2002年10月3日Dst指數(shù)繼續(xù)回升，如圖6e所示赤道異常普遍減弱．

圖6 2002年9月29日—10月3日反演結(jié)果Fig．6 Retrieved EDPs for Sep 29—Oct 3，2002

經(jīng)分析2002年9月29日到10月3日5°N 處的EDPs隨磁暴的變化情況發(fā)現(xiàn)，隨著磁暴的發(fā)生及其強(qiáng)弱變化，在峰值高度附近會(huì)出現(xiàn)電離層赤道異常，其變化程度和影響范圍在各經(jīng)度處不盡相同．

5 結(jié) 論

本研究利用全球遠(yuǎn)紫外成像光譜儀GUVI獲得的氧原子135．6nm 夜氣輝光譜輻射數(shù)據(jù)，建立了臨邊觀測(cè)模型，采用正則化與牛頓迭代法相結(jié)合的方法，計(jì)算得到了自由形式電離層電子密度；將自由形式電離層電子密度的峰值高度和峰值電子密度與Chapman型表達(dá)式相結(jié)合，得到了Chapman擬合電子密度，并將其與地基觀測(cè)數(shù)據(jù)及全解析性電離層模型FAIM 進(jìn)行了比較．結(jié)果表明，在海拔高度200km 以上本研究得到的電子密度與地基觀測(cè)數(shù)據(jù)符合得很好，但在200km 以下偏差較大，產(chǎn)生這一現(xiàn)象的原因有待進(jìn)一步研究．

利用本研究提出的算法反演得到了2002年9月29日到10月3日5°N 處的EDPs，初步分析電離層電子密度隨磁暴的變化情況發(fā)現(xiàn)，隨著磁暴的發(fā)生及其強(qiáng)弱變化，在峰值高度附近會(huì)出現(xiàn)電離層赤道異常，其變化程度和影響范圍在各經(jīng)度處不盡相同．致 謝 感謝美國(guó)霍普金斯大學(xué)和美國(guó)TIMED/GUVI衛(wèi)星的數(shù)據(jù)支持．

（References）

［1］ 吳健，郭兼善．空間天氣對(duì)通信和導(dǎo)航定位等無(wú)線電系統(tǒng)有什么影響？全球定位系統(tǒng)，1999，24（3－4）：1－4．

Wu J，Guo J S．What effects does space weather have on radio systems？GNSSWorldofChina（in Chinese），1999，24（3－4）：1－4．

［2］ Dymond K F，Thonnard S E，McCoy R P，et al．An optical remote sensing technique for determining nighttime F region electron density．RadioSci．，1997，32（5）：1985－1996．

［3］ Rajesh P K，Liu J Y，Hsu M L，et al．Ionospheric electron content and NmF2from nighttime OI 135．6nm intensity．J．Geophys．Res．，2011，116（A02313）：1－16．

［4］ Paxton L J，Morrison D，Strickland D J，et al．The use of far ultraviolet remote sensing to monitor space weather．Adv．SpaceRes．，2003，31（4）：813－818．

［5］ Wang J，Tang Y，Tang L J，et al．Retrieval of ionospheric O/N2based on FUV imaging data．SPIE，2009，7498：74982N－1－74982N－9．

［6］ 彭圣鋒，唐義，王靜等．遠(yuǎn)紫外光譜遙感反演熱層O/N2技術(shù)研究．光譜學(xué)與光譜分析，2012，32（5）：1296－1300．

Peng S F，Tang Y，Wang J，et al．Research on retrieving thermospheric O/N2from FUV remote sensing．Spectroscopy andSpectralAnalysis（in Chinese），2012，32（5）：1296－1300．

［7］ 王詠梅，付利平，王英鑒．星載遠(yuǎn)紫外極光/氣輝探測(cè)發(fā)展綜述．地球物理學(xué)進(jìn)展，2008，23（5）：1474－1479．

Wang Y M，F(xiàn)u L P，Wang Y J．Review of space－based FUV aurora/airglow observations．ProgressinGeophysics（in Chinese），2008，23（5）：1474－1479．

［8］ Paxton L J，Christensen A B，Morrison D，et al．GUVI：a hyperspectral imager for geospace．SPIE，2004，5660：228－240．

［9］ Jeffrey S D．Inversion of nighttime electron densities using satellite observations［Ph．D．thesis］．New Mexico：New Mexico Institute of Mining and Technology，2001．

［10］ Hanson W B．Radiative recombination of atomic oxygen ions in the nighttime F region．J．Geophys．Res．，1969，74（14）：3720－3722．

［11］ McDonald S E．Day to day and longitudinal variability of the nighttime low latitude terrestrial ionosphere［Ph．D．thesis］．Fairfax：George Mason University，2007．

［12］ DeMajistre R，Paxton L J，Morrison D，et al．Retrievals of nighttime electron density from Thermosphere Ionosphere Mesosphere Energetics and Dynamics （TIMED）mission Global Ultraviolet Imager （GUVI） measurements．J．Geophys．Res．，2004，109（A5）：1－14．

［13］ 吳杰，李明峰，余騰．測(cè)量數(shù)據(jù)處理中病態(tài)矩陣和正則化方法．大地測(cè)量與地球動(dòng)力學(xué)，2010，30（4）：102－105，108．

Wu J，Li M F，Yu T．Ill matrix and regularization method in surveying data processing．JournalofGeodesyandGeodynamics（in Chinese），2010，30（4）：102－105，108．

［14］ 王振杰，歐吉坤．用L－曲線法確定嶺估計(jì)中的嶺參數(shù)．武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），2004，29（3）：235－238．

Wang Z J，Ou J K．Determining the ridge parameter in a ridge estimation usingL－curve method．Geomaticsand InformationScienceofWuhanUniversity（in Chinese），2004，29（3）：235－238．

［15］ 蔣德明，董超華．大氣廓線物理反演的最優(yōu)化方法進(jìn)展．地球科學(xué)進(jìn)展，2010，25（2）：133－139．

Jiang D M，Dong C H．A review of optimal algorithm for physical retrieval of atmospheric profiles．AdvancesinEarth Science（in Chinese），2010，25（2）：133－139．

［16］ Solomon S C．Modeling of the thermosphere－ionosphere system．Phys．Chem．Earth，2000，25（5－6）：499－503．

［17］ 梅冰，萬(wàn)衛(wèi)星．基于Millstone Hill非相干散射雷達(dá)觀測(cè)的電離層電子濃度剖面的經(jīng)驗(yàn)正交函數(shù)分析．地球物理學(xué)報(bào)，2008，51（1）：10－16．

Mei B，Wan W X．An empirical orthogonal function（EOF）analysis of ionospheric electron density profiles based on the observation of incoherent scatter radar at Millstone Hill．ChineseJ．Geophys．（in Chinese），2008，51（1）：10－16．