引言

地球表面的辐射收支决定了地表的热状况和海洋的环流运动，是海气交换、陆面过程、海洋环流运动的能量来源，地表上行长波辐射(ULR)、地表下行长波辐射(DLR)是地表辐射收支4个参量的组成部分，长期以来是天气、气候模式的重要要素(孙帅等，2017；曹雯等，2016；张海宏等，2017；卢乃锰等，2017)，在农业、水文领域也有广泛的应用(邵小路等，2013；师春香等，2018)。关于DLR，早在20世纪30年代就开始对它计算研究，Brunt(1932)通过辐射图方式计算的DLR与地面气象观测资料做统计回归分析，建立起地面气温和水汽压计算DLR的经验函数，后来Brutsaert(1975)又基于较为精确的DLR计算和地面气温、水汽的回归，建立了DLR计算经验公式。这些公式的DLR计算都需要地面气象站的水汽和百叶箱温度作为输入，由于海洋上和高海拔地区测站的稀少，使得经验公式的应用受到局限。因此从20世纪80年代以来，科学家开始探索利用气象卫星资料反演DLR，如Darnell et al(1983)、Frouin et al(1988)借助NOAA卫星的TOVS资料和辐射传递方程计算DLR，这种方法被称为DLR的物理反演方法，但由于需要对庞大的卫星廓线资料做逐光谱的辐射率计算，耗费机时，因此国际上物理方法只用于研究，而不被业务产品处理采纳。1989年至今，陆续有统计反演方法研究出来，其反演模型总体分为2类：1类是通过建立卫星通道亮温与低层大气温度、大气柱总可降水量与低层大气比辐射率的回归关系，利用大气灰体辐射定律计算DLR，如Lee and Ellingson(1992)建立的NOAA/HIRS通道数据与DLR的统计回归模型；另1类是直接建立多个卫星相关通道辐射率与DLR的线性及非线性组合回归模式，如Wang and Liang(2009)、Tang and Li(2008)建立的从MODIS多通道辐射率计算DLR的回归模式。其中应用于如今Aqua、Terra卫星CERES仪器的DLR业务产品处理的是Gupta(1989)、Gupta et al(1992)最初提出的一套反演晴空、云天地面长波辐射通量统计模型，其天空完全云覆盖下的DLR反演模式通过建立云底温度、云底至地面的水汽含量与地面DLR的回归关系，把云底辐射与地表DLR联系起来；晴空模式则是通过近地层温度及水汽决定的低层大气比辐射率来反演DLR。相比于DLR，ULR的反演方法和模型较为简单，关于它的研究也较少，一般是通过地表温度和比辐射率用ULR物理公式进行计算，像CERES的ULR业务产品反演(Gupta et al，1997)；近年来Wang et al(2009)通过辐射传输和统计回归分析得到MODIS多通道辐射率与ULR的线性回归模式，也取得了较好的反演精度。

该文参考了CERES业务DLR反演模型，通过辐射传输正演和统计回归，建立了完全云覆盖下的地面DLR反演模式。应用模式和GFS资料以及FY-4A云检测产品反演出云覆盖地区的DLR，结合FY-4A晴空地面DLR业务产品，处理生成了2018年9月1日逐小时全天空DLR资料。反演资料经与CERES同类产品及与中国地面辐射基准站的观测相比，得到较好的精度。在下文中，重点介绍云天DLR反演模式的建立过程，云天ULR、DLR的处理方法，以及利用Aqua/CERES全天空DLR、ULR产品对FY-4A反演资料的验证试验。

1 数据资料

该文所涉及的资料有：全球大气廓线、全球预报系统(GFS)资料、FY-4A晴空DLR、ULR业务产品、Aqua卫星CERES仪器的DLR、ULR业务产品、中国地面辐射基准站DLR、ULR逐小时平均观测数据。其中178105条云大气廓线用于云天DLR反演模式建立过程中的辐射传输正演模拟，是由2010年1月、4月、7月和10月的ECMWF资料整理出22293条全球晴空大气廓线，随机地将云底从900~200 hPa每50 hPa取值，云底为900、850、800、750、700、650、600、550、500、450、400、350、300、250、200 hPa的廓线分别有4 467、3 982、7 997、7 730、9 960、9 124、11 335、10 732、13 550、12 260、13 550、12 260、13 550、12 260、17 532条；GFS资料用于云区DLR、ULR的反演；FY-4A晴空DLR、ULR业务产品用于与云区反演结果相结合、生成FY-4A观测范围内全天空的DLR、ULR资料；CERES的DLR、ULR业务产品，以及地面辐射基准站的DLR、ULR实测数据用于检验本文的反演结果。

2 云天DLR反演模式的建立 2.1 辐射传输方程及地表下行长波辐射率的模拟计算

完全云覆盖下的地表下行长波辐射率由如下公式计算：

$ {I_{v0}}(\theta) = B\left({v, {T_{{\rm{CB}}}}} \right){\tau _{v{\rm{CB}} }}(\theta) + \int_{{\tau _{v{\rm{CB}}}}}^1 B (v, T){{\rm{d}}_\tau } $

(1)

式中：I_v0是云覆盖下到达地表的长波辐射率，T_CB是云底温度，T是大气温度，v是波数，B(v, T_CB)、B(v, T)分别是温度为T_CB、T的普朗克函数，θ是当地天顶角，τ_vCB是云底到地面的光谱透过率(关于透过率的计算详见见=吴晓，1998)，τ是z高度上大气透过率。式(1)右边第一项是云底辐射的贡献，第二项是云底以下大气对地表下行辐射率的贡献。

晴空状况下的地表下行长波辐射率由如下公式计算：

$ {I_{v0}}(\theta) = \int_{{\tau _{zt}}}^1 {{B_v}} (v, T){\rm{d}}\tau $

(2)

式中：I_v0是晴空大气到达地表的长波辐射率，τ_zt是大气顶到地表的透过率，关于晴空地表下行辐射率的计算详见吴晓(2014)。

图 1和图 2分别是基于晴空和云廓线模拟的地表下行长波辐射率光谱计算结果示例。

2.2 云覆盖下地表下行长波辐射通量的模拟计算 2.2.1 云天地表总波段下行长波辐射率

基于以上辐射率的模拟，每条廓线的云天总波段地表下行长波辐射率由下式计算：

$ {I_0}(\theta) = \int_0^\infty {{I_{v0}}} (\theta){\rm{d}}v = \sum\limits_{i = 1}^{280} {{I_w}} (\theta)\Delta {v_{{i_i}}} $

(3)

式中：I₀为总波段下行长波辐射率；Δv_i=10 cm^-1；光谱区间0~∞取0~2 800 cm^-1共280个积分波段。数学上是将分段的辐射值进行积分计算，最后得到总波段的辐射率。

2.2.2 云天地表下行长波辐射通量密度的模拟计算

云天地表下行长波辐射通量密度是总波段辐射率对下半球立体角的积分，计算如下(Wark et al, 1962)：

$ DLR = \int_0^{2\pi } {\rm{d}} \varphi \int_{ - \frac{\pi }{2}}^0 {{I_0}} (\theta)\sin \theta \cos \theta {\rm{d}}\theta $

(4)

式中：φ是方位角，计算中假定I₀对方位角是各项同性的，即

$ DLR = 2\pi \int_{ - \frac{\pi }{2}}^0 {{I_0}} (\theta)\sin \theta \cos \theta {\rm{d}}\theta $

(5)

对于178 105条大气廓线，I₀有如下角度关系(Wark et al, 1962)：

$ \frac{{{I_0}(\theta)}}{{{I_0}(0)}} = 1 + a\theta + b{\theta ^2} + c{\theta ^3} $

(6)

式中：a、b、c是统计回归系数。式(6)的统计关系未免粗略，对于每一条廓线进一步定义(Wark et al, 1962)：

$ \gamma = \frac{{\left[ {\frac{{{I_0}(\theta)}}{{{I_0}(0)}} - 1} \right]}}{{a\theta + b{\theta ^2} + c{\theta ^3}}} $

(7)

对于178 105条廓线，γ有如下统计关系(Wark et al, 1962)

$ \gamma = \alpha + \beta {I_0}(0) $

(8)

将式(8)代入式(7)得到：

$ {I_0}(\theta) = {I_0}(0)\left\{ {1 + \left({a\theta + b{\theta ^2} + c{\theta ^3}} \right)\left[ {\alpha + \beta {I_0}(0)} \right]} \right\} $

(9)

将式(9)代入式(5)得到：

$ \begin{array}{l} \begin{array}{*{20}{c}} {DLR = {I_0}(0)\left[ {A + C{I_0}(0)} \right]}\\ {A = 2\pi \int_{ - \frac{\pi }{2}}^0 {\left[ {1 + \alpha \left({a\theta + b{\theta ^2} + c{\theta ^3}} \right)} \right]} \sin \theta \cos \theta {\rm{d}}\theta } \end{array}\\ C = 2\pi \int_{ - \frac{\pi }{2}}^0 \beta \left({a\theta + b{\theta ^2} + c{\theta ^3}} \right)\sin \theta \cos \theta {\rm{d}}\theta \\ \end{array} $

(10)

式中：A、C是由a、b、c、α、β计算出的常数。这样，通过式(10)即计算得到每条廓线的DLR。

2.3 云天DLR反演模型

Gupta et al(1992)将完全云覆盖下的DLR写为晴空地表长波辐射通量密度与云辐射强迫通量密度之和，下式为：

$ DLR = DL{R_{{\rm{clear }}}} + DL{R_{{\rm{cloud }}}} = {C_1} + {C_2} $

(11)

式中：DLR、DLR_clear(C₁)分别是云天和晴空的地表下行长波辐射通量密度(晴空DLR_clear的模拟计算详见吴晓(2014))，DLR_cloud(C₂)是云辐射强迫通量密度。Gupta et al(1992)通过敏感度试验，发现C₂与云底黑体辐射和云底到地面的水汽总量分别存在高度正和高度负相关，经过敏感度试验，提出如下模式：

$ {C_2} = \frac{{\sigma T_{{\rm{CB}}}^4}}{{{a_0} + {a_1}w + {a_2}{w^2} + {a_3}{w^3}}} $

(12)

式中：σ是斯蒂芬-玻耳兹曼常数；w是云底至地面的水汽总可降水量(单位：cm)；T_CB是云底温度；a₀、a₁、a₂、a₃是统计回归系数。

这里将式(12)改写为：

$ {C_2} = \frac{{\sigma T_{{\rm{CB}}}^4}}{\varepsilon } $

(13)

进一步整理有下式：

$ \varepsilon = \frac{{\sigma T_{{\rm{CB}}}^4}}{{{C_2}}} = \frac{{\sigma T_{{\rm{CB}}}^4}}{{DLR - DL{R_{{\rm{clear}}}}}} $

(14)

因此对于每一条廓线可由模拟的DLR与DLR_clear及云底温度计算得到一个云辐射强迫系数ε，178 105条廓线共178 105个ε，通过与廓线云底至地面的水汽总含量w统计回归。图 3是地面气压与云底气压相差200~300 hPa廓线的ε与w回归图，由图得到ε与w的回归关系式及系数：

$ \varepsilon = {a_0} + {a_1}w + {a_2}{w^2} + {a_3}{w^3} $

(15)

3 云天DLR和ULR反演方法

将以上模式结合GFS全球预报场资料，即可估算云覆盖地区的DLR。

3.1 云底至地面的大气可降水量计算

GFS资料提供了覆盖全球、数据分辨率为0.5° ×0.5°的大气温度、相对湿度、地面温度、地面气压等多种参数，其中大气垂直温湿度为31层(1000~0.1 hPa)，通过对相对湿度>90%判断为云区，得到对应的云底气压和云底温度，由下式计算云底到地面的大气总可降水量w_CB

$ {w_{{\rm{CB}}}} = \int_0^{{z_{{\rm{CB}}}}} {{w_i}} {\rm{d}}z $

(16)

式中：w_i是z高度上的绝对湿度，z_CB是云底高度，w_i由下式计算：

$ {w_i} = 0.217\frac{e}{T} $

(17)

式中：e是z高度上的水汽压(单位：hPa)；T是大气温度(单位：K)。

3.2 计算云辐射强迫通量密度

云辐射强迫通量密度C₂由下式计算：

$ {C_2} = \frac{{\sigma T_{{\rm{CB}}}^4}}{{{a_0} + {a_1}w + {a_2}{w^2} + {a_3}{w^3}}} $

(18)

式中回归系数见表 1。

3.3 云天DLR计算方法

云存在下的地表DLR由下式计算：

$ \begin{array}{*{20}{c}} {DFLUX = (1 - A) \times DFLU{X_{{\rm{clear }}}} + }\\ {A \times \left({DFLU{X_{{\rm{clear }}}} + {C_2}} \right)} \end{array} $

(19)

式中：右边第一项是晴空辐射贡献，第二项是云辐射贡献；A是云量，DFLUX是云天DLR，DFLUX_clear是晴空DLR，C₂是云辐射强迫通量密度。DFLUX_clear的反演模式如下：

$ {T_{\rm{e}}} = {k_1}{T_1} + {k_2}{T_2} + {k_3}{T_s} $

(20)

式中：T_e是低层大气有效辐射温度，T₁是与地面气压差150 hPa以内的低层大气平均温度，T₂是与地面气压差300~150 hPa的低层大气平均温度，T_s是地面温度；k₁=0.35，k₂=0.05, k₃=0.6。

$ DFLU{X_{{\rm{clear}}}} = {\varepsilon _{\rm{s}}}\sigma T_e^4 $

(21)

$ {\varepsilon _{\rm{s}}} = {b_0} + {b_1}{w_{\rm{c}}} + {b_2}\sqrt {{w_{\rm{c}}}} $

(22)

式中：b₀、b₁、b₂是统计回归系数，w_c是大气柱总可降水量，ε_s是低层大气比辐射率；关于DFLUX_clear的反演参见吴晓(2014)、吴晓等(2018)。这里需要说明的是，由于FY-4A云产品没有提供观测像元的云量，因此处理时将卫星确认的云像元云量做为100%，即A=1.0；同时将相对湿度>85%的GFS大气层湿度做调整，以确保利用此时的云区GFS廓线计算出晴空DFLUX_clear。

3.4 云天ULR计算方法

与DLR计算方法不一样，这里主要利用GFS资料计算云天ULR：

$ ULR = {\varepsilon _{{\rm{SF}}}}\sigma T_{\rm{s}}^4 + \left({1 - {\varepsilon _{{\rm{SF}}}}} \right)DLR $

(23)

式中：T_s是GFS地表温度，ε_SF是全球逐月的地表宽波段比辐射率气候数据；式(23)右边第一项是地表的向上发射辐射，第二项是地表反射的大气向下辐射。

4 FY-4A全天空DLR、ULR反演结果及其对比验证 4.1 FY-4A全天空DLR、ULR反演结果示例

本文中晴空地区的DLR、ULR估算直接采用国家卫星气象中心FY-4A业务系统输出的晴空DLR、ULR产品，产品通过FY-4A成像仪多个红外通道辐射率遥测数据经反演模式计算得到(吴晓，2014；吴晓等，2018)。结合以上云天ULR、DLR反演结果和FY-4A晴空ULR、DLR业务产品，处理得到2018年9月1日逐小时FY-4A观测范围内全天空状况下DLR、ULR资料。图 4是2018年9月1日18时(UTC)全圆盘ULR、DLR。

4.2 对比验证

本文主要采用美国Aqua卫星CERES仪器的全球瞬时DLR、ULR反演产品，以及中国地面辐射基准站的逐小时DLR、ULR观测资料，对反演的FY-4A全天空DLR、ULR资料进行对比验证。

4.2.1 FY-4A瞬时DLR资料与CERES的对比误差分析

Aqua是极轨卫星，过境时间为当地时02、14时左右，每天对地球上同一地点位置进行白天、夜间2次观测。它在每天的04—09 UTC, 15—23 UTC进入FY-4A卫星的观测范围，其星下点数据分辨率为28 km，约5~6条轨道覆盖完整FY-4A观测范围，取与每条轨道过境时间相差30 min以内的FY-4A数据与CERES轨道数据匹配，并将逐条轨道的数据拼接，形成50°S~50°N、55°~155°E、数据分辨率为0.5°×0.5°的两星匹配DLR、ULR格点场。图 5是匹配的2018年9月1日夜间FY-4A瞬时DLR、Aqua/CERES瞬时DLR。

CERES仪器的DLR产品，是一个逐轨道、逐像点瞬时视场反演产品，是在NASA再分析场资料MOA(meteorology ozone and aerosol)基础上，利用Gupta云天、晴空DLR反演模式和MOA的大气温湿度廓线，以及Aqua卫星云量产品计算得到的。图 6是匹配的2018年9月1日白天FY-4A与CERES瞬时DLR格点场，对白天、夜间两星DLR场做误差统计，结果为：RMSE=20.52 W·m^-2，R=0.9481，BIAS=3.3 W·m^-2(夜间格点场)；RMSE=25.58 W·m^-2，R=0.9096，BIAS=5.4 W·m^-2 (白天)。匹配格点场的差值如图 7，由图 7可见50°S~50°N、55°~155°E范围内大部分地区，FY-4A DLR与CERES的偏差在30 W·m^-2以内，对比误差大于30 W·m^-2均在大陆晴空地区，这是因为在晴空地区FY-4A与CERES的DLR反演方式不同，FY-4A晴空地区DLR反演利用的是卫星通道遥测数据，而CERES的DLR反演则是利用NASA再分析资料MOA的温湿度数据，目前再分析资料对陆地表面温度的模拟偏低实况(特别是在白天)，在某些地区达6 K以上，这使得利用MOA温度经反演模式估算的DLR偏低于实况，反映在大陆晴空地区FY-4A的DLR高于CERES，差值图上呈大的正偏差。图 8是夜间匹配数据的散点图。误差的另一个原因则是两星云区DLR反演采用的再分析资料不同(FY-4A采用GFS、CERES采用MOA)，由于再分资料与实况、以及不同再分析资料(如GFS、MOA)之间存在一些差异，如温湿度廓线、对云模拟的不尽相同，也造成FY-4A DLR与CERES对比的误差。当然，DLR反演模式的不同也是误差的来源。

4.2.2 FY-4A瞬时ULR资料与CERES的对比及误差分析

对2018年9月1日的FY-4A实时ULR与CERES产品做时空匹配，得到范围在50°S~50°N、55°~155°E、分辨率为0.5°×0.5°的两星ULR格点场。图 9和图 10分别是夜间、白天的两星瞬时匹配格点场。图 9中, RMSE=10.97 W·m^-2，R=0.9762，BIAS=-3.3 W·m^-2；图 10中, RMSE=19.97 W·m^-2，R=0.9283，BIAS=5.0 W·m^-2。图 11是夜间两者对比散点图。误差来源首先是再分析资料的陆地表面温度在白天明显偏低于实况多达6 K以上，而在夜间相对正确，CERES的ULR估算完全依赖于再分析资料的地表温度，而FY-4A的晴空ULR反演采用的是实时卫星通道亮温，这反映在两种ULR资料对比上有误差夜间相对小于白天。

4.2.3 FY-4A DLR、ULR与河南中国地面辐射基准站观测数据的比较

目前我国有3个地面辐射基准站，使用较为精准的辐射仪连续观测地面和大气的辐射。这里用北京时2018年9月1日08时到2日07时河南站(34.073°N、113.905°E)的实时观测与反演资料作对比(图 12、图 13)，图中DLR最大偏差约为40 W·m^-2、ULR最大偏差为20 W·m^-2；由图 13可见在9月2日07时反演DLR与地面实测相差最大、约为40 W·m^-2，这可能是卫星视场内有云，但云检测为晴空，因此造成反演DLR偏低。

5 结论

本文通过红外辐射传输方程以及全球大气廓线，模拟计算了有云存在下逐条廓线的地表下行长波辐射通量密度(DLR)，基于DLR与廓线云底黑体辐射通量以及云底至地面水汽总量的统计回归，建立了云辐射强迫通量密度与云底辐射的统计回归关系；通过借鉴Aqua、Terra卫星CERES仪器的DLR业务产品算法模型，建立了云条件下地表下行长波辐射通量密度的反演模式。应用模式和GFS资料处理了精度合理的云区DLR资料，同时利用GFS的地表温度也估算得到精度较好的云区ULR资料；而晴空地区的DLR、ULR沿用了FY-4A的业务反演产品。

全天空地表长波辐射通量DLR、ULR的反演，最困难、最不确定的是云天DLR反演。相对于近年来多种晴空反演模式的提出，云天DLR的反演方法为数不多，大多都是利用数值和预报场的温湿廓线来反演DLR，如CERES的DLR业务产品处理，由于云的因素，如再分析资料对云模拟的不确定性，使得反演云天DLR与实况存在一些误差；同样，由于云的影响，利用卫星通道遥测数据无法反演得到云天ULR，只能借助预报场资料，可见DLR、ULR反演精度的提高都有赖于精度更高的再分析场资料；另外进一步研究发展反演误差小的DLR、ULR反演模式，也是提高地表长波辐射通量反演的有效途径。