1 引言

大气辐射传输过程在气候模拟、大气遥感反演及大气光化学过程中起着非常重要的作用。而辐射方案的精度差异会对大气数值模式模拟效果产生一定影响(杨学胜等，2009；许尤福等，2011；荆现文等，2012；张梦等，2014；Zhang H et al，2014)。辐射传输的物理过程通常由辐射传输方程来描述，而该方程属于积分-微分方程，很难获得解析解，因此，已经有大量研究发展了一些近似方案来求解大气辐射传输方程，如二流近似方案(King et al，1986；Kylling et al，1995；Fu et al，1997；邱金桓，1998；Qiu，1999；Li et al，2005；Lu et al，2009)和四流近似方案(Liou et al，1988；Stamnes et al，1988；Li et al，1996；Zhang F et al，2013a, 2013b, 2016；Zhang H et al，2015)。二流近似是指用两个自由度来表示辐射强度在球空间的分布特征。多年来，由于其形式简单、计算速度快，被广泛应用于参数化气候模式的辐射传输过程。但这种方法有明显的缺陷和不足。King等(1986)对各种单层的二流近似解的精度进行了评估，指出在不同的光学厚度、太阳天顶角和单次散射反照率的情况下，采用二流近似方法计算得到的大气透射率和反射率的相对误差高达15%—20%。Lu等(2009)指出二流近似可以导致10%以上的云顶加热率误差，而云顶加热率影响着云的发生、发展。随着气候模式时空分辨率的提高和物理过程的改进，其对一个高精度、高速度的辐射算法的需求显得越来越迫切。张峰(2013)及Zhang F等(2013a, 2013b, 2016)基于辐射传输的4个不变性原理，建立了四流累加辐射传输方案。该方案越来越受到重视，并能显著提高对气溶胶辐射强迫的评估精度，因此在气候模式中得到广泛应用(Zhang H et al，2015；Li et al，2015)，在辐射传输计算过程中，各种近似方案都需要对散射相函数进行相应阶数的Legendre函数展开。例如，二流近似方案需要对散射相函数进行二阶Legendre函数展开，四流近似方案需要对散射相函数进行四阶Legendre函数展开。目前，无论是单柱的辐射模式，还是气候模式中的辐射方案，其四流近似方案通常采用Henyey-Greenstein(HG)方案(Henyey et al，1941)对散射相函数进行Legendre函数展开(Zhang F et al，2013a；Zhang H et al，2015)。Zhang F等(2013a)对新建的四流球谐函数展开累加方案的精度进行了较为全面的评估，认为该方案能显著提高辐射通量和加热率的计算精度，但以上结论同样是基于HG方案得到的。而在HG方案中，高阶展开采用的仅仅是非对称因子g的幂次方。这种单参数近似不能准确地表示相函数的特征，尤其是无法表述相函数的后向散射特性(Kattawar，1975)，因此，HG方案会直接影响到辐射方案本身的计算精度(Li et al，2015)。为了改进这种近似方法，Kattawar(1975)提出了双Henyey-Greenstein(DHG)方案，是指用两个HG相函数的线性组合来表示散射相函数的变化特征，即用3参数来表示相函数。相比单参数近似的HG方案，该方案具有更大的自由度去刻画散射相函数各向异性的特性。DHG方案提出已经有40多年，但该方案在辐射传输计算中应用甚少。此外，相比HG方案，该方案在辐射计算中能多大程度的提高计算精度，也均未见有相关研究。本研究对DHG方案进行了较为系统的研究，提出了改进的DHG方案(MDHG)，并评估了HG、DHG和MDHG这3种方案的计算精度。

2 相函数基本理论

散射相函数是散射角的函数。对于球形粒子，可以利用米散射理论计算出其散射相函数；对于非球形粒子，可以利用T矩阵(Mishchenko et al，1994)、几何光学(Yang et al，2007)和Bridging(Zhao et al，2006)技术等方法计算出其散射相函数。为了便于辐射传输计算，文中将相函数用已知的数学函数表达，采用勒让德多项式展开的相函数可以写成(Henyey et al，1941)

(1)

式中，，P_n(cosΘ)为勒让德多项式的相函数展开项，散射角的余弦cosΘ为

(2)

式中，μ为光的入射方向，μ′为光的散射方向。对式(1)进行方位角平均可得

(3)

式中，P_n(cosΘ)可以展开为

(4)

根据利用式(3)—(4)，可得

(5)

2.1 Henyey-Greenstein近似

目前，应用最为广泛的相函数是Henyey等(1941)提出的Henyey-Greenstein(HG)相函数。该方案采用非对称因子g表示相函数

(6)

式(6)可以展成为

(7)

根据式(3)、(5)和(7)可得方位角平均的HG近似相函数

(8)

式中，ω_n=(2n+1)gⁿ，并将这种方案记为HG。在该方案中，相函数的高阶展开ω_n仅用非对称因子g表示，因此，在辐射传输计算中得到广泛应用(Lu et al，2009；Liou et al，1988；Stamnes et al，1988；Li et al，1996；Zhang F et al，2013a；Zhang H et al，2015)。

2.2 双Henyey-Greenstein近似

双Henyey-Greenstein(DHG)(Kattawar，1975)近似是指由两个HG相函数的线性组合来表示相函数的散射特性。相比单参数的HG方案而言，DHG方案有3个待定参数，因此，有更多的自由度来表现相函数的真实特征。DHG方案采用两个HG方案的权重组合来计算相函数

(9)

根据上式，方位角平均的DHG近似相函数可表示为

(10)

根据式(5)，可得到式(10)中3个参数a、g₁和g₂满足的方程

(11)

(12)

(13)

式中，。

根据以上约束方程，Kattawar(1975)给出了a、g₁和g₂的解

(14)

(15)

(16)

为方便表述，将Kattawar(1975)提出的以上方案记为DHG。该方案的本质是采用相函数的展开项ω₁、ω₂、ω₃来参数化相函数。在一些特定情况下，DHG方案得到的相函数存在异常峰值，而这并不符合相函数的真实特征(Kattawar，1975)。因此，该方案在辐射传输计算中几乎没有得到应用。为了解决这一问题，本研究采用ω₁，ω₂两项的相函数展开系数和相函数的整体特征来求解DHG方案的3个参数。具体步骤如下：

(1) 根据式(11)和(12)可得到包含有未知数a的方程的解

(17)

(18)

(2) 权重a从0开始以0.001为间隔循环至2，利用式(17)、(18)能得到2000种不同的a、g₁和g₂组合；

(3) 从以上组合中挑选出满足|g₁|＞|g₂|和|g₁|＜1的a、g₁和g₂组合。在|g₁|和|g₂|均小于1时，使得相函数的高阶展开项中包含的g₁ⁿ和g₂ⁿ在n趋向于无穷后仍收敛。

(4) 以上每一种组合均能利用式(9)构造相函数，把构造的相函数和真实的相函数(根据粒子的谱分布和已有的粒子散射理论，就能计算出真实的相函数。文中计算得到真实的相函数从0°开始以0.3°为间隔，共包含601个值)进行比较，求得相对误差最小的一组a、g₁和g₂组合，这种方案记为MDHG。

3 计算结果比较

大气中气溶胶和云的强散射过程主要发生在短波波段。而在短波波段中，可见光波段的能量最强，通常采用波长0.55 μm来表示可见光波段的平均情况。因此，文中计算波长0.55 μm处的气溶胶和水云粒子的散射相函数，并分别利用HG、DHG和MDHG方案对其进行参数化，进一步讨论这几种参数化下相函数的精度。此外，利用四流球谐函数展开累加方案(Zhang F et al，2013a)比较分析以上3种不同相函数参数化方案下计算的反射率、透射率或吸收率的精度，将128流离散纵标的计算结果作为精确值(Stamnes et al，1988)。四流球谐函数展开累加方案的主要思想是用4个球谐函数来表示辐射强度在角度空间的分布特征，采用ω₄来进行δ调整(Wiscombe，1977)。该方案能显著改善辐射通量和加热率的计算，并已在北京气候中心气候系统模式(BCC_AGCM2.0.1)和加拿大气候模拟与分析中心的气候模式(CanAM4)中得到应用(Zhang H et al，2015；Li et al，2015)。

3.1 沙尘气溶胶

对于沙尘气溶胶，文中考虑随机取向的椭球粒子，并采用对数正态分布来描述沙尘粒子的分布。该分布函数曾为世界气象组织推荐(Sokolik et al，1998)

(19)

式中，r为椭球沙尘等效为等体积球体的半径，r_g为粒子群的模态半径，σ_g为标准偏差，N为气溶胶粒子的数浓度(cm^-3)。其有效半径为r_e=，有效方差为v_e=exp(ln²σ_g)-1。此处选用核模态和粗模态这两种模态，其粒子群模态半径r_g分别为0.07和1.90 μm，标准偏差σ_g分别为1.95和2.15。椭球粒子的纵横比为1.7(Fu et al，2009)。沙尘气溶胶的复折射指数为m=1.53+0.008i。利用T矩阵(Mishchenko et al，1994)和几何光学近似(Yang et al，2007)，计算得到了这两种模态在波长0.55 μm处的单次散射反照率和不对称因子：粗模态沙尘为ω=0.63和g=0.81，核模态沙尘为ω=0.95和g=0.68。表 1分别给出了核模态和粗模态沙尘在DHG和MDHG方案中的a、g₁、g₂和ω₄的值。

图 1a、b分别给出了核模态和粗模态沙尘气溶胶相函数的精确值(Exact)以及HG、DHG和MDHG方案的结果。可以看出，粗模态沙尘具有很强的前向峰值，其后向也有一个小的峰值。HG方案的相函数较平滑，但其在前向和后向并不接近真值，尤其没有表现出前向的强峰值。DHG方案的相函数在散射角为54°—170°时比其他两种方案更接近真值，但其在170°后出现了相函数为负数的情况，因此该方案没有体现出粗模态沙尘的后向峰值特征。本研究提出的MDHG方案的相函数在散射最强的前向几乎与精确值重合，对前向强峰值的模拟精度较高。此外，MDHG方案在后向也较HG和DHG方案更接近于真值。和粗模态沙尘相比，核模态沙尘没有很强的前向峰值。HG、DHG和MDHG这3种方案的相函数在前向均与真值较为接近，精度相近。而在后向，这3种方案的相函数有较大差异。DHG方案的相函数在后向出现了骤降，与真值差异最大，精度最差。HG方案的相函数虽然稳定，在后向的精度高于DHG方案，但与真值相比仍有较大的误差。而MDHG方案的相函数表现出了相函数的后向特征，更接近真实的相函数。

此外，还比较了不同相函数参数化方案下的四流球谐函数展开累加方案的反射率和吸收率的精度。核模态沙尘气溶胶的反射率和吸收率在不同方案下的精度如图 2a—f所示。从反射率的误差结果可以看出，DHG方案相比其他两种方案误差更大，尤其在光学厚度小于0.2且μ₀为0.5—0.7和0.9—1.0时相对误差可达5%。而在此范围内，HG和MDHG方案的相对误差约为2%。当光学厚度大于1时，MDHG方案中相对误差小于1%的范围比HG和DHG方案大。从吸收率的误差结果来看，3种方案精度相当，MDHG方案在μ₀为0.4—0.5且光学厚度大于3时，精度比HG和DHG方案略高。故总体而言，MDHG方案比HG和DHG方案精度更高。

粗模态沙尘气溶胶在不同方案下的反射率和吸收率的精度如图 2g—l所示。从反射率的误差结果来看，HG方案精度最低，在μ₀为0.1—0.4和0.8—1.0时误差都大于10%。DHG方案误差最小，仅当在μ₀小于0.1时相对误差大于10%。MDHG方案在μ₀小于0.2时以及光学厚度小于2且μ₀为0.8—1.0时的相对误差皆大于10%，但其相对误差大于10%的范围小于HG方案，故MDHG方案的精度高于HG方案。从吸收率的误差结果可以看出，这3种方案的精度都较高。HG方案在μ₀为0.05—0.2且光学厚度在0.1—1时误差大于5%，而其他两种方案均小于5%。此外，HG方案在光学厚度大于2且μ₀为0.25—0.45时相对误差大于1%，而其他两种方案均接近于0。相比之下，3种方案中HG方案精度最低。DHG方案在μ₀小于0.1时，误差小于MDHG方案。故整体而言，DHG方案相对误差最小，精度略优于MDHG方案，HG方案误差最大。

3.2 水云

水云的粒子尺度分布接近于Γ分布(Pruppacher et al，1996)

(20)

式中，A、α和β为常数，r为水云粒子的半径，有效半径为，有效方差为v_e=。

在计算云的光学特性的参数化试验中，为了方便计算，选用水云的有效方差为v_e=0.172(Fu et al，1997)，有效半径为r_e=5.89 μm(Fu et al，1997)。该水云可以代表典型的低云特征。水云的复折射指数为m=1.33+1.79×10^-9i。计算得到波长0.55 μm处水云的ω和g分别为1和0.85。表 2中给出了在水云情况下，DHG和MDHG方案得到的a、g₁、g₂和ω₄的值。

由图 1c可以看出，水云的HG方案的相函数没有很好地反映出相函数的真实特征，特别是相函数的后向特征。DHG方案在0°—132°的散射角下精度较高，但在后向出现了异常增大，甚至超过了真实相函数的前向最大值。在3种方案中，MDHG方案的相函数是整体上最接近于真值的，尤其是其后向精度高于其他两种方案，能很好地反应出真实相函数的后向特征。

在波长0.55 μm处，水云的吸收作用非常小，因此只讨论其反射率和透射率的情况。图 3a—f为3种方案在水云情况下的反射率和透射率的精度比较。MDHG方案在μ₀为0.4—1.0时相对误差小于2%，优于其他两种方案。这3种方案的透射率结果基本相同，在μ₀为0.4—1.0时误差几乎为0。故总体来看，MDHG方案相对误差最小，DHG方案其次，这两种方案整体都要好于HG方案。

3.3 霾

对于霾粒子，采用了Haze L的相函数散射模型，其非对称因子g=0.80(Garcia et al，1985)。从图 1d可以看出，HG、DHG和MDHG方案的相函数都几乎重合，均较为接近真实相函数。然而如表 3所示，将MDHG方案得到的最优a、g₁、g₂和ω₄与DHG方案得到的值进行比较，发现DHG方案中，g₂远大于g₁。此外，DHG方案中ω₄出现了较大负值，这导致了光学厚度在进行δ调整后反而增大，从而影响之后的反射率和透射率的计算精度，增大反射率和透射率的相对误差。而MDHG方案得到的g₁、g₂和ω₄均非常合理。

在计算HG、DHG和MDHG方案的反射率和透射率时，令单次散射反照率为1，从3种方案的反射率和透射率的精度(图 4a—f)可以看出，DHG方案下的反射率在光学厚度小于1时误差大于20%，精度明显低于其他两种方案；同样其透射率的相对误差也是3种方案中最大的。HG和MDHG方案的精度相当，这两种方案的透射率在μ₀为0.5—1.0时的相对误差几乎为0。

4 结论

本研究对DHG方案进行改进，提出了MDHG方案，并将HG、DHG和MDHG这3种方案的相函数进行了系统的对比。由相函数的比较结果发现：HG方案的结果虽然稳定但不能表现出相函数的前向或后向的峰值特征，例如该方案无法准确地描述粗模态沙尘气溶胶相函数的前向峰值和后向峰值以及水云相函数的后向峰值。DHG方案的相函数在后向易出现异常峰值或负值，例如该方案描述的水云相函数出现异常峰值以及粗模态沙尘气溶胶相函数出现负值。整体而言，MDHG方案精度最高，能表现出相函数的前、后向峰值的特征。此外，还利用四流球谐函数展开累加方案，进一步比较了不同相函数参数化方案下的反射率、透射率和吸收率的计算精度。由反射率、透射率和吸收率的计算结果发现：MDHG方案整体更为精确，而HG方案精度较低，例如在水云情况下其误差大于DHG和MDHG方案。DHG方案在沙尘气溶胶和水云情况下计算的结果都比较精确，但是对霾的计算误差较大。

综合比较相函数以及反射率、透射率和吸收率的计算精度，MDHG方案得到的相函数最能表现其前向和后向特征，且该方案计算的反射率、透射率和吸收率的精度也总体优于HG和DHG方案。因此，MDHG方案是理想的相函数参数化方案。

附录

在无时间变化，弹性、独立散射的近似条件下，方位角平均的平面平行的辐射传输方程为

(A1)

(A2)

(A3)

式中，I(τ, μ)是指介质中光学厚度为τ，天顶角的方向余弦为μ的辐射强度。F₀为大气顶垂直于太阳光方向的辐射通量，τ₀是该介质的总光学厚度，其中方位角平均的相函数为

(A4)

球谐函数谱展开的目的就是将辐射强度对角度的依赖进行分离，因此，假设散射强度为

(A5)

取N=3，将辐射强度的在角度空间分布用4个自由度来表示，即为四流球谐函数展开方法，根据辐射传输方程可得

(A6)

(A7)

(A8)

(A9)

上式中，。以上方程可以用矩阵形式表示

(A10)

利用一阶线性微分方程组基本理论，求得以上方程的解

(A11)

式中，参数f₀=1/μ₀，G=[C₁ D₁ C₂ D₂]^T，e₁=e^-k₁τ，e₂=e^-k₂τ，e₃=e^{-k₁(τ₀-τ)}，e₄=e^{-k₂(τ₀-τ)}，R_{1, 2}=-a₀/k_{1, 2}，，其他相关参数k_{1, 2}和η_{0, 1, 2, 3}相关定义见Zhang F等(2013a)和和张峰(2013)。G=[C₁ D₁ C₂ D₂]^T可由边界条件式(A2)—(A3)得到。根据式(A5)并取N=3，可得向上、向下的辐射通量为

(A12)

(A13)

因此，该层介质的反射率r(τ₀, μ₀)、透射率t(τ₀, μ₀)、吸收率a(τ₀, μ₀)为

(A14)

(A15)

(A16)

对于计算多层非均匀介质联合辐射通量，可采用累加法进行计算。首先，从第1层开始不断向下迭代，得到第1层至第k层联合的透射矩阵T_{1, k}(μ₀)和反射矩阵R_{1, k}^*

(A17)

(A18)

式中，为第1层至N层累积的光学厚度，为2×2的单位矩阵。其次，从地表(第N层)开始不断向上迭代，得到第N至第k层联合的反射矩阵R_{k, N}(μ₀)和R_{k, N}

(A19)

(A20)

式中，，r_s为朗伯地表反照率。在k+1层次处(第k层的下边界)向上和向下的无量纲的辐射强度和D_k+1(μ₀)=为

(A21)

(A22)

因此，在k+1层(第k层的下边界)处向上和向下的辐射通量分别

(A23)

(A24)

大气层顶向上和向下的辐射通量为

(A25)

(A26)

式中，u_k+1、d_k+1和r_{1, N}为矩阵U_k+1(μ₀)、D_k+1(μ₀)和R_{1, N}(μ₀)的第1个元素。关于T_{1, k}(μ₀)、R_{1, k}^*、R_{k, N}(μ₀)、R、U_k+1(μ₀)和D_k+1(μ₀)的相关定义见Zhang F等(2013a)和和张峰(2013)。