基于氧氣A波段發(fā)射譜線臨近空間大氣溫度的反演及分析

楊曉君，王后茂, 李葉飛，王詠梅，胡秀清

1. 中國科學(xué)院國家空間科學(xué)中心, 空間環(huán)境探測研究室，北京 100190 2. 天基空間環(huán)境探測北京市重點(diǎn)實(shí)驗(yàn)室，北京 100190 3. 中國科學(xué)院空間環(huán)境態(tài)勢感知技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100190 4. 中國科學(xué)院大學(xué)天文與空間科學(xué)學(xué)院，北京 100049 5. 上海衛(wèi)星工程研究所，上海 200240 6. 中國氣象局國家衛(wèi)星氣象中心，北京 100081

引 言

地球大氣臨近空間區(qū)域的物理和化學(xué)特性研究需要準(zhǔn)確的了解該區(qū)域的溫度分布，在臨近空間區(qū)域發(fā)生的許多過程都具有明顯的溫度依賴性。此外，臨近空間大氣參數(shù)很大程度上也影響著空間天氣預(yù)報以及航天器發(fā)射與再入軌過程的安全性預(yù)估等[1]，因此對臨近空間大氣的觀測與研究具有重要的科學(xué)意義和應(yīng)用價值。雖然火箭和地面測量已經(jīng)可以獲得準(zhǔn)確的溫度廓線，但很難用它們建立一個全球的溫度分布數(shù)據(jù)庫，星載探測可以在較短時間內(nèi)進(jìn)行全球或者區(qū)域的大范圍觀測。美國空間物理研究實(shí)驗(yàn)室基于搭載在URAS衛(wèi)星上的高分辨率多普勒成像儀HRDI(The high-resolution Doppler imager)通過觀測氧氣A波段的發(fā)射譜線進(jìn)行溫度的測量[2], 采用基于擾動理論的最優(yōu)估計(jì)法對80～100 km高度的大氣進(jìn)行反演，利用三條譜線協(xié)同反演得到誤差約為7 K的溫度廓線，首次提供了該高度范圍全球溫度測量結(jié)果。之后, 2001年搭載在TIMED衛(wèi)星上的TIDI(TIMED Doppler Interferometer)多普勒成像儀基于氧氣A波段發(fā)射譜線對60～300 km高度的臨近空間大氣溫度進(jìn)行了遙感觀測，基于譜線的觀測值采用約束非線性最小二乘擬合得到溫度廓線[3]。2001年2月20日，搭載在Odin衛(wèi)星上的光譜儀和紅外成像系統(tǒng)OSIRIS(optical spectrograph and infraRed imaging system)對氧氣A帶的發(fā)射光譜進(jìn)行了探測，將90～110 km內(nèi)的氧氣A波段光譜與模擬光譜進(jìn)行逐像素比較，用最小二乘法使波段譜差的平方和最小，從而得到中層-低熱層區(qū)域的溫度，在90 km附近的溫度反演精度約為±2 K，在較高高度的精度為±6 K[4]。先進(jìn)地球觀測衛(wèi)星ADEOS上搭載的大氣光譜儀ILAS測量了包括以762 nm為中心的分子氧A波段在內(nèi)的753～784 nm波段的大氣吸收光譜，通過衛(wèi)星太陽掩星對氧氣A波段吸收光譜的測量，采用層析法獲得了整個平流層的溫度分布[5]，與儀器功能相關(guān)的誤差是系統(tǒng)不確定性中最大的誤差源, 與所考慮的系統(tǒng)不確定性相關(guān)的溫度的均方根誤差估計(jì)為4 K。全球高分辨率熱層成像邁克爾遜干涉儀MIGHTI(Michelson Interferometer for Global High-resolution Thermospheric Imaging)搭載在NASA的電離層連接探測器ICON上于2019年10月11日發(fā)射[6]，MIGHTI被用來測量地球熱層的中性溫度[7]，反演方法是使用離散波長測量A波段的光譜形狀，從波段的旋轉(zhuǎn)包絡(luò)線推斷出環(huán)境溫度，在白天90～105 km之間的測量精度為1K，而在夜間，測量精度在90 km時為1 K，在105 km時增加到了3 K[8]。

以上儀器選取的譜線都不同，而影響反演結(jié)果非常重要的一點(diǎn)就是譜線的選擇。我們在氧氣A帶所有譜線中對比不同譜線特性，尋找譜線選擇的判斷條件，然后選取合適的譜線基于最優(yōu)估計(jì)的方法進(jìn)行溫度的反演，并進(jìn)一步探討了譜線線強(qiáng)對反演結(jié)果的影響。

1 正演模型

產(chǎn)生氧氣A波段的激發(fā)態(tài)氧氣分子是由一系列復(fù)雜的光化學(xué)反應(yīng)產(chǎn)生的，當(dāng)激發(fā)態(tài)氧氣分子自發(fā)輻射產(chǎn)生氧氣A波段時會發(fā)出光子，波數(shù)ν處的特定轉(zhuǎn)動譜線發(fā)出的光子數(shù)可以表示為fDη，D是多普勒線型，η是體發(fā)射率，f是譜線強(qiáng)度，是溫度的函數(shù)

(1)

式(1)中，fref是在參考溫度Tref=296°K時的譜線強(qiáng)度，E′是轉(zhuǎn)動譜線的上能級能量。式(2)是以vc為譜線中心的多普勒線型

(2)

(3)

其中，k是玻爾茲曼常數(shù)，c是光速，m是氧氣分子量。

綜上所述，衛(wèi)星觀測到的光譜亮度可以從切點(diǎn)高度Zt沿視線進(jìn)行路徑積分得到

(4)

其中，s是視線路徑，n為氧氣數(shù)密度，σ為吸收截面。由式(4)可得，光譜輻射主要由氧氣體發(fā)射率和溫度決定，而光譜輻射對溫度的依賴主要是通過自吸收截面等的作用[9]。

為了更真實(shí)的模擬星載測量，引入儀器函數(shù)，因此模擬的儀器計(jì)數(shù)值用式(5)獲得

(5)

式(5)中，I為卷積的歸一化儀器函數(shù)，是半高寬(FWHM)為20 cm-1的高斯函數(shù)。

2 反演模型

2.1 擾動理論

衛(wèi)星的觀測結(jié)果可表示為

C=KT+ε

(6)

式(6)中，C是衛(wèi)星觀測的輻射強(qiáng)度向量；T是待反演的溫度向量；K是描述大氣物理過程正演模型的權(quán)函數(shù)矩陣；ε是測量誤差向量。

首先，分別模擬臨邊輻射強(qiáng)度測量值和用體發(fā)射率廓線、溫度廓線初值計(jì)算得到的正演模型值。模擬過程中的譜線線強(qiáng)、上能級能量等參數(shù)均來自2004年的HITRAN光譜數(shù)據(jù)庫[10]，而大氣分子濃度和溫度廓線數(shù)據(jù)均來自MISIS-00模型[11]，體發(fā)射率廓線是由光化學(xué)模型利用溫度各組分?jǐn)?shù)密度廓線計(jì)算得到[12]。測量值Cmeans和用初始溫度、體發(fā)射率廓線計(jì)算得到的正演模型值C00之間的差值可以用線性積分算子表示

本文筆者將微課應(yīng)用到成人繼續(xù)教育中，將信息化教學(xué)手段融入了傳統(tǒng)教學(xué)中，從而形成了線上線下的混合式教學(xué)模式。這種混合式教學(xué)模式是一種優(yōu)秀的教學(xué)模式，它綜合了MOOC的優(yōu)勢，彌補(bǔ)了MOOC缺乏管理機(jī)制的缺陷，利用現(xiàn)有大量的MOOC資源，降低了微課制作的工作量，突現(xiàn)了“互聯(lián)網(wǎng)+”時代的優(yōu)勢。

(7)

其中z為切點(diǎn)高度(km)； 矩陣中C和T下標(biāo)為視線切點(diǎn)高度(km)； 矩陣K下標(biāo)為所在行視線積分的高度和依次改變溫度值所對應(yīng)的高度； ΔT=Ttrue-T00為真實(shí)溫度值與初始溫度值的差值；Cmeans是輸入溫度為T00+ΔT的模擬值；C00是輸入溫度為T00的模擬值。

以1 km為格點(diǎn)對溫度廓線加入擾動，每改變一個格點(diǎn)就產(chǎn)生一條新的擾動廓線ΔT/T00，51條廓線就可以用ΔT/T00組成的矢量Tf來表示。為了得到更高的精確度，式(7)中的積分與正演模型中的一樣，積分步長都為1 km。將模擬臨邊掃描的所有測量計(jì)數(shù)值作為C矢量，這樣就得到了一個矩陣方程ΔC=KTf用來求解未知量Tf。

2.2 反演算法

理論上，溫度可以直接利用氧氣A波段的觀測輻射光譜通過求解輻射傳輸方程來獲得，然而事實(shí)上求解這一方程是不可能的，這主要是因?yàn)槠浣馐遣晃ㄒ磺也环€(wěn)定的； 同時由于觀測誤差的存在，使得簡單的單步反演法在實(shí)際應(yīng)用時容易導(dǎo)致較大的反演誤差，抗干擾能力比較低?；谪惾~斯理論的最優(yōu)估計(jì)法是一種求解有噪聲反演問題的非常好的手段，用實(shí)時探測數(shù)據(jù)替代非待反演參數(shù)，并利用先驗(yàn)信息對反演結(jié)果進(jìn)行約束和修正，得到待反演參數(shù)最合理的解，提高了反演精度和抗干擾能力。

在2.1擾動理論的基礎(chǔ)上，選擇了基于貝葉斯理論的最優(yōu)估計(jì)法進(jìn)行溫度的反演。假定測量誤差和先驗(yàn)信息誤差服從高斯分布，將求解輻射傳輸方程轉(zhuǎn)化為使概率函數(shù)最大化的最優(yōu)化數(shù)學(xué)問題，求解這個概率問題就可以獲得溫度的最優(yōu)解。

將T看做一個多維隨機(jī)變量，其概率密度函數(shù)為

(8)

條件概率密度函數(shù)可表示為

(9)

其中，Sa是描述T0不確定性的先驗(yàn)信息誤差協(xié)方差矩陣； 對角矩陣Sε是從觀測值得到的測量噪聲協(xié)方差矩陣，是描述C和正演模型不確定性的觀測誤差協(xié)方差矩陣。假設(shè)測量誤差的總體方差和探測器上的計(jì)數(shù)值一樣服從泊松分布。先驗(yàn)估計(jì)協(xié)方差矩陣Sa可以模擬為

Sa=αTJ

(10)

式(10)中，z是反演的網(wǎng)格高度，αT=103K2。先驗(yàn)估計(jì)協(xié)方差矩陣表示了一個預(yù)期的均方根變化約為5%(～10 K)的先驗(yàn)溫度廓線。

根據(jù)貝葉斯定理，最有可能的大氣狀態(tài)剖面是使式(11)最大化[13]

P(Tf|C)=P(Tf)P(C|Tf)/P(C)

(11)

其中，P(Tf)是溫度參數(shù)分?jǐn)?shù)變化的先驗(yàn)概率分布；P(C)觀測數(shù)據(jù)的先驗(yàn)概率分布，其為常數(shù)；P(C|Tf)給定溫度參數(shù)分?jǐn)?shù)變化條件下具體觀測數(shù)據(jù)的似然概率；P(Tf|C)是組合先驗(yàn)信息和似然概率得到的溫度參數(shù)分?jǐn)?shù)變化后驗(yàn)概率。

將式(8)和式(9)相乘，最有可能的T值是質(zhì)式(12)最小化

(12)

式(12)最小化時的T的估計(jì)值為

(13)

其中，T0=0, 所以方程可化為

(14)

再由Test求得溫度廓線的最優(yōu)估計(jì)值Ttrue。

3 反演結(jié)果及討論

3.1 正演模型

為了驗(yàn)證上述反演方法的可行性和精確度，利用構(gòu)建的正演模型來模擬實(shí)測值進(jìn)行溫度的反演研究。首先利用一組溫度輸入值T00作為初始值來獲得一組正演模擬值，然后在T00的基礎(chǔ)上加上一組隨機(jī)溫度變量得到新的溫度廓線Ttrue=T00+ΔT，再由新的溫度廓線模擬一組實(shí)測輻射值，最后由該兩組輻射模擬測量值組合進(jìn)行溫度反演，將反演得到的結(jié)果與初始溫度輸入值進(jìn)行比較。圖1就是正演模型所用的溫度廓線初始值(T00)和模擬的實(shí)際測量值所用的溫度廓線(即溫度真實(shí)值：Ttrue)。

圖1 溫度廓線初始值和實(shí)際測量值對比

3.2 無噪聲時的反演結(jié)果分析

基于HITRAN數(shù)據(jù)庫中氧氣A帶(759～767 nm)內(nèi)包含的幾十條譜線參數(shù)，我們將所有的譜線進(jìn)行了輻射模擬與溫度反演分析，結(jié)果表明不同譜線的反演結(jié)果精度不同，精度較好的譜線主要集中在759～760和764～766 nm這兩個波段范圍，圖2為氧氣A帶中無噪聲情況下60～110 km范圍內(nèi)的溫度反演精度較高的譜線。如圖2所示，ΔT=Ttrue-T00為真實(shí)溫度值與初始溫度值的差值，其他折線為利用不同譜線反演ΔT的結(jié)果，17條氧氣A波段的溫度反演結(jié)果在60～110 km高度范圍內(nèi)效果較好，平均反演偏差為4.1 K。

圖2 無噪聲條件下譜線在60～110 km范圍內(nèi)的 溫度反演結(jié)果對比

但有些譜線的反演誤差在80 km以上較小，80 km以下卻較大，如圖3所示，764.28和764.17 nm為圖2中反演精度較高的譜線，另外兩條為反演精度較差的譜線(761.72和761.25 nm)。將他們的反演結(jié)果進(jìn)行進(jìn)一步的對比分析，如圖3所示，80 km以上四條譜線顯示了較好的反演精度，平均反演誤差<5 K； 而在80 km以下，761.72和761.25 nm的反演結(jié)果與真實(shí)值存在較大的偏差，平均偏差達(dá)到了34.9 K。原因如下： 由式(4)可以看出，溫度通過影響線強(qiáng)和自吸收兩部分來影響輻射強(qiáng)度，且溫度對它們的影響正好相反，權(quán)函數(shù)就是用來表示溫度對輻射強(qiáng)度影響大小的，而反演結(jié)果的差距可從其權(quán)函數(shù)中得到規(guī)律。圖4為兩條譜線761.25和764.28 nm的權(quán)函數(shù)，為方便對比各選取一條代表線標(biāo)記為紅色(71和80 km)，圖中每一點(diǎn)代表的都是對應(yīng)高度處的溫度變動對衛(wèi)星輻射強(qiáng)度觀測的影響與溫度的擾動量的比值，即ΔC/ΔT，溫度擾動ΔT設(shè)置為5 K，由圖4(a)可以看出，當(dāng)溫度對自吸收的影響所占比重大于對線強(qiáng)的影響時，ΔC/ΔT會發(fā)生正負(fù)翻轉(zhuǎn)，反演精度會變差； 從圖4(b)可以看出，影響因素占主導(dǎo)地位的是線強(qiáng)，溫度主要通過對線強(qiáng)的影響來改變輻射強(qiáng)度，其對自吸收的影響比重較小，ΔC/ΔT未發(fā)生正負(fù)翻轉(zhuǎn)，反演效果比較好。同時，我們將其他譜線的反演結(jié)果進(jìn)行比較得出相同的結(jié)論，因此可以根據(jù)這個規(guī)律，從氧氣A波段眾多譜線中篩選適合用于溫度反演的譜線，即溫度對自吸收影響相對較小的譜線。

圖3 無噪聲條件下四條譜線在60～110 km范圍內(nèi)的 溫度反演結(jié)果對比

圖4 權(quán)函數(shù)矩陣Fig.4 Weight function matrix

3.3 添加噪聲后的反演結(jié)果分析

在實(shí)際的觀測中，除了上述考慮的因素外，還需要考慮信噪比的影響。因此，本文將圖3中的四條譜線模擬值加入隨機(jī)噪聲并進(jìn)行了反演及結(jié)果分析，結(jié)果如圖5所示。764.28和764.17 nm兩條譜線的反演結(jié)果與真實(shí)值基本吻合，說明噪聲對這兩條譜線的影響非常小，該譜線的輻射強(qiáng)度能夠滿足溫度反演的要求。而761.72和761.25 nm這兩條譜線在80 km以下的反演偏差變得更大，最大偏差達(dá)到幾百K，原因是這兩條譜線本身在80 km以下反演誤差就非常大，微小擾動就會造成較大的誤差，所以加入噪聲后使得反演偏差變得更大。

圖5 加入噪聲情況下四條譜線在60～110 km范圍內(nèi)的溫度反演誤差

為了進(jìn)一步分析噪聲對反演結(jié)果的影響，本文選擇了四條較強(qiáng)譜線和一條弱譜線，如表1所示，其中762.2 nm為弱線，反演結(jié)果如圖6所示。圖6是762.2 nm無噪聲情況和加入噪聲時的反演結(jié)果比較，藍(lán)線為無噪聲情況下的結(jié)果，紅線為加入噪聲的結(jié)果。無噪聲情況下，弱線762.2 nm與強(qiáng)線764.28和764.17 nm的反演結(jié)果相對一致且都比較好； 加入噪聲后，由于本身輻射強(qiáng)度弱，所以信噪比低，反演效果變得比較差。原因是弱線容易受噪聲的干擾，抗干擾能力低。因此，譜線選擇時，線強(qiáng)也是另一個重要的判斷依據(jù)，強(qiáng)線更有利于提高反演精度。

表1 模型中用到的譜線參數(shù)Table 1 Constants of rotational lines of model

圖6 譜線762.2 nm在無噪聲情況和加入噪聲時60～110 km范圍內(nèi)的溫度反演誤差

為了進(jìn)一步確定線強(qiáng)對反演結(jié)果的影響，我們選擇通過增加線強(qiáng)的方式進(jìn)行反演精度分析，圖7為分別將線強(qiáng)增加1～5個量級后得到的反演結(jié)果。由圖可得，80 km以上，隨著線強(qiáng)量級的增大，反演精度增加，反演結(jié)果更接近初始輸入值； 線強(qiáng)增加三個量級時，反演精度可以達(dá)到5 K。由此得出： 弱線762.2 nm在信噪比足夠大，也就是線強(qiáng)達(dá)到10-26的量級時，也可用來反演80 km以上的溫度并獲得精確的反演結(jié)果。

圖7 增加譜線762.2 nm線強(qiáng)量級對反演結(jié)果的影響

4 結(jié) 論

在臨近空間大氣溫度反演中，不同的譜線反演精度不同，基于貝葉斯反演理論選取了不同的譜線進(jìn)行一系列反演試驗(yàn)及結(jié)果分析，并基于權(quán)函數(shù)和譜線線強(qiáng)與反演精度差異的關(guān)系進(jìn)行了分析。

具體結(jié)果如下：

(1)當(dāng)溫度對自吸收的影響所占比重大于對線強(qiáng)的影響時，權(quán)函數(shù)會發(fā)生正負(fù)翻轉(zhuǎn)，反演精度變差，原因是自吸收降低了輻射強(qiáng)度對溫度的靈敏度，因此，權(quán)函數(shù)變化規(guī)律可以作為譜線選擇的判斷依據(jù)。

(2)在有噪聲的情況下，強(qiáng)線比弱線的抗干擾能力更強(qiáng)，反演精度較高，更適合用于溫度的反演，所以線強(qiáng)也是譜線選擇的另一個重要的依據(jù)。根據(jù)分析結(jié)果，當(dāng)譜線線強(qiáng)達(dá)到10-26的量級時，譜線的輻射強(qiáng)度較強(qiáng)，此時譜線可用于80 km以上的溫度反演。