L曲线方法在地基红外高光谱反演温度廓线中的应用

黄 威，刘 磊，高太长，李书磊，胡 帅

解放军理工大学气象海洋学院, 江苏 南京 211101

L曲线方法在地基红外高光谱反演温度廓线中的应用

黄 威，刘 磊*，高太长，李书磊，胡 帅

解放军理工大学气象海洋学院, 江苏 南京 211101

利用地基红外高光谱辐射数据可以反演得到高时间分辨率的边界层大气温度廓线。目前的AERIoe最优化反演算法相比于传统的“剥洋葱”算法有较大的改进，且对初值的依赖程度较低。但AERIoe算法中正则化算子的选择对反演结果的稳定性和反演时间有重要影响。目前主要采用经验的方法选择正则化算子，迭代步数较多，耗费大量的计算时间。提出了利用L曲线方法代替经验法选取正则化算子的改进方案，以提高AERIoe方法的反演速度。改进后的算法通过绘制解范数和残余范数的二维曲线图，取其拐点作为最优的正则化参数，相比于传统的经验法有着更好的理论基础。采用2011年美国大气辐射测量计划中SGP站点的晴空大气红外辐射数据进行反演实验。结果表明，利用该方法得到的反演结果具有很好的稳定性、收敛性和精度。相比于经验的方法，利用L曲线方法获得的正则化算子反演温度廓线时的收敛速度更快，迭代步数较少，可以节约大量的计算时间；在反演精度方面，L曲线方法在边界层中上层的反演精度更高，1～3 km高度上温度廓线的RMSE值提高了大约0.2 K。

最优估计反演方法；L曲线；正则化算子; 高光谱

引 言

大气温度廓线在辐射传输过程、层结稳定度、降水过程以及云的形成和演变过程中具有重要的作用[1]。利用气球携带探空仪进行观测是获取大气温度廓线的重要手段，但是这种探测方式的时间分辨率较低，一般每天进行2～3次观测。地基红外高光谱仪器能够提供分辨率小于1 cm-1的下行红外辐射数据，可以用来反演得到高时间分辨率的大气温度廓线数据。相比于天基红外遥感反演，地基探测仪器受地面影响较小，在探测边界层大气温度廓线方面更有优势。

目前，利用地基红外高光谱仪反演边界层温度廓线的较成熟算法为“剥洋葱”算法(onion peeling)[2]。该算法已应用于美国威斯康辛大学研发的地基红外高光谱分辨率傅里叶变换光谱仪(atmospheric emitted radiance interferometer，AERI)的反演软件AERIprof中。但是受限于反演问题的病态特性，“剥洋葱”算法非常依赖于初值的准确性，尤其是最底层大气初值的准确程度[3]。针对“剥洋葱”方法的不足，Turner[3]提出了一种AERIoe最优化反演算法，该算法的本质是基于贝叶斯估计原理发展的一种高斯-牛顿迭代格式。为提高反演的稳定性，Turner在传统迭代格式的基础上加入了正则化参数，并基于经验的方法给出了每步迭代过程中的正则化参数值。但是这种选择方法需要的迭代步数较多，耗费了大量的计算时间；而且基于经验的方法具有较大的主观随意性，缺乏理论依据。

正则化参数在最优化反演中起着关键的作用，其选取的好坏直接影响到解的效果[4]。正则化参数的主要选取策略可以分为先验和后验两种。先验选取是指在计算正则化解之前先取定正则化参数，通常依赖于精确解的某种先验信息。因此该方法往往只有理论分析的价值，在实际应用中常常难以检验其成立的条件[5]。正则化参数的后验选取方法主要有偏差准则、广义交叉检验、L曲线等。其中L曲线法通过计算不同正则化算子的残余范数和正则化范数的二维曲线，通过寻找该曲线的拐点，能够形象直观地给出最优的正则化参数。相比于偏差准则，L曲线方法不需要观测数据的误差水平等信息；相比于广义交叉检验法，L曲线更容易计算[5]，且能更多地利用解范数的信息[6]。

重点研究了如何利用L曲线方法计算AERIoe算法中引入的正则化参数。首先介绍了使用的观测辐射数据及其预处理方法，然后对AERIoe算法、L曲线方法选取正则化参数方案和温度廓线反演流程进行了阐述，最后对L曲线方法反演结果的收敛性和计算精度进行了分析并与经验方法进行了比较。

1 使用数据及预处理方法介绍

采用了美国大气辐射测量计划(atmospheric radiation measurement，ARM)中南部大平原(south great plain，SGP)站点的数据开展研究。该站点位于北纬36°36′，西经97°29′，海拔高度320 m。下行高光谱辐射数据是由部署在该站点的AERI观测得到，其波长范围为520～3 020 cm-1。其中520～1 800 cm-1波段的红外辐射由碲铬汞(MCT)探测器获得，1 800～3 020 cm-1波段的下行红外辐射由铟化锑(InSb)探测器获得，光谱分辨率约0.5 cm-1。经过高温黑体和常温黑体定标后，AERI的探测精度可以达到1%[7]。为评估大气温度廓线的反演精度，由同站点释放的探空仪获取的温度数据作为真值，本工作所使用的均是17:30UTC的探空数据。

由于红外高光谱辐射数据含有噪声，在使用之前需要进行数据降噪等预处理。目前主要是利用主成分分析法(principal component analysis，PCA)进行降噪[8]。该方法是利用观测通道之间的相关性来降低观测数据中的噪声水平，从而提高观测数据的信噪比[9]。相比于常用的平滑方法，PCA降噪法能够保留观测数据中的有效信息，并且不影响数据的时间分辨率和光谱分辨率。由于PCA降噪方法是一个有损耗的过程，因此降噪好坏的关键是确定最优主成分的个数k，使得重构数据中的噪声水平显著降低，同时又保留观测数据中绝大部分有效信息[8]。采用经验公式法来确定最优主成分个数[8]，根据因素指标函数IND计算利用k个主成分降噪后保留在辐射数据中的噪声水平。具体公式为

(1)

其中，t表示样本个数，n是所使用的通道数，λ表示对观测数据做SVD分解后得到的特征值。所以，因素指标函数值越小，表明残留在观测数据中的噪声越少，此时计算的主成分个数即是最优的。

2 方法介绍

2.1 AERIoe反演算法

AERIoe反演算法基于最优估计方法，通过计算模拟辐射与观测辐射差值和廓线初始值之间的最小二乘关系，利用先验信息对解进行约束，并利用牛顿迭代法逐步逼近真解的最大似然估计。根据贝叶斯估计理论，在假设观测误差和背景场误差为高斯分布且二者相互独立条件下，大气状态的最优估计是使得目标函数J(X)取得极小值。

(2)

其中X是要反演的大气参数廓线，X0是大气初始廓线，Ym是观测辐射向量，Y(X)是利用辐射传输模式计算的辐射值，Se是观测误差协方差矩阵，Sa是背景协方差矩阵，γ是正则化参数。为反映辐射传输方程非线性的特性，在对辐射传输方程线性化的基础上引入牛顿非线性迭代，得到大气廓线的迭代解公式如下

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

(3)

其中Xn表示大气廓线第n次迭代的结果，Xn+1表示第n+1次迭代结果，Kn是利用第n次迭代结果Xn计算的雅克比矩阵，Y(Xn)表示利用Xn通过辐射传输模式计算得到的模拟辐射光谱。

在迭代公式中，正则化算子γ产生了平滑的效果，其作用主要是通过阻尼函数xLSQ对辐射传输方程的最小二乘解进行过滤。

(4)

其中σi，ui和vi分别表示对辐射传输方程做SVD分解得到的特征值和特征向量，b表示和观测数据相关的项，fi表示阻尼函数，其大小与正则化算子成反比，当σi=1有fi→0。显然，较小σi对b中的观测误差有放大作用，而阻尼函数fi则可以抑制这种放大效果，γ越大这种抑制效果越好；但是，正则化参数的引入会丢失一部分信息从而带来正则化误差，γ越大正则化误差也越大。

L曲线方法正是通过寻找对σi的抑制效果和正则化误差之间的平衡点来选择最优的正则化参数。该方法最早由Lawson和Hanson在1974年提出，并迅速在地球物理、遥感技术和生命科学等与反问题相关的领域获得了应用[6]。但目前尚未运用到地基高光谱反演算法中。

2.2 L曲线方法选择正则化算子

在目标函数J(X)中，正则化算子γ控制着解范数和残余范数的权重，因此正则化解中包含有解范数和残余范数的信息。上述分析可知，最优的正则化算子是使得产生的正则化误差与抑制效果取得平衡，该问题等价于寻找使得解范数和残余范数同时达到最小值时的正则化参数[6]。通过绘制出残余范数与正则化解范数随正则化参数变化间的二维曲线图，该曲线图像通常呈现出字母L的形状，垂直的部分表示正则化对解所造成的误差大小，水平部分表示观测辐射中的噪声对正则化解产生的影响，因此L曲线的拐点表示二者同时达到极小值的区域，此时图形的拐点即为正则化算子的最优值[10]。

在计算L曲线时，首先应将辐射传输方程线型化

AX=b

(5)

(6)

(7)

曲率κ取最大值时对应的γ即为最优的正则化参数。

2.3 利用L曲线法反演温度廓线的流程

利用L曲线法反演温度廓线的流程如图1所示。主要有以下6个步骤。

第1步，对晴空大气条件下的温度廓线求平均值得到初始廓线。

第2步，将晴空数据样本集中的温度廓线放入辐射传输模式中计算模拟辐射；之后，挑选出与探空仪时间匹配的观测辐射光谱，根据得到的模拟辐射光谱，利用多个样本资料计算二者的平均偏差，从而得到仪器的系统偏差光谱。

第3步： 将观测辐射减去系统偏差光谱，之后利用PCA做降噪处理。采用612～618，624～660，674～713和2 223～2 260 cm-1的波段来反演温度廓线。其中前三个波段由MCT探测器获取，第四个波段由InSb探测器观测得到，两个探测器的噪声水平并不一致，因此需要分开做降噪处理[8]。

图1 利用L曲线方法的反演流程图

第4步： 利用辐射传输模式计算初始廓线对应的模拟辐射光谱和雅克比矩阵Kn。其中，雅克比矩阵的计算是基于扰动法，具体过程见文献[12]： 对第i层的温度值增加和减少一个小量(0.5 K)，此时第i层的雅克比矩阵为

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

(8)

其中TB表示传输模式计算的亮温值。

第5步： 根据L曲线法确定正则化参数，具体过程参见2.2。

第6步： 将雅克比矩阵、初始廓线、降噪后的观测辐射以及正则化参数放入式(3)中进行迭代。其中观测误差协方差矩阵Se是在假定各通道之间的噪声不相关的条件下取为对角阵[13]；背景协方差矩阵Sa是通过该站点的2010年的晴空温度廓线数据，计算其协方差矩阵得到[3]。若迭代结果满足收敛条件，则输出廓线；否则，将调整后的廓线再次放入模式中，重复第4步及其以后的步骤。

3 反演结果分析

为检验L曲线方法在温度廓线反演中的应用表现，选用ARM计划中SGP站点的2011年的观测辐射资料进行反演实验。根据该站点总天空成像仪(total sky imager，TSI)获取的全天空云图，从2011年的观测辐射数据中挑选出9组具有代表性的数据进行反演。这9组数据涵盖了春、夏、秋、冬四个季节。由于观测数据中85%的温度廓线信息集中在2 km以下[3]，使得AERI反演的高度一般不超过3 km[14]。在反演实验中，将0.32～3 km的大气分为19层，正则化算子用L曲线的方法进行确定。

图2为2011年6月6日17:30UTC的温度廓线反演结果。其中黑色实线表示初始廓线，通过计算2010年晴空条件下温度廓线的平均值得到；实线表示探空仪的观测值，此处作为大气“真实”的状态对反演结果的好坏进行评估；带点缀的线条表示每步迭代的反演结果。可以看出，随着迭代次数的增加反演结果逐渐逼近探空廓线。从第2步开始，迭代解逐渐收敛，第3至7步的迭代解几乎完全重合。为分析L曲线方法在选取正则化算子方面的优劣，进一步对其在收敛性和反演精度等方面与经验方法的反演结果进行了对比，并对其反演的稳定性进行了分析。

图2 利用L曲线方法的温度廓线反演结果

3.1 收敛性和稳定性分析

AERIoe算法中采用了经验的方法选择正则化参数。其正则化算子γ取为固定的值[1 000 300 100 30 10 3 1]，其收敛性的判断是迭代至γ=1时判断解是否满足如下收敛条件

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

(9)

表1 不同迭代步骤的精度

其中，平均绝对偏差biasabsolute表示每一步反演的值与探空仪观测值在整个反演高度内的绝对偏差的平均值，用以下公式计算。Xsonde表示探空廓线，Xretrieval表示每一步迭代的反演结果。

(10)

从表1中可以看出，第3步和第4步解的平均偏差相差较小，均满足收敛条件(9)。因此，利用L曲线方法选取正则化参数，使得温度廓线的反演具有更快的收敛性。由于雅克比矩阵的计算时间较长，在每一次迭代过程中需要根据当时的迭代解重新计算雅克比矩阵，迭代次数的减少，大大缩短了温度廓线反演所需的计算时间。

由于使用的初始廓线是利用过去一年的晴空大气探空资料计算的平均值，所挑选的9个晴空样本有很长的时间跨度，这就不可避免的出现初始廓线与大气真值相差较大的情况。表2是不同初始绝对偏差情况下温度廓线的反演结果。其中，初始绝对偏差与迭代结果均是根据式(10)进行计算。

表2 不同初始绝对偏差条件下的收敛结果

从表中可以看出，在初始偏差15.28 K的情况下，解依然收敛，并且取得了0.64 K的反演精度；对于6月6日的反演结果，从图2中可以看出，其平均偏差较大是由边界层上层反演的廓线精度较低造成，由于观测数据中包含的边界层上层的信息较少，在初始偏差较大的情况下，其反演的精度就会变差。而此时第3至7步的迭代解几乎完全重合，表明该方法收敛性对初始廓线的依赖较低，具有较好的稳定性。

3.2 误差分析

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

(11)

(12)

图3是利用L曲线方法和经验方法的反演结果与探空廓线的平均偏差和均方根误差随高度变化。其中实线表示利用L曲线方法得到的结果，虚线则表示经验法的结果；蓝色线条表示平均偏差，红色为均方根误差。从平均偏差图像可以看出，除近地面位置外，无论是利用经验的方法还是L曲线方法计算γ，AERIoe方法的反演结果相比于探空仪的观测值均偏大。其中在1～2 km的高度上，L曲线方法偏差更小一些；从RMSE的垂直分布上来看，在1 km以上L曲线方法的反演结果要明显优于经验的方法，其RMSE值提高了约0.2 K。以上分析表明，利用L曲线方法计算正则化参数的反演精度要优于经验法的反演结果。

图3 L曲线方法和经验法得到的温度廓线反演的RMSE和平均偏差的垂直分布

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

4 结 论

在地基红外高光谱温度廓线遥感反演方面，结合正则化算子的最优化反演方法(AERIoe)是一种较新的方法。相比于“剥洋葱”算法，该方法具有对初始廓线依赖性较小，稳定性更高等诸多优点。利用L曲线方法代替经验方法选取正则化算子，采用2011年9组晴空观测数据对该方法的收敛性、稳定性和反演精度进行了分析。结果表明，利用L曲线方法选取正则化算子后，AERIoe的反演结果具有很好的稳定性；相比于经验的方法，改进后的方法在边界层上层的反演精度更高，且具有更快的收敛速度，能够节省大量的计算时间。

[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(地球物理学进展)，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(贵州大学学报·自然科学版)，2011，28(4): 29.

[12] GUAN Li(官 莉). Application of Satellite Borne Infrared Hyperctral Data(星载红外高光谱资料的应用). 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

国家自然科学基金项目(41575024)资助

黄 威，1992年生，解放军理工大学气象海洋学院硕士研究生 e-mail: huangwei_edu@sina.com *通讯联系人 e-mail: liuleidll@gmail.com

TP722.5

A

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