CN103792538B - 一种基于地基高光谱微波辐射计的大气廓线反演方法 - Google Patents

Abstract

本发明提供了一种基于地基高光谱微波辐射计的大气廓线反演方法，属于微波遥感技术领域。本方法包括以下步骤：设置地基微波辐射计的系统参数；获取背景参数以及背景场误差；建立辐射传输模型，确定大气辐射传输亮温；利用变分方法的迭代过程反演大气廓线。利用本发明可以通过地基高光谱微波辐射计得到精确反演的大气廓线，特别适用于多通道辐射计和先验数据不足的探测情况，能够完整地利用高光谱辐射计提供的大量频段信息，具有较强的适用性和较小的探测误差。

Description

一种基于地基高光谱微波辐射计的大气廓线反演方法

技术领域

本发明属于微波遥感技术领域，具体涉及一种基于地基高光谱微波辐射计的大气廓线反演方法。

背景技术

无源微波遥感传感器的最大通道数量大部分为数个或十数个，很少有几十个通道的微波辐射计，因为通道数量受限于现有技术。但是这些通道数量通常不能满足大气反演问题，尤其是在云和雨存在的情况下，由于待反演变量的数量增多使得病态反演问题严重，不能有效反演出待测参数，从而影响反演精度。

地基微波辐射计因其在低对流层具有更高的探测精度，并且能够对局部区域实现连续不间断地观测等优势，已成为大气廓线探测的重要工具。目前国际上比较成熟的推向市场的地基微波辐射计主要有美国Radiometrics公司的MP-3000系列微波辐射计以及前德国RPG公司RPG-xxxPRO系列微波辐射计。其无源微波遥感探测传感器探测通道数目大部分为5～30通道，但这些受限的通道数目不足以解决探测大气参数时的病态反演问题，特别是有云雨的情况下，由于待反演变量增多导致反演结果精度有限。随着微波器件加工技术与数字信号处理芯片技术的提高，微波探测传感器可以在一个频段内精细化分，从而使探测通道数目达到成百上千个。通道数目多，信息量就大，并且能够为探测大气参数提供更高的精度和垂直、水平分辨率，从而提高短时间小范围内的天气预测能力，改善数值天气的预报，特别是极端天气的预测能力。

但是，目前的大气廓线反演方法针对的微波辐射计通道数量较少。而随着通道数目的增多，相邻通道之间的相关性也随之增加，这就为大气廓线的反演带来麻烦。因为通常的大气参数反演方法大都是包含矩阵求逆的过程，而当通道数量多且相邻通道之间存在一定的相关性时，矩阵的求逆就变得十分困难，从而不能正确地反演相应大气参数。神经网络方法，对于通道数目大于20的微波辐射计，虽然能利用主元素分析法来降低维度，只选择一部分的探测通道进行计算，但是同时也减少了探测数据的信息量，失去了高光谱大量通道探测的意义，在一定程度上也减弱了高光谱通道数目多信息量大的优势，不能完全发挥高光谱辐射计的作用。

发明内容

本发明为充分利用高光谱高信息量的特点，不损失有效的探测信息，同时获得更精确的大气廓线反演结果，提出了一种基于地基高光谱微波辐射计的大气廓线反演方法。

本发明提供的基于地基高光谱微波辐射计的大气廓线反演方法，包括如下步骤：

步骤1：设置地基微波辐射计的系统参数，包括探测频段，通道数量以及各通道频率和辐射计的探测误差。

步骤2：获取背景参数以及背景场误差。背景参数表示为大气状态的平均值x_b，通过统计被反演地区探测日前具有共同特点的历史大气廓线探测结果得到，是大气廓线反演的先验数据。背景场误差反映了被反演地区大气廓线的准确性，表示为通过先验数据得到的背景场协方差矩阵B。

步骤3：建立辐射传输模型，确定大气辐射传输亮温。

所采用的辐射传输模型如下：

$T_B=\mathrm{\μ}{\∫}_0^d\mathrm{dz}{\mathrm{\κ}}_a\left(z\right)T\left(z\right)\exp (-\mathrm{\μ}{\∫}_0^z{\mathrm{dz}}^{\′}{\mathrm{\κ}}_a\left(z^{\′}\right))$

其中，T_B为大气辐射传输亮温，参数μ＝secθ，θ为探测角度，d表示地面到大气顶层的高度，κ_a(z)为高度z处的吸收率，T(z)为高度z处的大气温度。

所述辐射传输模型，以大气温湿度廓线为输入，通过MPM(Maximizer of theposteriori marginal,最大后验边缘概率)模型根据探测频率计算吸收率，再结合不同高度的温度廓线，计算微波辐射计接收的辐射传输亮温。

步骤4：利用变分方法的迭代过程反演大气廓线，分为以下几个步骤：

步骤4.1：建立辐射计测量的大气辐射传输亮温与待反演的大气状态的关系；

记x＝(x₁,x₂,…x_m)为描述大气状态的向量，为待反演的未知量，m为被反演地区研究范围内的高度层，x_m为m高度层的大气状态值；记y＝(y₁,y₂,…y_n)为辐射计测量的大气辐射传输亮温的向量，n为辐射计的通道数目，y_n为第n个通道的大气辐射传输亮温。

用前向模型F(x)表示y与x的关系：y＝F(x)+ε；ε为测量误差，是一个随机变量。

步骤4.2：确定计算待反演未知量的目标函数；

所述的目标函数J(x)如下：

J(x)＝[x-x_b]^TB^-1[x-x_b]+[y-F(x)]^TR^-1[y-F(x)]

通过下面方程寻找使得目标函数J(x)取得最小值时的x：

H(x)^TR^-1[y-F(x)]+B^-1[x-x_b]＝0

其中，雅克比矩阵包括辐射计通道的大气辐射传输亮温对每个大气状态的偏微分；R为测量误差ε的协方差矩阵；上角标T和-1分别表示矩阵转置和矩阵求逆。

步骤4.3：通过迭代求取目标函数的值，迭代公式为：

$x_{\left(i\right)}=x_{(i-1)}-{[B^{-1}+H_{(i-1)}^T{R}^{-1}{H}_{(i-1)}+\mathrm{\γI}]}^{-1}[H_{(i-1)}^T{R}^{-1}(y-F\left(x_{(i-1)}\right))-B^{-1}(x_{(i-1)}-x_b)]$

其中，i取正整数；x_(i)为大气状态向量第i次迭代的结果，x₍₀₎为初始迭代值；H_(i-1)为第i次迭代的雅克比矩阵，H₍₀₎根据初始迭代值x₍₀₎获得；I为单位对角阵，γ为levenberg-Marquardt参数并设定初始值，y为辐射计测量得到的大气辐射传输亮温向量。

判断迭代结果是否能接受：设当前迭代得到的大气状态向量的目标函数值为J(x_(i))，当满足条件J(x_(i-1))≥J(x_(i))时，接受本次迭代结果，否则不接受本次迭代结果；

当接受本次迭代结果时，恢复γ的值为初始值，继续步骤4.4；否则，将参数γ的值翻倍并重新计算大气状态向量x_(i)，继续判断迭代结果是否能接受；

步骤4.4：判断当前迭代是否达到如下收敛标准：

${\left[(F\left(x_{\left(i\right)}\right)-F\left(x_{(i-1)}\right))\right]}^T{S}_{\mathrm{\&delta;y}}^{-1}\left[(F\left(x_{\left(i\right)}\right)-F\left(x_{(i-1)}\right))\right]<<n$

其中，S_δy为测量值y与F(x_(i-1))的协方差矩阵；

当满足收敛标准时，结束迭代，输出最终迭代得到的大气状态向量的值；否则继续转步骤4.3进行迭代。

本发明的优点与积极效果在于：

（1）本发明方法用于地基高光谱微波辐射计，能够完整地利用高光谱辐射计多通道探测的数据，具有适应性强、计算准确的特点。

（2）本发明方法能够有效解决大气廓线反演中先验信息不足的问题，在先验信息较少，辐射计探测通道数量较多的大气廓线反演中，能完整地利用辐射计各通道的探测信息，获得较高的垂直分辨率，并具有较小的探测误差。

（3）本发明方法针对高光谱辐射计提出，能够广泛应用于通道数量很多的辐射计进行大气廓线反演，适用性较强。

附图说明

图1是本发明的大气廓线反演方法的流程示意图；

图2是本发明中湿度数据库廓线分布图；

图3是本发明中背景场协方差矩阵；

图4显示为2011年7月15日大气湿度反演的迭代过程；

图5是本发明中湿度廓线反演误差不同通道数量计算结果对比图；

图6是本发明中湿度廓线垂直分辨率不同通道数量计算结果对比图；

图7是利用本发明方法计算的BHU-K80与TP/WVP-3000湿度廓线相对误差的RMS计算结果对比图；

图8是利用本发明方法计算的BHU-K80与TP/WVP-3000湿度廓线垂直分辨率计算结果对比图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

为充分利用高光谱高信息量的特点，不损失有效的探测信息，同时获得更精确的大气廓线反演结果，本发明变分方法结合大气辐射传输正向模型和贝叶斯理论设计的大气廓线反演迭代计算方案，根据大量通道亮温反演大气廓线参数。针对高光谱辐射计，本发明方法能完整地利用其大量通道的探测信息，充分发挥高光谱辐射计探测的优势特点，获得准确的反演结果。

本发明提供的一种基于地基高光谱微波辐射计的大气廓线反演方法，是通过以下步骤实现的，流程如图1所示。

步骤一：设置地基高光谱微波辐射计的系统参数。

本发明实施例的地基微波辐射计的系统参数设置为：探测频段为K波段（18GHz-26GHz），从中以0.1GHz为间隔平均分配共80个通道，辐射计探测误差为0.5K。实际中，可根据辐射计具体系统参数情况，对地基微波辐射计的上述系统参数进行调整。在后面的具体实例的比较对照中，就给出了在不同的辐射计系统参数下，地基微波辐射计对大气廓线的反演结果精度对比，具体包括不同通道数量对比以及与TP/WVP-3000微波辐射计对比。

步骤二：获取背景参数以及背景场误差。

背景参数指被反演地区大气廓线的分布特点，它是大气廓线反演的先验数据，是对某一地区具有共同特点的大气廓线（如数年中的同一季节、或数天中的同一时刻）探测结果的统计规律计算。它在反演过程中用于给出廓线迭代的初始值和生成背景场协方差矩阵。

背景参数表示为大气状态的平均值x_b，也就是背景廓线。背景廓线一般选取探测日数天前的大气廓线。背景场误差通过一定量大气廓线反演的先验数据得到。背景参数为正向模型提供了计算基础，也为反演计算提供了背景场误差。背景场误差则反映了该大气廓线的准确性，也反映出大气廓线变化规律的统计特点，一般表示为一定数量先验信息（大气背景廓线）的协方差矩阵B。

本发明实施例的地基微波辐射计选用的大气数据来自数据库SeeBorV5。大气数据包含温度、水汽和压强信息。数据库SeeBorV5是晴朗的天空条件下进行高光谱和多光谱大气反演的全球培训数据库，广泛用于论证红外高光谱辐射计性能。在计算过程中，使用被测日第五天前的探空气球观测的大气廓线作为背景廓线，选取一定量的先验数据得到背景场误差。背景参数以及背景场协方差矩阵计算结果分别如图2和图3所示。图2中的横坐标表示为年份，背景数据库选择为从1991年1月1日到2013年3月15日，观测时间为每天上午十一点，纵坐标表示垂直高度，选择了从近地面到6km的数据进行分析，图中灰度值表示对数据库中水汽压取自然对数后的值。图3横纵坐标均为垂直海拔高度，灰度值代表数据库中相应高度层的相关系数。

步骤三：建立辐射传输模型，计算大气辐射传输亮温。

本发明实施例中选用地基微波辐射计标量辐射传输模型，如下式：

$T_B=\mathrm{\&mu;}{\&Integral;}_0^d{\mathrm{dz\&kappa;}}_a\left(z\right)T\left(z\right)\exp (-\mathrm{\&mu;}{\&Integral;}_0^z{\mathrm{dz}}^{\&prime;}{\mathrm{\&kappa;}}_a\left(z^{\&prime;}\right))---\left(1\right)$

其中，T_B为大气辐射传输亮温，参数μ＝secθ，θ为探测角度，d表示地面到大气顶层的高度，z表示地面到当前计算的大气层高度，κ_a(z)为高度z处的吸收率，T(z)为高度z处的大气温度。该方程以大气温湿度廓线为输入，通过MPM(Maximizer of the posteriorimarginal,最大后验边缘概率)模型，根据探测频率与大气廓线计算吸收率，本实施例的频段为K波段，主要为水汽吸收，因此这里的吸收率即为水汽的吸收率，再结合不同高度的温度廓线等大气状态，结合吸收率、辐射计参数计算微波辐射计接收的辐射传输亮温。输入的辐射计的系统参数和大气廓线（即背景参数）由步骤一和步骤二提供，通过辐射传输方程模拟辐射计的探测结果。

步骤三也是本发明大气廓线迭代反演的基础，为变分方法提供基本物理模型。由于迭代方法需要根据物理模型反演对象物理过程，需要使用本步骤中所使用的辐射传输模型进行计算，以求解下面步骤中的雅克比矩阵H(x)。

步骤四：利用变分方法的迭代过程反演大气廓线。图4为通过迭代进行的大气参数反演，可以发现，随着迭代次数的增加，反演的大气廓线能够逐渐接近目标廓线，受背景场干扰也越来越小。

步骤四包括以下步骤4.1～步骤4.4。

步骤4.1，建立辐射计测量的大气辐射传输亮温与待反演的大气状态的关系。

记x＝(x₁,x₂,…x_m)为描述大气状态的向量，待反演的未知量，如湿度廓线，其中m为大气状态参量的维度，表示研究范围内的高度层，x₁,x₂,…,x_m分别表示每个高度层的大气状态值；y＝(y₁,y₂,…y_n)为辐射计测量的大气辐射传输亮温，用于反演x，并假设其测量误差为ε，n为大气辐射传输亮温的维度，它可以代表辐射计的通道数目，y₁,y₂,…,y_n分别表示辐射计每个通道的大气辐射传输亮温。辐射计每个通道的大气辐射传输亮温根据公式（1）计算得到。

y与x的关系通常可以表示为前向模型F(x)：

y＝F(x)+ε (2)

向量y为包含误差的观测数据；F(x)是大气状态的向量化函数，称为前向模型。测量误差ε是一个随机变量，因此辐射计测量得到的大气辐射传输亮温也是一个随机变量。

步骤4.2，基于贝叶斯理论，确定计算待反演未知量的目标函数。

首先，假设ε服从高斯分布，则y的后验概率P(y|x)满足下式：

-2lnP(y|x)＝[y-F(x)]^TR^-1[y-F(x)]+C₁ (3)

其中，C₁是常数，R为测量误差的协方差矩阵，上角标T和-1分别表示矩阵转置和矩阵求逆。设对大气状态x进行长期观测，其测量值x相对于大气状态平均值x_b（背景廓线）的偏差也是一个随机变量，设这个随机变量的概率密度P(x)也是高斯分布，满足下式：

-2lnP(x)＝[x-x_b]^TB^-1[x-x_b]+C₂ (4)

其中C₂是常数，x_b是x的一个先验值，如气候平均值或者初始猜测值，B为背景场协方差矩阵。根据贝叶斯理论有：

$P\left(x\right|y)=\frac{P\left(y\right|x)P\left(x\right)}{P\left(y\right)}---\left(5\right)$

其中，P(x|y)是x的后验概率，P(x)是x的概率，P(y)是y的概率。

则通过以上各式有：

-2lnP(x|y)＝[x-x_b]^TB^-1[x-x_b] (6)

+[y-F(x)]^TR^-1[y-F(x)]+C₃

其中，C₃是常数。本发明方法中，待反演的状态量条件概率分布P(x|y)取得最大值时，状态量x即为反演值。所以设定目标函数J(x)如式(7)所示，则问题转化为找到目标函数J(x)＝-2lnP(x|y)取得最小值时的x，最常用的方法为另表示对函数求导。

J(x)＝[x-x_b]^ΤB^-1[x-x_b]+[y-F(x)]^ΤR^-1[y-F(x)] (7)

$\&dtri;J\left(x\right)=-{2[\&dtri;F\left(x\right)]}^T{R}^{-1}[y-F\left(x\right)]+{2B}^{-1}[x-x_b]---\left(8\right)$

记称为雅克比(Jacobian)矩阵或者核函数，它包括反演通道的大气辐射传输亮温值对每一个大气状态参量的偏微分其中y_i代表第i个辐射计通道的大气辐射传输亮温，x_j代表高度层j的大气状态参量。在实际计算时，通常对正向模型F(x)中的每个大气状态参量元素进行扰动，并查看辐射计观测值的变化来求得雅可比矩阵其中y_i是辐射计观测得到的大气辐射传输亮温，y_i,j是第j层大气状态量加入微小扰动Δx^*后通过正向模型计算得到的亮温。

通过下面方程寻找使得目标函数J(x)取得最小值时的x：

H(x)^TR^-1[y-F(x)]+B^-1[x-x_b]＝0 (9)

下面公式中的H(x)简写为H。

步骤4.3，利用迭代法求取目标函数的值。

大气廓线的反演问题可以转化为对式(9)求解，但是该方程组是病态的，未知数个数大于方程个数，无法求出其确定解，本发明采用迭代法求得近似解，基本迭代公式为：

$x_{\left(i\right)}=x_{(i-1)}-{[\&dtri;f\left(x_{(i-1)}\right)]}^{-1}f\left(x_{(i-1)}\right)---\left(10\right)$

x_(i)为大气状态（温度或者湿度廓线）向量的第i次迭代的结果，i取正整数。x₍₀₎为初始迭代值，一般选取与反演廓线观测地点和观测时间段一致的前某天的大气廓线值，例如要反演北京今天上午十一点的湿度廓线，初始值可以选取为北京昨天或者前天上午十一点的湿度廓线。中间变量 $f\left(x_{(i-1)}\right)=\&dtri;J\left(x_{(i-1)}\right),$ 于是有：

$x_{\left(i\right)}=x_{(i-1)}-{\left[{\&dtri;}^{2}J\left(x_{(i-1)}\right)\right]}^{-1}\&dtri;J\left(x_{(i-1)}\right)---\left(11\right)$

其中，令称为牛顿方向，一次迭代所产生的修正量称为一个牛顿步。目标函数的二阶导数通常称为海瑟矩阵即：

$G\left(x_{(i-1)}\right)=B^{-1}+H{\left(x_{(i-1)}\right)}^T{R}^{-1}H\left(x_{(i-1)}\right)-{[\&dtri;H\left(x_{(i-1)}\right)]}^T{R}^{-1}[y-H\left(x_{(i-1)}\right)]---\left(12\right)$

海瑟矩阵第二项中是一个非常复杂的量，很难计算。但是发现与R^-1[y-H(x_(i-1))]的乘积是一个很小的量，而且随着迭代求解的过程进行会越来越小，因此在反演算法中忽略此项。另外在迭代过程中可能会出现在某步迭代时目标函数上升，初始点与极小点较远时，产生的点列可能不收敛，或者海瑟矩阵奇异导致无法产生新的迭代点。为解决由于海瑟矩阵非正定或者奇异导致的迭代问题，对每次迭代的海瑟矩阵G(x_(i-1))进行相应的修正。令x_(i-1)为第i-1次迭代的大气状态向量，Ι为单位对角阵，只要保证γ足够大，G(x_(i-1))的特征值均为正数，这样就保证了海瑟矩阵G(x_(i-1))为正定同样也保证了其非奇异。γ称为levenberg-Marquardt参数，通过这种修正后的迭代方法称为拉凡格氏迭代法。其迭代公式为：

$x_{\left(i\right)}=x_{(i-1)}-{[B^{-1}+H_{(i-1)}^T{R}^{-1}{H}_{(i-1)}+\mathrm{\&gamma;I}]}^{-1}[H_{(i-1)}^T{R}^{-1}(y-F\left(x_{(i-1)}\right))-B^{-1}(x_{(i-1)}-x_b)]---\left(13\right)$

其中，H_(i-1)为第i-1次迭代的雅克比矩阵，H₍₀₎为根据初始迭代值x₍₀₎获得， x_b为步骤二中获取的初次的大气状态的平均值；y为微波辐射计测量的大气辐射传输亮温；F(x_(i-1))是用第i-1次迭代解通过正向模型计算得到的辐射传输亮温；R为辐射计观测值得到的大气辐射传输亮温测量误差的协方差矩阵；B为是背景场协方差矩阵。

根据当前迭代结果的目标函数值判断是否接受本次迭代。反演时γ设定相应的初始值，在迭代时根据目标函数J(x_(i))迭代前后的变化情况而适当调整。本发明实施例中，γ初始值设置为2。如果J(x_(i))增大，即J(x_(i-1))<J(x_(i))，说明新产生的迭代点不满足要求，则γ翻倍并重复当前迭代，利用公式(13)重新计算大气状态向量；如果J(x_(i))减小，即J(x_(i-1))≥J(x_(i))时，则接受这次迭代过程，γ恢复到初始值2。

当接受本次迭代结果时，恢复γ的值为初始值，继续步骤4.4；否则，继续将参数γ的值翻倍并重新计算大气状态向量x_(i)，继续迭代并判断直到迭代结果能接受。

步骤4.4，判断当前迭代是否达到收敛标准，若达到，则结束迭代，输出最终迭代得到的大气状态向量的值，否则，继续转步骤4.3迭代。

反演迭代时当满足如下收敛标准时迭代结束：

${\left[(F\left(x_{\left(i\right)}\right)-F\left(x_{(i-1)}\right))\right]}^T{S}_{\mathrm{\&delta;y}}^{-1}\left[(F\left(x_{\left(i\right)}\right)-F\left(x_{(i-1)}\right))\right]<<n---\left(14\right)$

其中S_δy为测量值y与F(x_(i-1))的协方差矩阵，n为测量值y的维度，即辐射计探测通道数目。

下面对本发明方法得到的大气状态向量的值进行误差计算。

垂直分辨率是评价反演结果的一种指标，指在空间中区别两点差异的能力，它反应了解决空间扰动的能力。而反演误差则直接反映了反演结果的精度。

垂直分辨率Δz为：

$\mathrm{\&Delta;z}=\&PartialD;z/\mathrm{diag}\left({\mathrm{AH}}^T{R}^{-1}H\right)---\left(15\right)$

其中A为反演误差：

$A={(H_{\left(K\right)}^T{R}^{-1}{H}_{\left(K\right)}+B^{-1})}^{-1}---\left(16\right)$

而为空间某高度层，K表示步骤4中总的迭代次数，H_(K)表示最后一次迭代所得到的雅可比矩阵。A的迹为每一层的信号自由度，而信号自由度可以理解为有效探测通道数目，A的倒数代表着每个自由度的高度层。一定有效探测通道数目的辐射计在探测某一特定高度的大气廓线时会受到附近其他高度的影响，其受影响高度大小即为垂直分辨率。垂直分辨率越高，敏感高度层越小，测量值更接近真实值。

为进一步说明本发明方法的探测优势，以不同通道数量的辐射计为对比，进行湿度廓线探测误差和探测垂直分辨率的分析，分别如图5和图6所示。可以发现，随着通道数量增加，大气参数反演结果更准确。

将本发明方法分别应用在TP/WVP-3000辐射计和本发明实施例中的辐射计（本发明实施例中的辐射计标记为BHU-K80）。图7是湿度廓线相对误差的均方根对比，图8是湿度廓线垂直分辨率的对比。对于每个相对误差计算方法为|x_true-x_retrieval|/x_true*100%，本发明实施例认为数据库中由湿度廓线为当前大气状态的真值，x_true表示大气状态的真值，x_retrieval表示利用辐射计测量计算得到的大气状态值。TP/WVP-3000辐射计水汽的探测通道是经过相应的通道优化设计后选取的，为：22.235GHz,23.035GHz,23.835GHz,26.235GHz以及30.000GHz，测量误差与本发明实施例的辐射计一致，均为0.5K，在同等条件下本发明实施例的80通道的微波高光谱辐射计较湿度探测通道数目为5的辐射计TP/WVP-3000，反演性能有了较大的提高。在近地表面，辐射计TP/WVP-3000的湿度廓线相对误差均方根为22%，而本发明实施例的辐射计有明显的了提高为12%，提高幅度为10%，同样的在海拔高度为5km时，由TP/WVP-3000的15%提高到10%。从图8中可以看出垂直分辨率随着高度的增加而减小，垂直海拔高度在0～1km范围内，辐射计TP/WVP-3000与本发明实施例的辐射计的垂直分辨率差别不大，但是在1～6km范围内，本发明实施例的微波高光谱辐射计具有更好的垂直分辨率，并且在高海拔处有着更明显的优势。根据上述模拟的结果，可以发现在相同的条件下，利用本方法实现高光谱辐射计的大气廓线反演，能够获得更精确的大气廓线，更好的湿度廓线反演性能与垂直分辨率。

Claims (3)

1.一种基于地基高光谱微波辐射计的大气廓线反演方法，其特征在于，包括如下步骤：

步骤1：设置地基高光谱微波辐射计的系统参数；所述系统参数包括探测频段，通道数量以及各通道频率和辐射计的探测误差；

步骤2：获取背景参数以及背景场误差；所述背景参数表示为大气状态的平均值x_b，通过统计被反演地区探测日前具有共同特点的历史大气廓线探测结果得到，是大气廓线反演的先验数据；所述背景场误差反映了被反演地区大气廓线的准确性，表示为通过先验数据得到的背景场协方差矩阵B；

步骤3：建立辐射传输模型，确定大气辐射传输亮温；

所采用的辐射传输模型如下：

$T_B=\mathrm{\&mu;}{\&Integral;}_0^d{\mathrm{dz\&kappa;}}_a\left(z\right)T\left(z\right)\exp (-\mathrm{\&mu;}{\&Integral;}_0^z{\mathrm{dz}}^{\&prime;}{\mathrm{\&kappa;}}_a(z^{\&prime;}\left)\right)$

其中，T_B为大气辐射传输亮温；参数μ＝secθ，θ为探测角度；d表示地面到大气顶层的高度；κ_a(z)为高度z处的吸收率；T(z)为高度z处的大气温度；

步骤4：利用变分方法的迭代过程反演大气廓线，包括如下子步骤：

步骤4.1：建立辐射计测量的大气辐射传输亮温与待反演的大气状态的关系；

记x＝(x₁,x₂,…x_m)为描述大气状态的向量，为待反演的未知量，m为被反演地区研究范围内的高度层，x_m为m高度层的大气状态值；记y＝(y₁,y₂,…y_n)为辐射计测量的大气辐射传输亮温的向量，n为辐射计的通道数目，y_n为第n个通道的大气辐射传输亮温；

采用前向模型F(x)表示y与x的关系：y＝F(x)+ε，ε为测量误差，是一个随机变量；

步骤4.2：确定计算待反演未知量的目标函数；

所述的目标函数J(x)如下：

J(x)＝[x-x_b]^ΤB^-1[x-x_b]+[y-F(x)]^ΤR^-1[y-F(x)]

通过下面方程寻找使得目标函数J(x)取得最小值时的x：

H^TR^-1[y-F(x)]+B^-1[x-x_b]＝0

其中，H为雅克比矩阵H(x)的简写，H(x)＝▽F(x)，包括辐射计通道的大气辐射传输亮温对每个大气状态的偏微分；R为测量误差ε的协方差矩阵；上角标T和-1分别表示矩阵转置和矩阵求逆；

步骤4.3：通过迭代求取目标函数的值，迭代公式为：

$x_{\left(i\right)}=x_{\left(i-1\right)}-{\&lsqb;B^{-1}+H_{\left(i-1\right)}^T{R}^{-1}{H}_{\left(i-1\right)}+\mathrm{\&gamma;}I\&rsqb;}^{-1}\&lsqb;H_{\left(i-1\right)}^T{R}^{-1}\left(y-F\left(x_{\left(i-1\right)}\right)\right)-B^{-1}\left(x_{\left(i-1\right)}-x_b\right)\&rsqb;$

其中，i取正整数；x_(i)为大气状态向量第i次迭代的结果，x₍₀₎为初始迭代值；H_(i-1)为第i-1次迭代的雅克比矩阵，H₍₀₎根据初始迭代值x₍₀₎获得；I为单位对角阵，γ为levenberg-Marquardt参数并设定了初始值；

判断迭代结果是否能接受：设当前迭代得到的大气状态向量的目标函数值为J(x_(i))，当满足条件J(x_(i-1))≥J(x_(i))时，接受本次迭代结果，否则不接受本次迭代结果；

当接受本次迭代结果时，恢复γ的值为初始值，继续步骤4.4；否则，将参数γ的值翻倍并重新计算大气状态向量x_(i)，继续判断迭代结果是否能接受；

其中，所述的迭代公式通过如下步骤确定：

(1)确定基本迭代公式为：

x_(i)＝x_(i-1)-[▽f(x_(i-1))]^-1f(x_(i-1))

(2)设定中间变量f(x_(i-1))＝▽J(x_(i-1))，得到：

x_(i)＝x_(i-1)-[▽²J(x_(i-1))]^-1▽J(x_(i-1))

(3)设定海瑟矩阵G(x_(i-1))＝▽²J(x_(i-1))，得到：

G(x_(i-1))＝B^-1+H(x_(i-1))^TR^-1H(x_(i-1))-[▽H(x_(i-1))]^TR^-1[y-H(x_(i-1))]

忽略项[▽H(x_(i-1))]^TR^-1[y-H(x_(i-1))]，

修正海瑟矩阵G(x_(i-1))，令G(x_(i-1))＝(▽²J(x_(i-1))+γΙ)，

进而获得所述的迭代公式；

步骤4.4：判断当前迭代是否达到如下收敛标准：

${\&lsqb;\left(F\right(x_{\left(i\right)})-F(x_{(i-1)}\left)\right)\&rsqb;}^T{S}_{\mathrm{\&delta;}y}^{-1}\&lsqb;\left(F\right(x_{\left(i\right)})-F(x_{(i-1)}\left)\right)\&rsqb;<<n$

其中，S_δy为测量值y与F(x_(i-1))的协方差矩阵；

当满足收敛标准时，结束迭代，输出最终迭代得到的大气状态向量的值；否则继续转步骤4.3进行迭代。

2.根据权利要求1所述的一种基于地基高光谱微波辐射计的大气廓线反演方法，其特征在于，步骤4.2中所述的目标函数，通过如下方法确定：

首先，设ε服从高斯分布，则观测数据y的概率分布P(y|x)满足下式：

-2ln P(y|x)＝[y-F(x)]^TR^-1[y-F(x)]+C₁

设测量的大气状态x相对于大气状态的平均值x_b的偏差是一个随机变量，设这个随机变量的概率密度P(x)也是高斯分布，满足下式：

-2ln P(x)＝[x-x_b]^ΤB^-1[x-x_b]+C₂

其中，C₁、C₂都是常数；

根据贝叶斯理论，得到下式：

$P\left(x\right|y)=\frac{P\left(y\right|x)P\left(x\right)}{P\left(y\right)}$

P(x|y)是x的后验概率，P(x)是x的概率，P(y)是y的概率，则进一步得到：

-2ln P(x|y)＝[x-x_b]^ΤB^-1[x-x_b]+[y-F(x)]^ΤR^-1[y-F(x)]+C₃

其中，C₃是常数；

当概率P(x|y)取得最大值时，状态量x为反演值，也就是目标函数J(x)＝-2lnP(x|y)，获取目标函数J(x)取得最小值时的x。

3.根据权利要求1所述的一种基于地基高光谱微波辐射计的大气廓线反演方法，其特征在于，步骤4.3中所述的参数γ，设定初始值为2。