应王敏 ，吴 骅 ※

（1. 中国科学院地理科学与资源研究所/资源与环境信息系统国家重点实验室，北京100101；2. 中国科学院大学/资源与环境学院，北京100049）

0 引言

地表短波净辐射（Net Surface Shortwave Radiation，NSSR），指地表吸收的短波辐射通量（0.3～5 μm）。地表短波净辐射是地表辐射平衡中的一个重要组分，是气候形成及气候变化的主要依据，影响着物质输送、地气能量交换、水汽流动、生物光合作用等［1-2］。根据地表短波净辐射遥感反演方法的不同，可以分为基于统计经验关系的方法、基于物理过程的方法、基于物理意义的参数化方法以及混合模型方法。

基于统计经验关系的方法是指不考虑太阳、大气与陆地之间复杂大气福射传输物理机制过程，直接建立地面短波辐射的测量值与遥感观测值的经验关系或统计关系，以此来估算地面的短波辐射物理量［3-4］。Pink 和Corio［5］通过研究发现，NOAA5 卫星的大气顶层观测数据与地表净辐射值具有很强的相关性，这种关系随后也被许多学者所证实［6-7］。经验统计法具有步骤简单的优点，但由于其受到复杂大气状况和复杂地表特征的限制，普适性较差。

基于物理过程的方法指根据大气辐射传输理论，给定大气和地表参数，通过求解复杂的辐射传输方程计算得到地表短波辐射［8］。根据太阳辐射经大气到达地表的传输过程，Joseph 等［9］提出了一种快速但精确的计算大气吸收—散射辐射通量的方法，Delta-Eddington 近似法，通过联合Dirac delta 函数及二流近似，解决了Eddington 近似处理高度不对称相函数精度不高的问题。Lenoble［10］讨论了不同的大气多次散射辐射传输方法，还比较模型的复杂性及精度。该方法的优点是具有明确的物理意义，计算精确。但其由于输入参数众多，计算复杂，难以直接应用于像元级运算。

基于物理意义的参数化方法则模拟了大气与地表间辐射传输的物理机制，包括大气反射、散射和吸收作用、地表大气交互作用等［11］，这些相互作用在地表短波辐射反演中进行的参数化构建［12］。例如，Pinker 和Ewing［13］对大气层中水汽、云、气溶胶的吸收与散射作用进行参数化处理，成功反演得到地表短波净辐射，并分析出太阳天顶角（SZA）对地表短波净辐射影响较大。基于物理意义的参数化方法具有一定的物理意义，能够较为精确地估算模拟结果，适用性较强。但其模型构建过程相对复杂，模型精度受大气地表输入参数的精确度和获取困难度的影响。此外，大部分参数化方案是基于比较理想的前提假设下提出的，部分参数的取值具有区域性和经验性。

混合模型法是将经验统计法和物理模型参数化法结合使用，既考虑了大气辐射传输方程的物理机制，又发挥了经验统计方法的间接操作和较高精度的优势。例如Tang 等［14］利用辐射传输模型MODTRAN 模拟生成各种复杂情况下的短波净辐射，通过统计回归的方式构建大气顶层观测数据与地表短波净辐射的非线性关系，以此来实现MODIS/TERRA 地表短波净辐射的遥感反演。该类方法具有一定物理意义，反演精度高，模型构建相对简洁等优势，是目前反演地表短波净辐射广泛应用的方法之一［15-17］。

文章基于混合模型法，试图在大量MODTRAN 模拟集上构建针对MODIS/AQUA 数据的地表短波净辐射反演模型，并利用具有不同下垫面情况的多站点长时间序列数据对反演模型进行综合全面的鲁棒性分析和精度检验，最终发展出具有一定泛化能力的适用于MODIS/AQUA 数据的高精度地表短波净辐射反演模型。

1 研究区与数据

1.1 研究区概况

地表辐射观测网（SURFRAD）是一个常用的地表辐射观测网络，由美国国家海洋和大气局提供（http：//www.srrb.noaa.gov/surfrad）。该地表辐射观测网络现有7 个观测站点，其空间分布情况如图1 所示，分别位于蒙大拿州（Fort Peck，MT；北纬48.31°，西经105.10°）、卡罗拉多州（Table Mountain，CO；北纬40.13°，西经105.24°）、伊利罗伊州（Bondville，IL；北纬40.05°，西经88.37°）、密西西比州（Goodwin Creek，MS；北纬34.25°，西经 89.87°）、宾夕法尼亚（Penn State，PA；北纬 40.72°，西经 77.93°）、内华达州（Desert Rock，NV；北纬36.62°，西经116.02°）和南达科他州（Sioux Falls，SD；北纬43.73°，西经96.62°）。站点位置选择考虑了美国的不同气候区特点、下垫面类型的异质性，相应的测量值能够有效地代表该区域的特征。SURFRAD 观测网络为美国本土的气候、数值预报及卫星遥感监测等研究提供准确、连续、长期的地面辐射观测数据，包括地表短波下行辐射、地表短波上行辐射、地表短波净辐射等，其采样的时间分辨率为1 min（2009 年1 月1 日之前的数据记录间隔为3 min）。此外，SURFRAD 还提供了质量控制（QC）、大气温度、相对湿度、风速、风向、气压等辅助变量。

图1 SURFRAD 地面实测站点分布Fig.1 The distribution of SURFRAD ground stations

1.2 MODIS/AQUA 数据

搭载于美国AQUA 太阳同步极轨卫星平台上的中分辨率成像光谱仪（MODIS），拥有36 个通道，覆盖了可见光—红外整个波谱［18］。MODIS/AQUA 在地方时下午13∶30 过境，提供了相应MYD 系列产品，空间分辨率一般为1 km。NASA 研究人员在原始辐亮度数据产品的基础上生成了很多种不同级别、不同类型的遥感应用产品（陆地产品、大气产品、海洋产品等），提供具有时空连续性的地表监测资料，为地表短波辐射的研究提供有效的手段，满足地球科学多方面的研究需求［19］。

该文所用到的产品包括MYD02，MYD03，MYD05 和MYD35，下载了过境研究区域7 个地面站点的2017 年全年遥感影像（https：//ladsweb.modaps.eosdis.nasa.gov/）。其中MYD02 产品可以提供通道辐射亮度信息，MYD03 产品提供了观测几何参数，MYD05 产品可提供大气可降水量变量，而MYD35 产品则提供相应的云掩膜信息，可用于判断晴空和有云条件。由于受仪器自身硬件的影响，MYD02 的第6 个通道存在部分数据缺失的现象，在大气层的遥感影像上表现出较多的条纹，因此该文仅利用了MYD02 的前5 个通道和第7 个通道来反演地表的短波净辐射。

1.3 MODTRAN 模拟数据

大气辐射传输模型MODTRAN 5 可以有效模拟各种典型地表、大气剖面、气溶胶、云、观测几何等复杂条件下的地表短波辐射。基于MODTRAN 的模拟就能获取不同地表和大气状况下的模拟数据，从而构建地表短波辐射反演模型。

设置模拟的波谱范围为300～5 000 nm，光谱分辨率为2 波数（cm-1）。为了充分模拟各种复杂情况，从MODTRAN 默认数据库里选择了9 种有代表性的地类波谱反射率数据作为输入，代表着不同下垫面类型，包括草地、湿地、土壤、森林、沙漠、城市、海水、新雪、海冰。在模拟时，采用了MODTRAN 5 提供的6 种标准大气廓线（热带大气，中纬度夏季大气，中纬度冬季大气，副极带夏季大气，副极带冬季大气以及美国1976 年标准大气US76）。同时，使用了3 种气溶胶类型（乡村气溶胶、海洋性气溶胶以及对流层气溶胶），相应的能见度设置为缺省值（能见度取值范围为10～35 km）。除了晴空天气以外，还考虑了3 种云模型，即层云、高层云和积云。考虑到一般MODIS 数据的最大观测天顶角小于 65°，故模拟了 6 种观测天顶角 VZA（0°，33.5°，44.4°，51.3°，56.2°，60°）。此外，考虑了36 种太阳天顶角SZA 变化，即0°～70°，间隔2°。将这些设置参数条件排列组合，总共模拟了139968 个样本量的模拟数据集。

MODTRAN 模型可以模拟相应情况下的波谱辐射亮度（存储在TP7 文件里）、多层高度对应的向上辐射通量和向下辐射通量（存储在FLX 文件里）等。地表短波净辐射（NSSR）可以通过将FLX 文件中最底层高度对应的向下辐射通量减去向上辐射通量得到；MODIS 通道辐射亮度（Li）可以将TP7 文件中波谱辐射亮度与MODIS 通道响应函数卷积得到。

2 地表短波净辐射反演方法

利用大气顶层TOA 向上辐射通量（Fu）与地表短波净辐射（NSSR）的非线性关系可有效反演NSSR［14］，该方法已经被广泛应用到各种卫星传感器，并得到了较好的反演结果。公式（1）是该方法的核心思想，NSSR 与Fu整体上呈现一种负相关的关系：

式（1）中，µs表示太阳天顶角SZA 的余弦值，D是地球与太阳的天文距离（与日期有关），E0是进入大气顶部的太阳辐照度（可以通过MODTRAN 内部文件获取，取1 368 W/m2）。此外，α′和β′是方程系数，可以通过公式（2）和公式（3）获得：

式（2）～（3）中，μs表示太阳天顶角SZA 的余弦值，w代表着大气可降水含量，a1～a7，x，y，z是方程待拟合系数，与天气情况有关（晴空或有云）。联合公式（1）、（2）、（3），基于上述139968 对NSSR 和Fu的样本集，结合MATLAB 和最小二乘法，可以拟合出特定天气情况下的适用系数a1～a7，x，y，z。

由于现实情况下，我们往往无法获取到大气顶层向上辐射通量（Fu）这个物理量，考虑到Fu与大气顶层宽通道反照率（r）有转换关系：

式（4）中，μs，D，E0的物理意义和公式（1）一致。因此若能获取到大气顶层宽通道反照率（r），则可以将其转换成大气顶层向上辐射通量（Fu），进而得到NSSR。但是，大部分卫星传感器只提供了窄波段的表观反射率，所以需将窄通道表观反射率向宽通道反照率进行转换，通过将传感器几个窄波段反射率进行线性组合进而获取宽波段反照率。以MODIS/AQUA 为例，宽波段反照率可以表示为：

式（5）～（7）中，i表示波段，ρi表示该波段的表观反射率，Li表示该波段的辐射亮度；μv表示卫星观测天顶角VZA 的余弦值；μs表示太阳天顶角SZA 的余弦值；D是地球与太阳的天文距离；方程（6）中的c1i，c2i，c3i，c4i是系数，需要通过最小二乘法进行拟合；Ei是该波段大气顶层的平均太阳辐照度，可以通过MODTRAN 内部文件和MODIS/AQUA 通道响应函数卷积得到，表1 列出了MODIS/AQUA 选定的6 个通道响应的平均太阳辐照度。

表1 MODIS/AQUA 选定的6 个通道大气顶层太阳辐照度Table 1 the TOA solar irradiance of the six selected MODIS/AQUA radiance bands

3 地表短波净辐射反演模型构建与验证

由上可知，构建地表短波辐射反演模型需要两个步骤：（1）窄通道表观反射率向宽通道反照率转换；（2）利用大气顶层向上辐射通量与地表短波净辐射的非线性关系反演NSSR。公式（5）～（7）是第一个步骤所需要拟合的公式，其拟合系数与太阳天顶角SZA有关，每一个VZA 对应一套系数（28 个系数）。这里共模拟了36 个SZA 值，故我们的系数查找表的维度为36*32，由于篇幅关系，这里并不作展示。

图2 展示了步骤一的拟合结果，横坐标表示MODTRAN 模拟出来的宽波段反射率，纵坐标是我们基于窄波段反射率模拟得到的宽波段反射率。可以发现，散点的样本数为13.996 8 万个，偏差Bias 为0，均方根误差RMSE 为0.011，决定系数R2为0.997，该结果保证了大气顶层窄波段反射率到大气顶层宽波段反照率转换的方法精度。图中散点几乎平均分布在1∶1 参考线附近，并且每类代表性地物具有特定的反照率范围。

图2 NSSR 反演模型构建步骤一——大气顶层宽波段反照率与估算值的散点Fig.2 Comparison of the estimated TOA broadband albedos with MODTRAN albedos

反演模型的第二个步骤，是先利用公式（4）将大气顶层宽波段反照率转换到大气顶层向上辐射通量，再通过公式（1）～（3）拟合地表短波净辐射与大气顶层向上辐射通量的非线性关系。步骤二中的方程系数与天气状况（晴空或有云）有关，即晴空或有云各自有一套方程系数，表2 是通过最小二乘法和MATLAB 拟合得到的系数结果。图3展示了NSSR 反演模型步骤二的结果，横坐标表示MODTRAN 模拟出来的地表短波净辐射，纵坐标是基于大气顶层向上辐射通量得到的地表短波净辐射。晴空和有云散点的总样本数为13.996 8 万个，偏差Bias 为0 W/m2，均方根误差RMSE 为7.29 W/m2，决定系数R2为0.999。晴空状况的NSSR 数值范围为0～1 200 W/m2，而有云情况的数值范围为0～400 W/m2，可以发现有云状况时反演结果相对比较离散。由此可见，基于MODTRAN模拟数据的地表短波净辐射反演模型的构建精度较高。

表2 NSSR 反演模型构建步骤二——方程（1-3）系数拟合统计Table 2 Fitting statistics of equation coefficients

图3 NSSR 反演模型构建二——地表短波净辐射和估算值的散点Fig.3 Comparison of the estimated NSSR with MODTRAN-modeled NSSR

基于MODTRAN 模拟数据构建出的地表短波辐射反演模型，需要经过地面站点实测数据来检验。进行地面站点验证前，为了减少异常值对数据验证结果的影响，将地面站点部分存在异常情况的数据进行了剔除。异常值的判断原则是，某瞬时测量值的前后5 min 内的标准差大于给定阈值，这样可以消除现场测量时随机噪声等带来的不良影响。以2017 年全年为研究时间，需要提取MODIS/AQUA 产品中距离相应地面站点最近的遥感影像像元值。提取的变量信息包括：MYD02 产品中的通道辐射亮度Li，MYD03 产品中的经纬度、太阳天顶角SZA、观测天顶角VZA 等观测几何参数，MYD05 产品中的大气可降水量w，MYD35 产品中的云掩膜信息。利用MODIS 云掩膜信息区分晴空和云两种情况，将Li、SZA、VZA、w等自变量输入到已构建好的NSSR 反演模型。将反演模型得到的NSSR 和地面实测NSSR 的进行比较，统计反演误差。

图4 是NSSR 反演模型的地面验证结果。图4（a）展示了站点散点图，可以看出，模型验证的总体偏差Bias 为-4.8 W/m2，均方根误差RMSE 为77.1 W/m2，决定系数R2为0.9106。经过对比探究，发现NSSR 反演模型在不同站点上的性能表现不同，其在“Bondville，IL”站点的验证结果最好，而在“Sioux Falls，SD”站点则误差较大（表3）。图4（b）则是模型误差分布直方图，可以发现超70%样本的反演误差低于50 W/m2，且误差符合正态分布。综上，所构建的地表短波净辐射反演模型在地面站点验证结果良好，具有一定泛化能力。

图4 NSSR 反演模型地面验证（a）散点图（b）误差直方图Fig.4 Validation of NSSR retrieval model using in situ measurements.

表3 地表短波净辐射在SURFRAD 站点的反演误差统计Table 3 Error statistics of NSSR retrieval at seven SURFRAD sites

4 结论

地表短波净辐射是地表辐射平衡中的一个重要组分，是许多全球、区域气候、水文、陆面过程模型中的重要参数，其高精度遥感反演对全球和区域气候变化、能量平衡、生态环境、地面模型构建、大气循环探究等领域具有相当重要的实际意义。该文以MODIS/AQUA 数据为驱动源，构建了具有一定泛化能力的高精度地表短波净辐射反演模型。

在反演模型构建过程中，利用大气传输模型MODTRAN 模拟了13.996 8 万个种复杂情况下的地表短波辐射。基于大量的模拟数据集，拟合地表短波净辐射与大气顶层观测值的非线性关系，该反演模型分为两个步骤：（1）窄通道表观反射率向宽通道反照率转换；（2）利用大气顶层向上辐射通量与地表短波净辐射的非线性关系反演NSSR。步骤一构建结果表明，其偏差Bias 为0，均方根误差RMSE 为0.011，决定系数R2为0.997。步骤二的偏差Bias 为0 W/m2，均方根误差RMSE 为7.29 W/m2，决定系数R2为0.999。误差统计结果说明地表短波净辐射反演模型的构建精度较高。

构建完地表短波净辐射反演模型后，利用地表辐射观测网（SURFRAD）7 个地面实测站点数据，结合MODIS/AQUA 遥感数据，对反演模型进行验证。验证结果表明，其总体偏差Bias 为-4.8 W/m2，均方根误差RMSE 为77.1 W/m2，决定系数R2为0.9106。研究发现，超70%样本的反演误差低于50 W/m2，在不同站点具有不同的性能表现。