沙尘天气下黑河流域大气长波辐射遥感估算

王子超， 王春磊， 马俊俊

（1.华北理工大学人工智能学院，河北 唐山 063210；2.自然资源部咨询研究中心，北京 100100）

作为地-气系统能量交换的重要表征因子，近地表大气长波辐射（Downward surface longwave radiation，DSLR）在能量收支和气象气候变化研究中不可或缺。DSLR 是整个垂直柱中大气散射、吸收和发射的结果，可以通过复杂的辐射传输模型（如MODTRAN）利用大气成分垂直分布状态、温湿度廓线、气溶胶和云参数等进行估算［1-3］。但在实际使用中，由于辐射传输模型的输入参数较多且难以获取，容易引起误差传递。且由于DSLR 一半以上来源于地表100 m 内的大气发射辐射，受近地表大气温度梯度影响很大。受制于仪器测量精度和局地地理信息状况，当研究者通过地面观测资料特别是大气温压数据参数化大气比辐射率来进行DSLR 估算时［4-7］，模型稳定性和普适性较差。近年，部分学者提出利用卫星反演的时空连续性和多层次性的地表定量物理参数替代建模时需要的部分参数估算DSLR［8］。也有学者基于辐射传输模拟数据建立大气层顶（Top of the atmosphere，TOA）辐亮度或亮温数据与DSLR 数据集，利用透过率和气压的关系代数式作为通道权重函数描述各层大气辐射对长波下行的贡献，建立DSLR 和传感器辐亮度数据之间的统计关系模型［9-10］。后者具有一定的物理基础，且简单快捷，算法具有较高的稳定性，适合遥感数据产品的业务化运行，但未考虑云和气溶胶类型对模型的影响。

沙尘气溶胶作为重要的辐射强迫因子，占对流层气溶胶总量的半数之多［11］，沙尘气溶胶影响下的长波辐射估算是一项具有挑战性的工作，沙尘光学性质和辐射传输机理的复杂性是主要误差来源。近年已有大量学者就沙尘气溶胶开展研究，茹建波等［12］尝试从气溶胶总的后向散射系数中提取沙尘气溶胶消光系数，谢艳清［13］、贾臣［14］、王联霞等［15］分别利用FY-4A、Landsat8、VIIRS 等开展星载传感器反演气溶胶光学厚度的研究，李丁等［16］利用S5P/TROPOM 数据对中国东部气溶胶单次散射反照率进行反演。相对成熟的气溶胶光学参数反演方法为研究其辐射效应奠定基础，Xin 等［17］利用便携式太阳光度计测量了中国西北4个沙漠半沙漠地区的沙尘气溶胶光学厚度和不对称参数，发现其光学特性和辐射特性有强烈的季节变化，气溶胶的TOA 辐射强迫在3.9～12 W·m-2左右；田磊等［18］利用地基观测资料计算分析了沙尘气溶胶对大气逆辐射的影响，当大气浑浊度小于0.3 时，大气逆辐射随大气浑浊度的增加而增加，超过临界值时大气逆辐射呈现相反的趋势；Sicard 等［19］基于辐射传输模型对沙尘气溶胶的长波辐射强迫进行研究，发现沙尘粒子的半径、折射率、分布状态和表面温度等对辐射强迫敏感性较高。前人的研究集中在对沙尘气溶胶参数的定量提取及其引起的辐射效应描述和分析上，对沙尘气溶胶条件下的DSLR遥感估算研究较少。

为此，本文通过分析不同地表参量和大气廓线场景下沙尘气溶胶变动对DSLR 和MODIS通道辐亮度的影响，利用辐射传输模拟数据建立TOA 辐亮度和DSLR 的统计关系模型，发展新的DSLR 遥感反演方法，最终利用黑河流域地面实测数据对模型进行验证和普适性评价，以提高沙尘气溶胶条件下DSLR遥感反演精度，并为相关研究提供参考。

1 数据与方法

1.1 研究区概况

本文选取黑河流域作为研究区，地理位置介于98°～102°E、38°～42°N 之间，面积约1.429×105km2。黑河流域上、中、下游区域涵盖草地、农田、湿地等多种地表类型，是我国陆地河流研究的重要基地。围绕该区域已开展大量的科学研究实验［20-21］，积累了大量水文气象观测网数据，既为流域生态水文研究奠定了很好的科学数据基础，亦可为本文研究沙尘气溶胶模型提供帮助。研究区范围和站点分布如图1所示。

图1 研究区范围和站点分布Fig.1 Distribution map of study area and sites

1.2 数据来源

1.2.1 模拟数据考虑到大气温湿度、风以及气溶胶参量等因素处于动态变化中，地面测量数据与卫星同步观测数据难以达到同步，本文利用辐射传输模型MODTRAN生产的数据集分析和评估沙尘气溶胶对DSLR 遥感的影响。其中，大气廓线来源于热力学初始分析资料（Thermodynamic Initial Guess Retrieval，TIGR）数据库，每条廓线含40 层温度、湿度、压力、高程等信息。考虑到晴空和沙尘爆发的场景因素，以相对湿度和连续湿度作为筛选条件选出946 条晴空大气廓线，大气水汽含量（Water vapor content，WVC）范围为0.05～6.15 g·cm-2，底层廓线温度范围为235～314 K，基本囊括了研究区内四季常见的廓线类型；发射率来自ASTER 地表发射率库，针对黑河流域地表覆盖类型特点，选择7 种主要地表覆盖类型进行模拟。气溶胶对辐射的影响取决于其空间分布、自身的物理化学性质（包括粒子尺度、谱分布、化学成分等）以及光学特性等（如光学厚度、单次散射反照率），沙尘气溶胶相应的信息取自云和气溶胶粒子光学特性软件（Optical properties of aerosols and clouds，OPAC）的参数库［22］。为增强普适性，基于底层大气温度对地表温度按一定区间变动，变动区间范围为-5～15 K，步长为5 K。

1.2.2 遥感数据本文利用Aqua 卫星上的MODIS数据产品来验证模型的可行性，采样时间点为2018年1 月至2019 年12 月。主要使用6 种类型的MODIS 产品，分别为：星上辐亮度产品MYD021KM、几何定位产品MYD03、水汽数据产品MYD05_L2、云掩膜产品MYD35_L2、地表温度和发射率产品MYD11_L2以及气溶胶产品MYD04(C6.1)。

MYD021KM 是1 级辐亮度产品，包含经过辐射校正大气顶部的反射率和辐射亮度数据。MYD03是地理定位的1 级产品，内含每个1 km 分辨率像元的地理坐标等信息，大地经纬度以及太阳天顶角和观测角、卫星观测天顶角（View zenith angle，VZA）和方位角、海陆掩膜等重要信息。

MYD05_L2 是大气可降水量2 级产品，它是在表面反射率线性光谱变化的假设下，基于多个不同吸收通道反演水汽量的加权平均，其中每个通道的值通过查找表法得到。MYD35_L2 是大气2 级标准数据产品，内容为云掩膜，可以用于云检测并进行晴天选取。MYD11_L2 是陆地2 级标准数据产品，其中包含1 km地表温度和发射率值。

MYD04(C6.1)是MODIS Collection 006 版气溶胶产品，该产品融合深蓝算法和暗目标算法的合成结果，并在气体吸收计算、瑞利散射计算、陆地/海洋判识、云检测和产品质控方面做了较多改进［23］，空间分辨率为3 km。

为增强模型验证的可信度，需要在时间和空间上将MODIS 观测数据与地面站进行匹配。空间上，为保障数据的有效性，匹配过程中站点对应像元周边至少需要有10 个有效像元。时间上，地面站点DSLR 数据为利用卫星过境前后15 min 的测量值双线性内插得到卫星过境时刻所对应的值，由于DSLR 在短时间内变化不大，该值近似代表观测站点所在MODIS像元对应的DSLR 值。为了消除云的影响，本文只选取晴空数据。为确保选择的像元达到要求，仅选择MODIS 温度产品MYD11_L2 质量控制集中标记为晴空以及MYD35_L2云掩膜数据判断为晴空的像元；另外对像元是否存在气溶胶进行判断，选择站点所在位置的中心像元及周边存在气溶胶的像元作为符合要求的目标像元。

1.2.3 站点数据本研究选取黑河流域2018年1月至2019 年12 月的站点数据对模型进行检验。数据来源于国家青藏高原科学数据中心提供的黑河流域生态-水文过程综合遥感联合试验站点［24］。以花寨子荒漠站、混合林站、黑河遥感站、张掖湿地站4个不同下垫面气象站（表1）数据进行验证。这些数据是由地面辐射计每隔10 min 观测一次所得，利用卫星过境前后15 min 的测量值双线性内插得到卫星过境时刻所对应的DSLR 测量值，该值近似代表观测站点所在MODIS像元对应的DSLR值。

表1 气象站基本信息Tab.1 Basic information of meteorological stations

1.3 研究方法

1.3.1 近地表大气长波辐射线性模型算法Wang等［10］基于透过率和气压的关系代数式作为通道权重函数描述各层大气辐射对长波下行的贡献，建立了大气长波辐射估算的混合模型法。该算法机理明确，能弥补物理模型中由于辐射传输过程过于复杂或大气输入数据难以获取所带来的缺陷，且直接利用TOA 辐亮度估算DSLR，算法计算效率高，稳定性和鲁棒性较好，适合业务化运行［9-10］。该模型形式如下：

式中：DSLRC为晴空条件下大气下行长波辐射估算值；a0～a7为依赖于VZA、WVC 以及地表高程的回归系数；L27～L29和L31～L34为MODIS 热红外通道测得的TOA 辐亮度。此模型可以解释大气长波辐射数据92%以上的变化，在晴空MODIS 的数据应用中，模型的标准误差约16 W·m-2，但其主要应用于稳定气溶胶类型且气溶胶含量较小的晴空大气环境中，在沙尘气溶胶天气环境中的应用能力还未得到检验。本研究拟以此模型为基础，增加相应的气溶胶参数，拓展模型的应用能力。

1.3.2 算法修正思路沙尘天气下DSLR 的反演难点在于难以量化气溶胶对大气的长波辐射强迫和由于气溶胶变化造成的通道TOA 辐亮度的变动。为研究沙尘气溶胶的长波辐射效应，本研究首先分析了4 种不同气溶胶类型引起的辐射强迫效应，在此基础上利用气溶胶光学厚度（Aerosol optical depth，AOD）对MODIS 第29、31、32 通道进行星上辐亮度校正、优化系数，构建沙尘气溶胶天气下近地表大气长波辐射遥感模型。

图2 为沙尘、乡村、城市和海洋4 种气溶胶存在时的气溶胶长波辐射强迫（Aerosol longwave radiative forcing，ALRF）均值偏差（Bias）和标准差（Standard deviation，STD）随AOD 变动的关系图。从图中可以看出，4 种气溶胶的ALRF 及其标准差随AOD的增大而增大，ALRF与AOD 之间近似呈线性关系。其中，海洋气溶胶的ALRF 最小，在AOD=1.0时平均为3.85 W·m-2，其STD 为1.1 W·m-2；沙尘气溶胶造成的ALRF 最大，随AOD 变化全部为正强迫，ALRF均值从2.09 W·m-2到29.8 W·m-2，其STD 变化为0.7～8.3 W·m-2。

图2 在不同AOD情况下4种气溶胶类型的ALRF均值及STDFig.2 Average ALRF and STD of four different aerosol types varies with the AOD

沙尘气溶胶条件下的MODIS 热红外TOA 通道辐亮度变化情况如图3 所示，TOA 通道辐亮度变化率的Bias、STD和AOD同样呈线性关系。整体来说，TOA辐亮度变化率Bias较小，但通道29、31和32这3个地表温度通道的STD 相对较大，分别达到3.51%、3.14%和3.32%，说明在改进DSLRC模型，使之能适用于沙尘气溶胶条件下的DSLR 估算时，利用AOD进行辐亮度校正并进行系数率定的思路是可行的。

图3 沙尘气溶胶条件下不同AOD情况下TOA通道辐亮度变化率及STDFig.3 TOA channel radiance change rate and STD of different AOD conditions under dust aerosol conditions

在此基础上，本文利用AOD 对MODIS 第29、31、32 通道进行星上辐亮度校正，并优化系数，沙尘气溶胶天气下近地表大气长波辐射遥感模型（DSLRA）形式为：

式中：a0～a7与b0～b2为由VZA、AOD 和WVC 的回归系数。为使算法更具代表性和减少输入参数本身带有的误差，采用在WVC 和AOD 相邻区间设置部分重叠区间并分组细化的方法来优化系数。其中，VZA分为6种情况：0、33.56°、44.42°、51.32°、56.25°、60.00°；WVC（g·cm-2）分12 组：［0.0,1.0］、［0.5,1.5］、［1.0, 2.0］、［1.5, 2.5］、［2.0, 3.0］、［2.5, 3.5］、［3.0,4.0］、［3.5, 4.5］、［4.0, 5.0］、［4.5, 5.5］、［5.0, 6.0］、［5.5, 6.5］；AOD 分3 组：［0.0, 0.4］、［0.3, 0.7］、［0.6,1.0］，依上分组，采用最小二乘法共计回归216 组系数，算法构建的总体流程如图4。

图4 修正算法反演DSLR的总体流程图Fig.4 Flow chart of DSLR inversion by using the improved model

2 结果与分析

2.1 敏感性分析

为评估输入变量的不确定性所导致的DSLR 反演误差，需要对模型进行敏感性分析。模型误差主要由MODIS仪器噪声和校准误差、AOD 的不确定性和模型固有误差所决定。模型总体误差e(ΔDSLR)可表示为：

根据MODIS 的等效噪声误差在通道27、28、33和34 为0.25 K，在通道29、31 和32 为0.05 K；对于WVC，由于MODIS 产品中其精度为0.5 g·cm-2［26］，在VZA=60°观测情况下本文得到模型的总体误差（图5）。

从图5 可以看出，VZA=60°观测情况下，由通道辐亮度L和AOD 不确定性及其造成的δ(LDSLR)和δ( AOD )随着AOD 的增大而增加，但δ(alg)随着AOD的增大而减小。由于算法本身误差在整个误差中占比大，e( ΔDSLR )随着AOD 增大出现先减小后增大的趋势，最小值出现在AOD=0.4 处。当AOD=1.0，e( ΔDSLR )达到10.7 W·m-2。在沙尘气溶胶条件下，精确的AOD 值对DSLR 反演非常重要。由于不同输入参数的误差可能相互抵消，实际误差可能小于图5中的理论估计值。从其他角度进行的研究也得出了类似的结论，且e( ΔDSLR )随着观测天顶角的增大而增加。当VZA=0，DSLR 均方误差最大为9.6 W·m-2。

图5 VZA=60°时不同AOD情况下的误差Fig.5 Errors under different AOD cases when VZA=60°

2.2 检验结果分析

最终筛选出在花寨子荒漠站、混合林站、黑河遥感站、张掖湿地站的匹配数据为28、33、28、29组。花寨子荒漠站AOD 的变化范围为0.06～0.83，WVC变化范围为0.12～2.84 g·cm-2；混合林站AOD 的变化范 围 为0.09～0.86，WVC 变 化 范 围 为0.29～3.09 g·cm-2；黑河遥感站AOD 的变化范围为0.09～0.89，WVC 变 化 范 围 为0.33～3.02 g·cm-2；张 掖 湿 地 站AOD 的 变 化 范 围 为0.09～0.87，WVC 变 化 范 围 为0.44～3.43 g·cm-2。图6 为 卫 星 过 境 时 两 种 模 型DSLR反演值与实测值散点图。

图6 沙尘气溶胶条件下两种模型的DSLR反演值与实测值散点图Fig.6 Scatter plots of DSLR inversion values and measured values of the two models under the condition of dust aerosol

在花寨子荒漠站，沙尘气溶胶天气多发生在春季，DSLRA模型估算的均方误差（RMSE）为20.3 W·m-2，偏差（Bias）为-5.8 W·m-2，误差范围为-47.5～31.9 W·m-2；在混合林站，沙尘气溶胶天气多发生在春季，DSLRA模型估算的RMSE 为19.5 W·m-2，Bias为-7.8 W·m-2，误差范围为-45.3～32.9 W·m-2；在黑河遥感站，沙尘气溶胶天气多发生在春季，DSLRA模型估算的RMSE 为17.1 W·m-2，Bias 为-1.8 W·m-2，误差范围为-36.5～29.3 W·m-2；在张掖湿地站，沙尘气溶胶天气多发生在春季，DSLRA模型估算的RMSE 为20.4 W·m-2，Bias 为-12.3 W·m-2，误差范围为-46.9～18.6 W·m-2。为了更好地评价两种模型在沙尘气溶胶天气下的应用能力，模型估算误差关于AOD和WVC的散点图如图7和图8所示。

图7 沙尘气溶胶条件下两种模型的DSLR反演误差与AOD散点图Fig.7 Scatter plots of DSLR inversion error and AOD for two models under dust aerosol conditions

图8 沙尘气溶胶条件下两种模型的DSLR反演误差与WVC散点图Fig.8 Scatter plots of DSLR inversion error and WVC for two models under dust aerosol conditions

DSLRC模型实际数据验证结果比模拟数据验证结果误差大，主要因为模型未考虑到ALRF 以及没有对沙尘气溶胶情况下受影响的星上辐亮度值进行必要的系数调整。沙尘气溶胶不只产生较强的ALRF，其对星上辐亮度影响也较大，特别是对通道29、31 和32，最大变动幅度达到3.51%。从图6 中可以得出改进的算法在轻度和重度沙尘气溶胶条件下表现效果都比较理想，由于气溶胶对星上辐亮度有影响并对其进行了误差修正，使得模型在一些高水汽的廓线中也能表现良好。

但也明显看到，模型误差在实测数据中比在模拟数据中大的多。原因是多方面的：在验证过程中卫星数据是整个区域成像时刻的观测结果，而地面站点仅仅是一个点上的观测数据，这个站点的测量存在系统误差以及随机误差；本文估算的沙尘天气下DSLR 是瞬时值，且研究是在假定DSLR 短时间不变化的前提下进行的，在采样点时间前后有风或云也会使DSLR 产生较大偏差。另外，有研究表明，土地覆盖类型特别是作物生长阶段的发射率变动对DSLR 反演精度也有显著影响［27］，但在实际模拟中无法考虑这种因素也是误差较大的一个因素。沙尘气溶胶颗粒表面不均匀性和角度效应也不容忽视，而在模拟过程中我们假定气溶胶为球形，但实际气溶胶散射特性为非球面，也会影响验证结果。

除了以上因素外，AOD、WVC 和TOA 辐亮度的不确定性也需要着重考虑。AOD 具有空间异质性，受地形、湿度、温度等的影响。在垂直分布上，AOD和WVC各层面上的散射和吸收均不同，其垂直廓线上的变动是造成许多气候模型预测不准确的重要因素。利用激光雷达探测大气参数的时空分布和描述气溶胶特性是提高模型精度可行的方法。但这几个站点没有测得对应时刻的垂直廓线数据，这对模型的验证及结果分析有一定影响。尽管OPAC提供了沙尘气溶胶的光学特性，但这与站点中实际的气溶胶光学特性还是有差异的，比如源地、气溶胶成分、粒径分布及其时空变化规律不同［28］，也会对模型检验的结果产生影响。另外，云天像素被误认为晴空像素会使TOA 辐射亮度大幅度降低也是本文不可忽略的问题。站点符合要求的数据不多，模型实地验证的精度可能受个别样本点的干扰。

3 结论

（1）该算法建立在以AOD校正部分敏感通道并进行系数率定的基础上，综合考虑了AOD、WVC 和VZA等因素，在模拟数据垂直观测中可将e( ΔDSLR)控制在9.6 W·m-2。敏感性分析结果表明，e( ΔDSLR)随着AOD 增大出现先减小后增大的趋势，最小值出现在AOD=0.4处。

（2）在黑河流域4 个站点的验证中，改进的模型RMSE 在17.1～20.4 W·m-2，Bias 介于-12.3～-1.8 W·m-2，基本满足实际应用需求，但站点检验整体误差比模拟数据中偏大。

总体来说，相比DSLRC模型，在WVC 和AOD 及观测条件已知的情况下，由于考虑到了ALRF，本文提出的模型能显著提高沙尘气溶胶条件下的DSLR估算精度。然而，限于目前的研究条件和数据，一些方面有待改进，尤其是气溶胶垂直分布剖面对DSLR 反演的影响有待进一步研究。未来的研究将考虑气溶胶廓线变化对DSLR 的影响，并利用激光雷达和太阳光度计数据对气溶胶进行分析研究，研究气溶胶的垂直分布廓线，定量分析廓线变化引起的气溶胶消光后向散射比变化及其带来的消光系数误差对DSLR 反演的影响；增加站点数据并开展有关气溶胶影响下DSLR 反演模型的普适性研究；深入开展土地覆盖类型变化对DSLR 反演精度影响的研究。