L曲線方法在地基紅外高光譜反演溫度廓線中的應(yīng)用

黃 威，劉 磊，高太長，李書磊，胡 帥

解放軍理工大學(xué)氣象海洋學(xué)院, 江蘇 南京 211101

L曲線方法在地基紅外高光譜反演溫度廓線中的應(yīng)用

黃 威，劉 磊*，高太長，李書磊，胡 帥

解放軍理工大學(xué)氣象海洋學(xué)院, 江蘇 南京 211101

利用地基紅外高光譜輻射數(shù)據(jù)可以反演得到高時間分辨率的邊界層大氣溫度廓線。目前的AERIoe最優(yōu)化反演算法相比于傳統(tǒng)的“剝洋蔥”算法有較大的改進，且對初值的依賴程度較低。但AERIoe算法中正則化算子的選擇對反演結(jié)果的穩(wěn)定性和反演時間有重要影響。目前主要采用經(jīng)驗的方法選擇正則化算子，迭代步數(shù)較多，耗費大量的計算時間。提出了利用L曲線方法代替經(jīng)驗法選取正則化算子的改進方案，以提高AERIoe方法的反演速度。改進后的算法通過繪制解范數(shù)和殘余范數(shù)的二維曲線圖，取其拐點作為最優(yōu)的正則化參數(shù)，相比于傳統(tǒng)的經(jīng)驗法有著更好的理論基礎(chǔ)。采用2011年美國大氣輻射測量計劃中SGP站點的晴空大氣紅外輻射數(shù)據(jù)進行反演實驗。結(jié)果表明，利用該方法得到的反演結(jié)果具有很好的穩(wěn)定性、收斂性和精度。相比于經(jīng)驗的方法，利用L曲線方法獲得的正則化算子反演溫度廓線時的收斂速度更快，迭代步數(shù)較少，可以節(jié)約大量的計算時間；在反演精度方面，L曲線方法在邊界層中上層的反演精度更高，1～3 km高度上溫度廓線的RMSE值提高了大約0.2 K。

最優(yōu)估計反演方法；L曲線；正則化算子; 高光譜

引 言

大氣溫度廓線在輻射傳輸過程、層結(jié)穩(wěn)定度、降水過程以及云的形成和演變過程中具有重要的作用[1]。利用氣球攜帶探空儀進行觀測是獲取大氣溫度廓線的重要手段，但是這種探測方式的時間分辨率較低，一般每天進行2～3次觀測。地基紅外高光譜儀器能夠提供分辨率小于1 cm-1的下行紅外輻射數(shù)據(jù)，可以用來反演得到高時間分辨率的大氣溫度廓線數(shù)據(jù)。相比于天基紅外遙感反演，地基探測儀器受地面影響較小，在探測邊界層大氣溫度廓線方面更有優(yōu)勢。

目前，利用地基紅外高光譜儀反演邊界層溫度廓線的較成熟算法為“剝洋蔥”算法(onion peeling)[2]。該算法已應(yīng)用于美國威斯康辛大學(xué)研發(fā)的地基紅外高光譜分辨率傅里葉變換光譜儀(atmospheric emitted radiance interferometer，AERI)的反演軟件AERIprof中。但是受限于反演問題的病態(tài)特性，“剝洋蔥”算法非常依賴于初值的準(zhǔn)確性，尤其是最底層大氣初值的準(zhǔn)確程度[3]。針對“剝洋蔥”方法的不足，Turner[3]提出了一種AERIoe最優(yōu)化反演算法，該算法的本質(zhì)是基于貝葉斯估計原理發(fā)展的一種高斯-牛頓迭代格式。為提高反演的穩(wěn)定性，Turner在傳統(tǒng)迭代格式的基礎(chǔ)上加入了正則化參數(shù)，并基于經(jīng)驗的方法給出了每步迭代過程中的正則化參數(shù)值。但是這種選擇方法需要的迭代步數(shù)較多，耗費了大量的計算時間；而且基于經(jīng)驗的方法具有較大的主觀隨意性，缺乏理論依據(jù)。

正則化參數(shù)在最優(yōu)化反演中起著關(guān)鍵的作用，其選取的好壞直接影響到解的效果[4]。正則化參數(shù)的主要選取策略可以分為先驗和后驗兩種。先驗選取是指在計算正則化解之前先取定正則化參數(shù)，通常依賴于精確解的某種先驗信息。因此該方法往往只有理論分析的價值，在實際應(yīng)用中常常難以檢驗其成立的條件[5]。正則化參數(shù)的后驗選取方法主要有偏差準(zhǔn)則、廣義交叉檢驗、L曲線等。其中L曲線法通過計算不同正則化算子的殘余范數(shù)和正則化范數(shù)的二維曲線，通過尋找該曲線的拐點，能夠形象直觀地給出最優(yōu)的正則化參數(shù)。相比于偏差準(zhǔn)則，L曲線方法不需要觀測數(shù)據(jù)的誤差水平等信息；相比于廣義交叉檢驗法，L曲線更容易計算[5]，且能更多地利用解范數(shù)的信息[6]。

重點研究了如何利用L曲線方法計算AERIoe算法中引入的正則化參數(shù)。首先介紹了使用的觀測輻射數(shù)據(jù)及其預(yù)處理方法，然后對AERIoe算法、L曲線方法選取正則化參數(shù)方案和溫度廓線反演流程進行了闡述，最后對L曲線方法反演結(jié)果的收斂性和計算精度進行了分析并與經(jīng)驗方法進行了比較。

1 使用數(shù)據(jù)及預(yù)處理方法介紹

采用了美國大氣輻射測量計劃(atmospheric radiation measurement，ARM)中南部大平原(south great plain，SGP)站點的數(shù)據(jù)開展研究。該站點位于北緯36°36′，西經(jīng)97°29′，海拔高度320 m。下行高光譜輻射數(shù)據(jù)是由部署在該站點的AERI觀測得到，其波長范圍為520～3 020 cm-1。其中520～1 800 cm-1波段的紅外輻射由碲鉻汞(MCT)探測器獲得，1 800～3 020 cm-1波段的下行紅外輻射由銦化銻(InSb)探測器獲得，光譜分辨率約0.5 cm-1。經(jīng)過高溫黑體和常溫黑體定標(biāo)后，AERI的探測精度可以達(dá)到1%[7]。為評估大氣溫度廓線的反演精度，由同站點釋放的探空儀獲取的溫度數(shù)據(jù)作為真值，本工作所使用的均是17:30UTC的探空數(shù)據(jù)。

由于紅外高光譜輻射數(shù)據(jù)含有噪聲，在使用之前需要進行數(shù)據(jù)降噪等預(yù)處理。目前主要是利用主成分分析法(principal component analysis，PCA)進行降噪[8]。該方法是利用觀測通道之間的相關(guān)性來降低觀測數(shù)據(jù)中的噪聲水平，從而提高觀測數(shù)據(jù)的信噪比[9]。相比于常用的平滑方法，PCA降噪法能夠保留觀測數(shù)據(jù)中的有效信息，并且不影響數(shù)據(jù)的時間分辨率和光譜分辨率。由于PCA降噪方法是一個有損耗的過程，因此降噪好壞的關(guān)鍵是確定最優(yōu)主成分的個數(shù)k，使得重構(gòu)數(shù)據(jù)中的噪聲水平顯著降低，同時又保留觀測數(shù)據(jù)中絕大部分有效信息[8]。采用經(jīng)驗公式法來確定最優(yōu)主成分個數(shù)[8]，根據(jù)因素指標(biāo)函數(shù)IND計算利用k個主成分降噪后保留在輻射數(shù)據(jù)中的噪聲水平。具體公式為

(1)

其中，t表示樣本個數(shù)，n是所使用的通道數(shù)，λ表示對觀測數(shù)據(jù)做SVD分解后得到的特征值。所以，因素指標(biāo)函數(shù)值越小，表明殘留在觀測數(shù)據(jù)中的噪聲越少，此時計算的主成分個數(shù)即是最優(yōu)的。

2 方法介紹

2.1 AERIoe反演算法

AERIoe反演算法基于最優(yōu)估計方法，通過計算模擬輻射與觀測輻射差值和廓線初始值之間的最小二乘關(guān)系，利用先驗信息對解進行約束，并利用牛頓迭代法逐步逼近真解的最大似然估計。根據(jù)貝葉斯估計理論，在假設(shè)觀測誤差和背景場誤差為高斯分布且二者相互獨立條件下，大氣狀態(tài)的最優(yōu)估計是使得目標(biāo)函數(shù)J(X)取得極小值。

(2)

其中X是要反演的大氣參數(shù)廓線，X0是大氣初始廓線，Ym是觀測輻射向量，Y(X)是利用輻射傳輸模式計算的輻射值，Se是觀測誤差協(xié)方差矩陣，Sa是背景協(xié)方差矩陣，γ是正則化參數(shù)。為反映輻射傳輸方程非線性的特性，在對輻射傳輸方程線性化的基礎(chǔ)上引入牛頓非線性迭代，得到大氣廓線的迭代解公式如下

(Ym-Y(Xn)+Kn(Xn-X0))

(3)

其中Xn表示大氣廓線第n次迭代的結(jié)果，Xn+1表示第n+1次迭代結(jié)果，Kn是利用第n次迭代結(jié)果Xn計算的雅克比矩陣，Y(Xn)表示利用Xn通過輻射傳輸模式計算得到的模擬輻射光譜。

在迭代公式中，正則化算子γ產(chǎn)生了平滑的效果，其作用主要是通過阻尼函數(shù)xLSQ對輻射傳輸方程的最小二乘解進行過濾。

(4)

其中σi，ui和vi分別表示對輻射傳輸方程做SVD分解得到的特征值和特征向量，b表示和觀測數(shù)據(jù)相關(guān)的項，fi表示阻尼函數(shù)，其大小與正則化算子成反比，當(dāng)σi=1有fi→0。顯然，較小σi對b中的觀測誤差有放大作用，而阻尼函數(shù)fi則可以抑制這種放大效果，γ越大這種抑制效果越好；但是，正則化參數(shù)的引入會丟失一部分信息從而帶來正則化誤差，γ越大正則化誤差也越大。

L曲線方法正是通過尋找對σi的抑制效果和正則化誤差之間的平衡點來選擇最優(yōu)的正則化參數(shù)。該方法最早由Lawson和Hanson在1974年提出，并迅速在地球物理、遙感技術(shù)和生命科學(xué)等與反問題相關(guān)的領(lǐng)域獲得了應(yīng)用[6]。但目前尚未運用到地基高光譜反演算法中。

2.2 L曲線方法選擇正則化算子

在目標(biāo)函數(shù)J(X)中，正則化算子γ控制著解范數(shù)和殘余范數(shù)的權(quán)重，因此正則化解中包含有解范數(shù)和殘余范數(shù)的信息。上述分析可知，最優(yōu)的正則化算子是使得產(chǎn)生的正則化誤差與抑制效果取得平衡，該問題等價于尋找使得解范數(shù)和殘余范數(shù)同時達(dá)到最小值時的正則化參數(shù)[6]。通過繪制出殘余范數(shù)與正則化解范數(shù)隨正則化參數(shù)變化間的二維曲線圖，該曲線圖像通常呈現(xiàn)出字母L的形狀，垂直的部分表示正則化對解所造成的誤差大小，水平部分表示觀測輻射中的噪聲對正則化解產(chǎn)生的影響，因此L曲線的拐點表示二者同時達(dá)到極小值的區(qū)域，此時圖形的拐點即為正則化算子的最優(yōu)值[10]。

在計算L曲線時，首先應(yīng)將輻射傳輸方程線型化

AX=b

(5)

(6)

(7)

曲率κ取最大值時對應(yīng)的γ即為最優(yōu)的正則化參數(shù)。

2.3 利用L曲線法反演溫度廓線的流程

利用L曲線法反演溫度廓線的流程如圖1所示。主要有以下6個步驟。

第1步，對晴空大氣條件下的溫度廓線求平均值得到初始廓線。

第2步，將晴空數(shù)據(jù)樣本集中的溫度廓線放入輻射傳輸模式中計算模擬輻射；之后，挑選出與探空儀時間匹配的觀測輻射光譜，根據(jù)得到的模擬輻射光譜，利用多個樣本資料計算二者的平均偏差，從而得到儀器的系統(tǒng)偏差光譜。

第3步： 將觀測輻射減去系統(tǒng)偏差光譜，之后利用PCA做降噪處理。采用612～618，624～660，674～713和2 223～2 260 cm-1的波段來反演溫度廓線。其中前三個波段由MCT探測器獲取，第四個波段由InSb探測器觀測得到，兩個探測器的噪聲水平并不一致，因此需要分開做降噪處理[8]。

圖1 利用L曲線方法的反演流程圖

第4步： 利用輻射傳輸模式計算初始廓線對應(yīng)的模擬輻射光譜和雅克比矩陣Kn。其中，雅克比矩陣的計算是基于擾動法，具體過程見文獻(xiàn)[12]： 對第i層的溫度值增加和減少一個小量(0.5 K)，此時第i層的雅克比矩陣為

KT, i=TB(Ti+0.5)-TB(Ti-0.5)

(8)

其中TB表示傳輸模式計算的亮溫值。

第5步： 根據(jù)L曲線法確定正則化參數(shù)，具體過程參見2.2。

第6步： 將雅克比矩陣、初始廓線、降噪后的觀測輻射以及正則化參數(shù)放入式(3)中進行迭代。其中觀測誤差協(xié)方差矩陣Se是在假定各通道之間的噪聲不相關(guān)的條件下取為對角陣[13]；背景協(xié)方差矩陣Sa是通過該站點的2010年的晴空溫度廓線數(shù)據(jù)，計算其協(xié)方差矩陣得到[3]。若迭代結(jié)果滿足收斂條件，則輸出廓線；否則，將調(diào)整后的廓線再次放入模式中，重復(fù)第4步及其以后的步驟。

3 反演結(jié)果分析

為檢驗L曲線方法在溫度廓線反演中的應(yīng)用表現(xiàn)，選用ARM計劃中SGP站點的2011年的觀測輻射資料進行反演實驗。根據(jù)該站點總天空成像儀(total sky imager，TSI)獲取的全天空云圖，從2011年的觀測輻射數(shù)據(jù)中挑選出9組具有代表性的數(shù)據(jù)進行反演。這9組數(shù)據(jù)涵蓋了春、夏、秋、冬四個季節(jié)。由于觀測數(shù)據(jù)中85%的溫度廓線信息集中在2 km以下[3]，使得AERI反演的高度一般不超過3 km[14]。在反演實驗中，將0.32～3 km的大氣分為19層，正則化算子用L曲線的方法進行確定。

圖2為2011年6月6日17:30UTC的溫度廓線反演結(jié)果。其中黑色實線表示初始廓線，通過計算2010年晴空條件下溫度廓線的平均值得到；實線表示探空儀的觀測值，此處作為大氣“真實”的狀態(tài)對反演結(jié)果的好壞進行評估；帶點綴的線條表示每步迭代的反演結(jié)果?？梢钥闯觯S著迭代次數(shù)的增加反演結(jié)果逐漸逼近探空廓線。從第2步開始，迭代解逐漸收斂，第3至7步的迭代解幾乎完全重合。為分析L曲線方法在選取正則化算子方面的優(yōu)劣，進一步對其在收斂性和反演精度等方面與經(jīng)驗方法的反演結(jié)果進行了對比，并對其反演的穩(wěn)定性進行了分析。

圖2 利用L曲線方法的溫度廓線反演結(jié)果

3.1 收斂性和穩(wěn)定性分析

AERIoe算法中采用了經(jīng)驗的方法選擇正則化參數(shù)。其正則化算子γ取為固定的值[1 000 300 100 30 10 3 1]，其收斂性的判斷是迭代至γ=1時判斷解是否滿足如下收斂條件

(Xn-Xn+1)TS-1(Xn-Xn+1)

(9)

表1 不同迭代步驟的精度

其中，平均絕對偏差biasabsolute表示每一步反演的值與探空儀觀測值在整個反演高度內(nèi)的絕對偏差的平均值，用以下公式計算。Xsonde表示探空廓線，Xretrieval表示每一步迭代的反演結(jié)果。

(10)

從表1中可以看出，第3步和第4步解的平均偏差相差較小，均滿足收斂條件(9)。因此，利用L曲線方法選取正則化參數(shù)，使得溫度廓線的反演具有更快的收斂性。由于雅克比矩陣的計算時間較長，在每一次迭代過程中需要根據(jù)當(dāng)時的迭代解重新計算雅克比矩陣，迭代次數(shù)的減少，大大縮短了溫度廓線反演所需的計算時間。

由于使用的初始廓線是利用過去一年的晴空大氣探空資料計算的平均值，所挑選的9個晴空樣本有很長的時間跨度，這就不可避免的出現(xiàn)初始廓線與大氣真值相差較大的情況。表2是不同初始絕對偏差情況下溫度廓線的反演結(jié)果。其中，初始絕對偏差與迭代結(jié)果均是根據(jù)式(10)進行計算。

表2 不同初始絕對偏差條件下的收斂結(jié)果

從表中可以看出，在初始偏差15.28 K的情況下，解依然收斂，并且取得了0.64 K的反演精度；對于6月6日的反演結(jié)果，從圖2中可以看出，其平均偏差較大是由邊界層上層反演的廓線精度較低造成，由于觀測數(shù)據(jù)中包含的邊界層上層的信息較少，在初始偏差較大的情況下，其反演的精度就會變差。而此時第3至7步的迭代解幾乎完全重合，表明該方法收斂性對初始廓線的依賴較低，具有較好的穩(wěn)定性。

3.2 誤差分析

為表征在各高度上的反演精度，引入平均偏差bias和均方根誤差(RMSE)，具體計算公式如下所示

(11)

(12)

圖3是利用L曲線方法和經(jīng)驗方法的反演結(jié)果與探空廓線的平均偏差和均方根誤差隨高度變化。其中實線表示利用L曲線方法得到的結(jié)果，虛線則表示經(jīng)驗法的結(jié)果；藍(lán)色線條表示平均偏差，紅色為均方根誤差。從平均偏差圖像可以看出，除近地面位置外，無論是利用經(jīng)驗的方法還是L曲線方法計算γ，AERIoe方法的反演結(jié)果相比于探空儀的觀測值均偏大。其中在1～2 km的高度上，L曲線方法偏差更小一些；從RMSE的垂直分布上來看，在1 km以上L曲線方法的反演結(jié)果要明顯優(yōu)于經(jīng)驗的方法，其RMSE值提高了約0.2 K。以上分析表明，利用L曲線方法計算正則化參數(shù)的反演精度要優(yōu)于經(jīng)驗法的反演結(jié)果。

圖3 L曲線方法和經(jīng)驗法得到的溫度廓線反演的RMSE和平均偏差的垂直分布

Fig.3 The vertical distribution of the RMSE and bias of the retrieval temperature by L-curve method and empirical method

4 結(jié) 論

在地基紅外高光譜溫度廓線遙感反演方面，結(jié)合正則化算子的最優(yōu)化反演方法(AERIoe)是一種較新的方法。相比于“剝洋蔥”算法，該方法具有對初始廓線依賴性較小，穩(wěn)定性更高等諸多優(yōu)點。利用L曲線方法代替經(jīng)驗方法選取正則化算子，采用2011年9組晴空觀測數(shù)據(jù)對該方法的收斂性、穩(wěn)定性和反演精度進行了分析。結(jié)果表明，利用L曲線方法選取正則化算子后，AERIoe的反演結(jié)果具有很好的穩(wěn)定性；相比于經(jīng)驗的方法，改進后的方法在邊界層上層的反演精度更高，且具有更快的收斂速度，能夠節(jié)省大量的計算時間。

[1] Volker Wulfmeyer R M H D. Rev. Geophys.，2015, 53(3)： 819.

[2] Smith W L, Feltz W F, Knuteson R O, et al. Journal of Atmospheric and Oceanic Technology，1999，16(2-3): 323.

[3] Turner D D, L?hnert U. Journal of Applied Meteorology and Climatology，2014，53(3): 752.

[4] MA Tao, CHEN Long-wei, WU Mei-ping, et al(馬 濤，陳龍偉，吳美平，等). Progress in Geophysics(地球物理學(xué)進展)，2013，28(5)： 2485.

[5] Hansen P C, O’Leary D P, SIAM Journal of Scientific Computing, 1993, 14(6): 1487.

[6] Hansen P C. SIAM Review，1992, 34(4): 561.

[7] Knuteson R O, Revercomb H E, Best F A, et al. Journal of Atmospheric and Oceanic Technology，2004，21(12): 1777.

[8] Turner D D, Knuteson R O, Revercomb H E, et al. Journal of Atmospheric and Oceanic Technology，2006，23(9): 1223.

[9] Antonelli P, Revercomb H E, Sromovsky L A, et al. Journal of Geophysical Research, 2004(D23).

[10] Hansen P C, Jensen T K, Rodriguez G. Journal of Computational and Applied Mathematics，2007, 198(2): 483.

[11] FAN Qian, FANG Xu-hua, FAN Juan(范 千，方緒華，范 娟). Journal of Guizhou University·Natural Sciences(貴州大學(xué)學(xué)報·自然科學(xué)版)，2011，28(4): 29.

[12] GUAN Li(官 莉). Application of Satellite Borne Infrared Hyperctral Data(星載紅外高光譜資料的應(yīng)用). Beijing: China Meterological Press(北京: 氣象出版社), 2007.

[13] Merrelli A, Turner D D. EN，2012，(4): 510.

[14] Feltz W F, Smith W L, Knuteson R O, et al. Journal of Applied Meteorology，1998，37(9): 857.

(Received Feb. 22, 2016; accepted May 16, 2016)

*Corresponding author

The Application of the L-Curve Method in the Retrieval of Temperature Profiles Using Ground-Based Hyper-Spectral Infrared Radiance

HUANG Wei, LIU Lei*, GAO Tai-chang, LI Shu-lei, HU Shuai

College of Meteorology and Oceanography, the PLA University of Science and Technology, Nanjing 211101, China

The thermodynamic profiles of Planetary Boundary Layer could be retrieved by using ground-based hyper-spectral infrared radiance. The AERIoe algorithm has a better performance at the dependency of initial profiles than the “onion peeling” method which was originally applied in the Atmospheric Emitted Radiance Interferometer. The regularization parameter is the key to the AERIoe algorithm, and the strategy for choosing the regularization parameter in the retrieval algorithm is based on the empirical method, which requires too much time for computation while the empirical method needs many iteration steps. A L-curve criterion is proposed to calculate the regularization parameter in AERIoe algorithm. The L-curve criterion is based on a log-log plot of corresponding values of the residual and solution norms, and the optimal regularization parameter corresponds to a point on the curve near the “corner” of the L-shaped region. Therefore, the L-curve criterion has better theoretical basis than the traditional empirical method. The result of retrieval experiment using the observed data collected at the SGP site of the year 2011 shows that, the L-curve method has a good performance in terms of stability, convergence and accuracy of the retrieval. Compared with empirical method, L-curve algorithm converges more quickly which saves much computation time when retrieving the temperature profiles. When considering the retrieval accuracy, the L-curve method has a better behavior at the middle and top of the boundary layer, with an improvement of 0.2 K of RMSE at the altitude of 1～3 km than the empirical method. Therefore, the L-curve algorithm has a better performance compared with the empirical method when choosing the regularization parameter in the retrieval of temperature profiles using the ground-based hyper-spectral infrared radiance.

Optimal inverse methods; L-curve; Regularization parameter; Hyper-spectral

2016-02-22，

2016-05-16

國家自然科學(xué)基金項目(41575024)資助

黃 威，1992年生，解放軍理工大學(xué)氣象海洋學(xué)院碩士研究生 e-mail: huangwei_edu@sina.com *通訊聯(lián)系人 e-mail: liuleidll@gmail.com

TP722.5

A

10.3964/j.issn.1000-0593(2016)11-3620-05