一种基于辐射通量密度的三维云层辐射计算方法与流程

本发明涉及大气辐射计算领域，尤其是一种基于辐射通量密度的三维云层辐射计算方法。

背景技术：

云是自然界中十分复杂的物理现象之一，在大气辐射传输过程中产生着较大的影响。由于云在三维空间内的分布的不均匀性，使得三维云层场景的辐射传输问题成为大气辐射传输领域中的难题。为了研究有云大气中的辐射传输问题，需要充分考虑云的三维结构和辐射特性。

三维大气辐射传输的精确计算，相对应的是成本的大幅度提高。其中，三维空间内多次散射的精确计算很大程度上影响着最终计算成本。通常可以通过多次迭代方法求解多次散射源函数，复杂的多次迭代过程势必会导致较高的成本和较慢的运算速度。因此，合理地简化多次散射源函数的计算过程来降低所需成本并提高运算速度，成为大气辐射传输领域亟需研究讨论的问题。

技术实现要素：

为了解决上述背景技术中的问题，本发明旨在提供一种基于辐射通量密度的三维云层辐射计算方法，能够在保持较高计算精度的前提下，加快计算三维云层辐射源函数，降低计算成本，并且可以很好地应用于三维多云场景。

为解决上述技术问题，本发明提供一种基于辐射通量密度的三维云层辐射计算方法，包括如下步骤：

(1)利用高斯方向和大气属性构建水平辐射通量密度传输模型；

(2)利用爱丁顿近似法和严格单次散射计算式计算三维云层基础辐射场；

(3)考虑邻云辐射效应对水平传输方向上的三维云层基础辐射场的影响，将其分为邻云影响辐射增强和邻云影响辐射减弱，计算由邻云辐射效应影响所产生的水平辐射通量密度变化；

(4)将水平辐射通量密度变化映射回高斯方向，再将三维云层基础辐射场与变化辐射强度结合得到更新后三维云层辐射强度，并由其计算多次散射源函数；

(5)结合由初始条件计算的单次散射源函数与热源函数，得到三维云层辐射源函数。

优选的，步骤(1)中，利用高斯方向和大气属性构建水平辐射通量密度传输模型，具体步骤如下：

(1)将三维云层场景分为所需大小尺度的网格单元，每个单元都有其相邻的四个单元(如图2)，利用大气属性文件，将大气属性根据空间位置存储在每个单元内；

(2)基于离散坐标下的高斯积分，将三维云层辐射场离散到已设置个数的高斯方向上，且同一方位角的高斯方向总和构成一个“三棱柱”空间(如图3)；

(3)将高斯方向上的三维云层辐射场映射到水平辐射传输方向上。

优选的，步骤(3)中，考虑邻云辐射效应对水平传输方向上的三维云层基础辐射场的影响，将其分为邻云影响辐射增强和邻云影响辐射减弱，具体步骤如下：

(1)基于高斯方向上的由爱丁顿近似法和严格单次散射计算式计算得到的三维云层基础辐射场，将高斯方向上的辐射强度映射到水平辐射通量密度传输方向上；

(2)将邻云辐射效应分为邻云影响辐射增强和邻云影响辐射减弱即i⁺(τ；ω′ts)和i^-(τ；ω′ts)，得到水平辐射通量密度的增强值或减弱值；

(3)考虑太阳天顶角变化导致的辐射误差，对水平辐射通量密度的增强值和减弱值进行修正。

优选的，步骤(3)中，计算由邻云辐射效应影响所产生的水平辐射通量密度变化，具体步骤如下：

(1)选择三维云层场景内第一个单元a，再搜索与其水平方向相邻的四个单元b1、b2、b3和b4；

(2)将基础辐射场映射到水平传输方向，计算映射系数ψts；

ψts＝cos(δφts)sin(θts)

其中θts高斯方向的天顶角；δφts高斯方向在xy平面上的投影与其所处ts空间内水平坐标轴方向的平面夹角，δφts的取值范围是-45°＜δφts＜45°。

(3)计算水平辐射通量密度在单元a与其水平相邻单元b之间的反射比例，即反照率r；

其中，τ为光学厚度，ω为单次散射反照率，g为不对称因子。

(4)计算单元a的水平辐射通量密度受邻云辐射影响的变化i⁺(τ；ω′ts)和i^-(τ；ω′ts)；

其中i⁺(τ；ω′t1)是ω′t1方向水平辐射通量密度的增加值；i^-(τ；ω′t1)是ω′t1方向水平辐射通量密度的减弱值；ω′t2方向是与ω′t1相同高斯方向相反方位的方向(如图3)；r1是水平辐射通量密度在a与b1的交界处发生反射的比例；r2是水平辐射通量密度在a与b2的交界处发生反射的比例。

(5)计算单元a由于邻云辐射效应影响所产生的辐射射通量密度变化ice(τ；ω′ts)，计算修正系数ace；

ace＝exp(sinθ0-sinθe)-1

ice(τ；ω′ts)＝(1+ace)i⁺(τ；ω′ts)-(1-ace)i^-(τ；ω′ts)

其中θ0为太阳天顶角；θe为二流近似法计算辐射通量所使用的高斯方向天顶角；ice(τ；ω′ts)代表某一单元处ω′ts方向由邻云辐射效应所产生的水平辐射通量密度变化值，i⁺(τ；ω′ts)为ω′ts方向水平辐射通量密度的增加值，i^-(τ；ω′ts)为ω′ts方向水平辐射通量密度的减弱值，ace为弥补太阳方向变化对模型带来误差的修正系数。

(6)重复步骤1-5，直到选则到三维云层场景内最后一个单元a，并且单元a的ice(τ；ω′ts)都计算完成为止。

优选的，步骤(4)中，将水平辐射通量密度变化映射回高斯方向，再将三维云层基础辐射场与变化辐射强度结合得到更新后三维云层辐射强度，并由其计算多次散射源函数，具体步骤如下：

(1)在三维云层场景内选取第一个单元a；

(2)计算从水平通量密度方向到高斯方向的分配权重p(ω′ts；ωs)和高斯积分权重ats(ω′ts)，并由此计算高斯方向变化辐射强度i'ce(τ；ω′)：

ice'(τ；ω′)＝p(ω′ts；ωs)ats(ω′ts)ice(τ；ω′)

(3)利用爱丁顿近似法计算和严格单次散射计算式计算得到单元a的三维云层基础辐射场，结合变化辐射强度得到更新后三维云层辐射强度：

其中iems(τ；ω′)为由爱丁顿近似计算得到三维云层近似多次散射强度，iss(τ；ω′)为由严格单次散射计算式计算得到的三维云层单次散射强度，i”(τ；ω′)为以三维云层基础辐射场i′(τ；ω′)作为输入计算得到的更新后三维云层辐射强度；

(4)利用更新后三维云层辐射亮度计算单元a处的多次散射源函数；

(5)重复步骤1-4，直到选取了三维云层场景内最后一个单元作为a点，且该处的多次散射源函数项计算完成为止。

本发明的有益效果为：本发明可在保持较高计算精度的前提下，加快计算三维云层辐射源函数，降低计算成本，并且可以很好地应用于三维多云场景。

附图说明

图1为本发明的方法流程示意图；

图2为本发明的水平方向单元划分示意图；

图3为本发明的高斯方向构成的三棱柱空间示意图；

图4为本发明的应用于所设云场景1随太阳天顶角变化的改进辐射源函数计算法和爱丁顿近似法与真实辐射值的相对误差结果图；

图5为本发明的应用于所设云场景2随太阳天顶角变化的改进辐射源函数计算法和爱丁顿近似法与真实辐射值的相对误差结果图。

具体实施方式

现将结合附图对本发明的技术方案进行完整的描述。以下描述仅仅是本发明的一部分实施案例而已，并非全部。基于本发明中的实施案例，本领域技术人员在没有作出创造性劳动的前提下所获得的所有其他实施案例，都属于本发明的权利保护范围之内。

如图1所示，一种基于辐射通量密度的三维云层辐射计算方法，包括如下步骤：

(1)利用高斯方向和大气属性构建水平辐射通量密度传输模型；

(2)利用爱丁顿近似法和严格单次散射计算式计算三维云层基础辐射场；

(3)考虑邻云辐射效应对水平传输方向上的三维云层基础辐射场的影响，将其分为邻云影响辐射增强和邻云影响辐射减弱，计算由邻云辐射效应影响所产生的水平辐射通量密度变化；

(4)将水平辐射通量密度变化映射回高斯方向，再将三维云层基础辐射场与变化辐射强度结合得到更新后三维云层辐射强度，并由其计算多次散射源函数；

(5)结合由初始条件计算的单次散射源函数项与热源函数项，得到三维云层辐射源函数。

1.邻云辐射效应

首先将三维云层场景分为所需大小尺度的网格单元，每个单元都有其相邻四个单元(如图2)，再由已知的云层含水量分布及米氏散射结果生成大气属性文件，计算对应网格位置处三维云层的反照率、消光系数和相函数，再将以上三类参量的大气属性存储在它们所对应位置的单元中。

通过爱丁顿近似计算得到多次散射项、严格单次散射计算式得到的单次散射项之后求得三维云层基础辐射场，再基于邻云辐射效应，计算得到三维云层变化辐射强度，结合三维云层基础辐射场和变化辐射强度，再通过辐射强度到源函数的转换，可以实现快速计算三维云层场景辐射源函数。

邻云辐射效应是由于光和热等辐射量在云层与邻云之间互相作用和交换而产生的辐射场分布特性改变，该效应可以准确而全面地对非均匀性较强的云层辐射进行描述。故通过计算与辐射相关参量的具体变化数值，可以对邻云辐射效应进行量化表达。三维云层空间内每个单元处都储存着由属性文件计算得到的辐射通量密度，考虑邻云辐射效应对辐射通量密度造成的影响，在网格单元内由影响造成的辐射通量密度变化分为邻云影响辐射增强和邻云影响辐射减弱(如图1)：

ice(τ；ω′ts)＝(1+ace)i⁺(τ；ω′ts)-(1-ace)i^-(τ；ω′ts)

其中ice(τ；ω′ts)代表某一单元处ω′ts方向由邻云辐射效应所产生的水平辐射通量密度变化值，i⁺(τ；ω′ts)为ω′ts方向水平辐射通量密度的增加值，i^-(τ；ω′ts)为ω′ts方向水平辐射通量密度的减弱值，ace为弥补太阳方向变化对模型带来误差的修正系数。

1.1映射到水平方向

对于三维云层辐射问题，邻云辐射效应可以表现为光束在云层之间传播路径的改变以及由此而产生的辐射能量变化，由邻云辐射效应所产生的水平辐射通量密度变化，可以计算多次散射源函数的近似值。辐射通量密度的方向为空间内的已设置高斯方向，分析邻云影响产生的水平辐射通量密度变化，需要将高斯方向辐射映射到水平方向上来，由高斯方向到水平方向的变换系数ψts：

ψts＝cos(δφts)sin(θts)

其中θts为高斯方向的天顶角；δφts为高斯方向在xy平面上的投影与其所处ts空间内水平坐标轴方向的平面夹角，δφts的取值范围是-45°＜δφts＜45°。

1.2邻云辐射影响增强和减弱

首先计算三维云层空间的基础辐射场，再找到第一个单元a，由其基础辐射强度及其相邻的单元b处的基础辐射强度计算得到受邻云辐射影响增强的辐射通量密度和减弱的辐射通量密度：

其中i⁺(τ；ω′t1)是ω′t1方向水平辐射通量密度的增加值；i-(τ；ω′t1)是ω′t1方向水平辐射通量密度的减弱值；ω′t2方向是与ω′t1相同高斯方向相反方位的方向(如图3)；r1是水平辐射通量密度在a与b1的交界处发生反射的比例；r2是水平辐射通量密度在a与b2的交界处发生反射的比例。

2辐射通量密度变化值计算

2.1修正系数计算

考虑太阳方向变化对模型带来的影响，从而弥补方法模型中存在的输入性误差。修正系数计算式如下：

ace＝exp(sinθ0-sinθe)-1

其中θ0为太阳天顶角；θe为二流近似法计算辐射通量所使用的高斯方向天顶角。

由于高斯积分的标准计算式可以写为：

其中f(γ)被积分函数；γ积分变量；aj高斯积分权重；γj高斯积分点。且在爱丁顿近似法等二流近似方法中，m的取值为1，此时高斯积分权重a1的取值为0.5773503，其为该高斯方向天顶角的余弦，因此可以求得(6)式中θe的数值约为54.74°。

2.2反照率r模型

找到三维空间网格内某一单元a，当在水平方向尺度上对邻云辐射效应进行分析时，考虑与其水平相邻的单元b1、b2、b3和b4，在三维空间含云区域单元交界处会发生反射，反照率r表示单元交界处发生反射的比例。当τ在0-20时，反照率r^*的计算参考邻云间反射率计算方法和二流近似中的反照率计算方法得到：

其中，τ为光学厚度；ω为单次散射反照率；g为不对称因子。

当τ不断增大时，r^*将逐渐趋近于1，此时无论ω和g取值如何，反照率的计算结果都会发生“饱和”，r^*的计算结果将会严重偏离真实值，针对该情况，对r^*进行修正得到适用于任意光学厚度的反照率r：

不对称因子可以描述前向散射和后向散射的程度，其计算方法为：

其中，ω为散射角；p(ω)为散射相函数。

光学厚度描述辐射传输过程中经过云层或大气吸收和散射产生的总衰弱：

其中，k为消光系数；s为路径距离。

2.3权重系数计算

将水平辐射通量密度分配至相应的高斯方向处，以实现该高斯方向上辐射强度变化的更新。辐射场辐射场从水平通量密度方向到高斯方向的分配权重为p(ω′ts；ωs)ats(ω′ts)。

其中，p(ω′ts；ωs)为水平方向到高斯方向的夹角作为散射角方式的散射相函数；ats(ω′ts)为高斯积分权重。

2.4映射回高斯方向

根据得到的邻云辐射影响增强和邻云辐射影响减弱，计算单元a处的水平方向变化辐射通量密度：

ice(τ；ω′ts)＝(1+ace)i⁺(τ；ω′ts)-(1-ace)i^-(τ；ω′ts)

其中ice(τ；ω′ts)代表某一单元处ω′ts方向由邻云辐射效应所产生的水平辐射通量密度变化值，i⁺(τ；ω′ts)为ω′ts方向水平辐射通量密度的增加值，i^-(τ；ω′ts)为ω′ts方向水平辐射通量密度的减弱值。

高斯方向上辐射通量密度更新为：

ice'(τ；ω′)＝p(ω′ts；ωs)ats(ω′ts)ice(τ；ω′)

于是，可以计算得到单元a处由于邻云辐射效应影响得到的变化辐射通量密度ice'，接着选取新的a点，重复1与2所述计算过程，直到选取了三维云层场景内最后一个单元作为a点，即所有单元的ice′都计算完成为止。

3.三维云层源函数计算

在所有a点的ice′都计算完成后，将三维云层基础辐射场与其结合，得到更新后辐射强度，并利用其计算多次散射源函数，进而结合单次散射源函数项与热源函数项，得到三维云层空间总辐射源函数。

3.1辐射强度累加

多次散射源函数近似计算模型以基础辐射场i′(τ；ω′)作为输入，基础辐射场由通过严格单次散射计算式得到的单次散射强度和利用爱丁顿近似法得到的多次散射强度构成。接下来，计算单元a处的更新后辐射强度：

i′(τ；ω′)＝iss(τ；ω′)+iems(τ；ω′)

i″(τ；ω′)＝i′(τ；ω′)+ice'(τ；ω′)

其中iss(τ；ω′)通过严格单次散射计算式得到的单次散射强度；iems(τ；ω′)通过爱丁顿近似法计算得到的多次散射强度。

3.2多次散射函数项计算

利用邻云辐射效应影响的变化辐射亮度更新后的辐射强度i”来计算单元a处的多次散射函数jms(τ；ω)：

其中p(ω；ω′)为三维云层空间单元a处的散射相函数。

3.3辐射源函数计算

在多次散射源函数jms计算完成后，将其与单次散射源函数jss和热源函数jt结合，得到辐射源函数j：

其中ω为反照率；f0为从辐射源产生并直接照射在单元上的辐射亮度；τ为大气光学厚度；b为普朗克公式；t为温度。

选取新的单元a，重复3.1、3.2与3.3所述计算过程，直到选取了三维云层场景内最后一个单元a，且该处的总辐射源函数计算完成为止。

4.辐射计算方法应用

为了验证本文提出的基于辐射通量密度的三维云层辐射计算方法，选用了真实云场景进行实验。该云场景为i3rcphaseⅱ的层积云场，网格尺寸为64*64*16，空间分辨率为55m*55m*25m，所处空间范围为x:0-3.52km,y:0-3.52km,z:10.1-11.6km。为测试方法在云层场景中的适用情况，截取了其中5*5*5的网格点数且云覆盖率较高的两块区域进行实验。考虑实验计算成本，将高斯方向离散数设为4，方位角离散数设为4，即在16个高斯方向上描述整个辐射场。

基于实验环境，将本文改进辐射源函数计算方法与爱丁顿近似法、shdom法进行比较分析。爱丁顿近似法是二流近似法的一种，其一维计算模式不考虑水平辐射变化，本方法在其基础上考虑水平辐射交换，与其比较可以明显观察到本方法的改进效果；而shdom方法虽然耗费成本较高，但计算结果趋于真实值，可以将其视为自然场景真实辐射值，与其比较能够分析本方法的误差。

在截取的云场景1中，云覆盖率为93.3％，在层积云场所处空间范围为x:1.98-2.2km,y:2.2-2.42km,z:10.9-11.3km。在截取的云场景2中，云覆盖率为69.44％，在层积云场所处空间范围为x:2.585-2.805km,y:2.805-3.025km,z:10.7-11.1km。通过三种方法的辐射计算，得到在不同太阳天顶角情况下爱丁顿近似法和改进辐射源函数计算法同真实辐射值的相对误差结果图。

结果中展示在上行空间本文辐射计算方法的平均相对误差明显低于爱丁顿近似法，云层场景1的平均精度提升为22.5％，云层场景2的平均精度提升为15％。对于下行空间而言，下行空间中存在太阳直射方向的前向方向，前向散射中包含了云层中的大多数单次散射能量，故下行空间内的辐射值大多高于上行空间，因此无论是爱丁顿近似还是本文方法，计算得到的的相对误差都较小，爱丁顿近似法与本方法皆趋近于真实辐射情况，本文方法的平均相对误差基本低于爱丁顿近似法，云层场景1的平均精度提升为2.8％，云层场景2的平均精度提升为1％。

在整个三维云层空间内，云层场景1的平均相对误差为8.75％，云层场景2的平均相对误差为11％，即在两种云场景中，方法误差皆在可接受范围内，可以满足一般应用需求。

本发明可在保持较高计算精度的前提下，加快计算三维云层辐射源函数，降低计算成本，并且可以很好地应用于三维多云场景。

以上实施例仅供说明本发明之用，而非对本发明的限制，有关技术领域的技术人员，在不脱离本发明的精神和范围的情况下，所作出各种变换或变型，均属于本发明的范畴。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前下，本发明还会有各种变化和改进，本发明要求保护范围由所附的权利要求书、说明书及其等效物界定。