基于蒙特卡洛算法的气溶胶多次散射模拟方法及系统与流程

本发明涉及卫星大气遥感技术领域，更为具体地，涉及一种基于蒙特卡洛算法的气溶胶多次散射模拟方法及系统。

背景技术：

随着我国工业的快速发展，污染物排放和城市悬浮物的大量增加使得部分地区发生灰霾天气的频率迅速增加，这种高气溶胶颗粒物浓度的大气不仅影响了空气质量、能见度以及人体健康，还加剧了大气的多次散射效应。因此，准确模拟气溶胶的多次散射过程是计算大气辐射传输精度的前提保证。

大气辐射传输是一门涵盖天文物理、应用光学和物理、行星空间科学与大气科学的既古老又获得蓬勃发展的综合学科。现有的计算大气辐射传输的方法有逐次散射法、球谐函数方法、离散纵坐标法、倍加-累加法、二通量法、四通量法、四流近似解等。

离散纵坐标法，是将辐射强度在空间方向上的变化进行离散，并且假设空间某一立体角内辐射强度是一个定值，从而将涉及角度离散的复杂积分形式或者积分-微分形式的辐射传输方程简化为简单的微分或者偏微分形式的线性微分方程。离散纵坐标法可以很好地使辐射传输方程与其它输运方程进行耦合求解、以及可以很方便的处理各向异性散射等。

倍加-累加法，是已知两个介质层的反射和透射性质，根据两层总的反射和透射性质通过两层之间的连续反射过程得到大气辐射传输方程。对于一个厚的均匀层，可以将其划分为2n个完全形同的薄层，通过倍加法可以迅速的获得整层的反射和透射性质。

球谐函数方法与离散纵坐标法类似，其是将辐射传输方程离散化，不同点在于该方法将辐射强度展开成球谐函数。球谐函数的二流形式即为Eddington近似。

逐次散射法，是通过迭代的方式求解辐射传输问题。理论上，逐次散射法 不仅可以用于平行大气，而且也可用于各种不均匀结构。

多通量理论可导出简洁准确的公式，计算简单并保证精度，常用于求解输运问题。在四通量理论中用前向和反向流动的准直射束以及前向和反向流动的扩散通量共四个量来代表强度流。

虽然上述计算大气辐射传输的算法都各有优点，但在解决一些复杂几何形状的辐射传输问题方面还是存在一定的局限性。

技术实现要素：

鉴于上述问题，本发明的目的是提供一种基于蒙特卡洛算法的气溶胶多次散射模拟方法及系统，以解决一些复杂几何形状的大气辐射传输问题。

根据本发明的一个方面，提供一种基于蒙特卡洛算法的气溶胶多次散射模拟方法，包括：

初始化光子权重、光子在大气顶层的随机位置和发射方向；

抽样获取光子在发生下一次碰撞前在大气中的传输距离；

根据所获取的传输距离和气溶胶散射参数库中的消光系数获取光子在物理空间中的传输距离，并根据所获取的光子在物理空间中的传输距离确定光子是否发生碰撞，气溶胶散射参数库中的参数包含每种气溶胶组分的散射相函数、单次散射反照率和消光系数；其中，

在光子发生碰撞时，根据光子的碰撞类型确定经过碰撞后的光子的新权重；其中，光子的碰撞类型包括散射碰撞和地表反射碰撞；

将光子的新权重与预设阈值进行比较；其中，

在光子的新权重不小于预设阈值时，计算光子发生碰撞后的新的传输方向，并重新抽样获取光子在发生下一次碰撞前在大气中的传输距离；

在光子的新权重小于预设阈值或者光子未发生碰撞时，将追踪的光子个数与预设的光子个数进行比较；其中，

在追踪的光子个数小于预设的光子个数时，重新初始化光子权重、光子在大气顶层的随机位置和发射方向；

在追踪的光子个数与预设的光子个数相同时，进行光子误差统计，获取大气透过率和大气顶层平均反射率。

根据本发明的另一方面，提供一种基于蒙特卡洛算法的气溶胶多次散射模拟系统，包括：

初始化单元，用于初始化光子权重、光子在大气顶层的随机位置和发射方向；

传输距离获取单元，用于抽样获取光子在发生下一次碰撞前在大气中的传输距离，根据所获取的距离和气溶胶散射参数库中的消光系数获取光子在物理空间中的传输距离，其中，气溶胶散射参数库中的参数包含每种气溶胶组分的散射相函数、单次散射反照率和消光系数；

光子新权重获取单元，用于在根据传输距离获取单元获取的光子在物理空间中的传输距离确定光子发生碰撞后，根据光子的碰撞类型获取经过碰撞后的光子的新权重；其中，光子的碰撞类型包括散射碰撞和地表反射碰撞；

光子传输方向获取单元，用于计算光子发生碰撞后的新的传输方向，并重新抽样获取光子在发生下一次碰撞前在大气中的传输距离；其中，当光子新权重获取单元获取的光子的新权重比预设阈值大时，光子传输方向获取单元获取光子发生碰撞后的新的传输方向；

比较单元，用于将追踪的光子个数与预设的光子个数进行比较；其中，当光子的新权重小于预设阈值或者光子未发生碰撞时，比较单元将追踪的光子个数与预设的光子个数进行比较；其中，当追踪的光子个数小于预设的光子个数时，初始化单元重新初始化光子权重、光子在大气顶层的随机位置和发射方向；

统计单元，用于进行光子误差统计，获取大气透过率和大气顶层平均反射率；其中，当光子新权重获取单元获取的光子的新权重比预设阈值小，且追踪的光子个数达到预设的光子个数时，统计单元进行光子误差统计。

利用上述根据本发明的基于蒙特卡洛算法的气溶胶多次散射模拟方法及系统，不仅能够处理任意单次散射反照率和各向异性很强的散射相函数，还能够处理任意几何形状下的大气辐射传输，以及对其他大气辐射传输算法的计算结果进行验证，进而准确地对气溶胶的多次散射过程进行模拟。

为了实现上述以及相关目的，本发明的一个或多个方面包括后面将详细 说明并在权利要求中特别指出的特征。下面的说明以及附图详细说明了本发明的某些示例性方面。然而，这些方面指示的仅仅是可使用本发明的原理的各种方式中的一些方式。此外，本发明旨在包括所有这些方面以及它们的等同物。

附图说明

通过参考以下结合附图的说明及权利要求书的内容，并且随着对本发明的更全面理解，本发明的其它目的及结果将更加明白及易于理解。在附图中：

图1为根据本发明实施例的基于蒙特卡洛算法的气溶胶多次散射模拟方法的第一流程示意图；

图2为根据本发明实施例的基于蒙特卡洛算法的气溶胶多次散射模拟方法的第二流程示意图；

图3为根据本发明实施例的基于蒙特卡洛算法的气溶胶多次散射模拟系统的逻辑结构示意图。

在所有附图中相同的标号指示相似或相应的特征或功能。

具体实施方式

以下将结合附图对本发明的具体实施例进行详细描述。

针对前述现有的计算大气辐射传输的算法在解决复杂几何形状的辐射传输问题方面存在局限性的问题，本发明以蒙特卡洛概率抽样为基础，把光子的传输过程转化为随机问题来求解，通过本发明能够较为准确地对大气气溶胶多次散射过程进行模拟，有利于进一步提高在高颗粒物浓度条件下的大气辐射传输计算精度。

在对本发明进行说明之前，先对蒙特卡洛算法进行一个初步介绍。

蒙特卡洛算法是将光子的散射看成一个随机过程，以辐射传输方程为依据，并直接模拟辐射传输方程。它是一种随机模拟方法，即将散射过程当成光子和介质中的碰撞过程，两次碰撞之间光子在介质所走的距离与消光系数有关，碰撞后光子将改变前进方向，散射角由相函数确定，对大量光子“行为”跟踪并进行统计就可得具体问题的结果。蒙特卡洛算法的优点是：其可以解 决一些复杂几何的辐射传输问题或作为其他辐射传输算法的结果验证，其能够处理任意几何形状下的辐射传输问题，也能够处理任意单次散射反照率和各向异性很强的散射相函数，而现有的其他辐射传输算法在这方面都具有一定的局限性。

为了说明本发明提供的基于蒙特卡洛算法的气溶胶多次散射模拟方法，图1示出了根据本发明实施例的基于蒙特卡洛算法的气溶胶多次散射模拟方法的第一流程。

如图1所示，本发明提供的基于蒙特卡洛算法的气溶胶多次散射模拟方法包括：

S110：计算并建立气溶胶散射参数库；其中，该气溶胶散射参数库中的参数包含每种气溶胶组分的散射相函数、单次散射反照率和消光系数。

S120：初始化光子权重、光子在大气顶层的随机位置和发射方向。

S130：抽样获取光子在发生下一次碰撞前在大气中的传输距离。

S140：根据所获取的传输距离和气溶胶散射参数库中的消光系数获取光子在物理空间中的传输距离，并根据所获取的光子在物理空间中的传输距离确定光子是否发生碰撞。

S150：当光子发生碰撞时，根据光子的碰撞类型确定经过碰撞后的光子的新权重；其中，光子的碰撞类型包括散射碰撞和地表反射碰撞。

S160：将光子的新权重与预设阈值进行比较；其中，当光子的新权重不小于预设阈值时，进入步骤S170；当光子的新权重小于预设阈值或者光子未发生碰撞时，进入步骤S180。

S170：计算光子发生碰撞后的新的传输方向，并返回S130。

S180：将追踪的光子个数与预设的光子个数进行比较；

S190：当追踪的光子个数小于预设的光子个数时，返回S120；当追踪的光子个数与预设的光子个数相同时，进行光子误差统计，获取大气透过率和大气顶层平均反射率。

为了更为详细的对本发明进行说明，图2示出了根据本发明实施例的基于蒙特卡洛算法的气溶胶多次散射模拟方法的第二流程。如图2所示：

S1：计算并建立气溶胶散射参数库。

具体地，基于Mie算法，获取不同组分、不同粒径大小的气溶胶的散射 相函数、单次散射反照率和消光截面；将所获取的不同组分、不同粒径大小的气溶胶的散射相函数、单次散射反照率和消光截面与粒子谱分布函数进行积分，建立包含每种气溶胶组分的散射相函数、单次散射反照率和消光系数的散射参数库。

例如：利用Mie散射方法，计算不同半径为0.01～15微米的黑碳、硫酸盐、有机碳、沙尘、海盐等气溶胶粒子的散射相函数、单次散射反照率和消光截面；将上述计算结果与正态分布的粒子谱分布函数进行积分，建立包含黑碳、硫酸盐、有机碳、沙尘、海盐气溶胶的散射相函数、单次散射反照率和消光系数的散射参数库。

S2：初始化光子权重、光子在大气顶层的随机位置、发射方向。

具体地，将光子的初始权重设为1，将光子在大气顶层的初始位置设置为在大气顶层x–y平面上的任意随机位置，基于当地的时间和经纬度，由太阳天顶角和太阳方位角确定光子的发射方向。

例如，在假设太阳是在均一照射的前提下，对于由太阳光发射的光子，将其初始位置设为在大气顶层x–y平面上的任意随机位置。即：x₀＝ξ_x(x_max-x_min)+x_min；y₀＝ξ_y(y_max-y_min)+y_min；z₀＝z_max；

其中，ξ_x和ξ_y均为在[0，1]范围内均匀分布的随机数，[x_min，x_max]、[y_min，y_max]、[z_min，z_max]分别为x、y和z方向的大气单元格的取值范围，其中x_min＝0km、x_max＝10km、y_min＝0km、y_max＝10km、z_min＝0km、z_max＝12km。

至于光子的发射方向，则可以由太阳天顶角30°和太阳方位角150°来进行描述，光子的发射方向与当地时间和经纬度有关，而不是一个随机过程。

S3：模拟光子在下一次发生散射碰撞前在大气中的传输距离，并判断光子是否逃逸出计算区域。

具体地，随机抽样获取光子在下一次发生散射碰撞前在大气中的传输距离所对应的光学厚度τ，即τ＝-ln(1-ξ)，其中，ξ为在[0，1]范围内均匀分布的随机数；基于步骤S1中计算得到的气溶胶消光系数β，获取光子在物理空间中的传输距离S，即S＝τ/β；判断S是否大于光子当前位置沿传输方向距计算区域边界的距离(即预设的计算区域)。

其中，当获取的光子在物理空间中的传输距离S小于预设的计算区域时， 说明光子未逃逸出预设的计算区域，此时光子发生散射碰撞，进入步骤S4；当获取的光子在物理空间中的传输距离S不小于预设的计算区域时，说明光子已逃逸出预设的计算区域，进入步骤S5。

S4：将光子移动至碰撞位置，并判断光子的碰撞类型，然后进入S6。

具体地，将光子移动至发生散射碰撞的位置(x_new，y_new，z_new)；建立黑碳、硫酸盐、有机碳、沙尘、海盐气溶胶的消光系数的计概率查找表；由在[0，1]范围内均匀分布的随机数ξ，在已有的消光系数累积概率查找表中，依次查找光子的散射碰撞类型，直至确定随机数ξ所对应的散射碰撞类型，以确定与光子发生散射碰撞的气溶胶组分。

S5：判断光子是否发生地表反射碰撞，若没有发生地表反射碰撞，则光子在预设的计算区域中的传输结束，进入步骤S9；若发生地表反射碰撞，则进入步骤S6。

S6：计算光子经过碰撞后的新权重。

具体地，若光子发生散射碰撞，那么光子的新权重为当前权重值与所选择碰撞的气溶胶组分的单次散射反照率之积；若光子发生地表反射，那么光子的新权重为当前权重值与地表反照率之积。

需要说明的是，地表反照率为0.3；上述所选择碰撞的气溶胶组分即为确定的与光子发生散射碰撞的气溶胶组分。

S7：判断光子的新权重是否小于设定阈值。其中，若光子的新权重小于设定阈值，则光子传输在计算区域中的传输结束，进入步骤S9；若光子的新权重不小于设定阈值，则进入S8。

S8：计算光子发生碰撞后新的传输方向，并返回步骤S3。

例如，设定阈值为0.5，那么当光子的新权重小于0.5时，选择一随机数ξ与光子的新权重进行比较，若光子的权重不大于该随机数ξ，则该光子传输终结，停止对该光子的追踪；若光子的权重大于该随机数ξ，则将光子的权重重置为1，并继续其传输过程；若光子发生散射碰撞，计算选择碰撞的气溶胶组分散射相函数的累计概率，并由在[0，1]范围内均匀分布的随机数ξ抽样选择散射角；若光子发生地表反射碰撞，则反射角余弦值为反射方位角为2px₂，其中，x₁和x₂均在[0，1]范围内均匀分布的随机数。

其中，上述随机数均为由梅森旋转算法生成的在[0，1]范围内均匀分布的伪随机数，上述累积概率为概率密度分布函数的积分。

S9：判断追踪的光子个数是否达到初始设定的光子个数。例如，可以将初始设定的光子个数设为10000000；当追踪的光子个数达到该设定的光子个数时，进入步骤S10，否则进入步骤S2。

S10：进行光子误差统计，并计算大气透过率和大气顶层平均反射率。

与上述方法相对应，本发明提供一种基于蒙特卡洛算法的气溶胶多次散射模拟系统。图3示出了根据本发明实施例的基于蒙特卡洛算法的气溶胶多次散射模拟系统的逻辑结构。

如图3所示，本发明的提供的基于蒙特卡洛算法的气溶胶多次散射模拟系统300包括气溶胶散射参数库建立单元310、初始化单元320、传输距离获取单元330、光子新权重获取单元340、光子传输方向获取单元350、比较单元360、统计单元370。

其中，气溶胶散射参数库建立单元310建立气溶胶散射参数库；其中，气溶胶散射参数库中的参数包含每种气溶胶组分的散射相函数、单次散射反照率和消光系数。

具体地，基于Mie算法，获取不同组分、不同粒径大小的气溶胶的散射相函数、单次散射反照率和消光截面；将所获取的不同组分、不同粒径大小的气溶胶的散射相函数、单次散射反照率和消光截面与粒子谱分布函数进行积分，建立包含每种气溶胶组分的散射相函数、单次散射反照率和消光系数的散射参数库。

初始化单元320用于初始化光子权重、光子在大气顶层的随机位置和发射方向。

具体地，设定光子初始权重为1，设定光子在大气顶层的初始位置为在大气顶层x–y平面上的任意随机位置；基于当地的时间和经纬度，由太阳天顶角和太阳方位角确定光子的发射方向。

传输距离获取单元330用于抽样获取光子在发生下一次碰撞前在大气中的传输距离，根据所获取的距离和气溶胶散射参数库中的消光系数获取光子在物理空间中的传输距离。

光子新权重获取单元340用于在根据传输距离获取单元获取的光子在物理 空间中的传输距离确定光子发生碰撞后，根据光子的碰撞类型获取经过碰撞后的光子的新权重；其中，光子的碰撞类型包括散射碰撞和地表反射碰撞。

具体地，根据传输距离获取单元330所获取的光子在物理空间中的传输距离确定光子是否发生碰撞，当获取的光子在物理空间中的传输距离小于预设的计算区域时，光子则发生散射碰撞，将光子移动至发生散射碰撞的位置；根据参数库中的消光系数建立所有气溶胶组分的消光系数累计概率查找表；由在[0，1]范围内均匀分布的随机数ξ，在已有的累计概率查找表，依次查找光子的散射碰撞类型，直至确定随机数ξ所对应的散射碰撞类型，以确定与光子发生散射碰撞的气溶胶组分；将与光子发生碰撞的气溶胶组分的单次散射反照率和光子当前的权重之积作为光子的新权重。

当获取的光子在物理空间中的传输距离不小于预设的计算区域，且发生地表反射碰撞时，将光子移动至所述计算区域的边界位置；将光子的当前权重与地表反照率之积作为光子的新权重。

光子传输方向获取单元350用于计算光子发生碰撞后的新的传输方向，并重新抽样获取光子在发生下一次碰撞前在大气中的传输距离；其中，当光子新权重获取单元获取的光子的新权重不比预设阈值小时，光子传输方向获取单元350获取光子发生碰撞后的新的传输方向。

具体地，当光子的新权重小于0.5时，选择一随机数与光子的新权重进行比较；其中，当光子的新权重不大于所选择的随机数时，光子传输结束，停止对光子的追踪，判断追踪的光子的个数是否达到预设的光子个数；当光子的新权重大于所选择的随机数时，将光子的权重重置为1，光子继续传输。

在光子继续传输的过程中，当光子发生散射碰撞时，基于气溶胶散射参数库中的与所确定的散射碰撞的气溶胶组分相对应的散射相函数，获取与光子相碰撞的气溶胶组分的散射相函数的累计概率，并由在[0，1]范围内均匀分布的随机数ξ抽样选择散射角；当光子发生地表反射碰撞时，则反射角余弦值为反射方位角为2px₂；其中，x₁和ξ₂均为在[0，1]范围内均匀分布的随机数。其中，上述随机数均为由梅森旋转算法生成的在[0，1]范围内均匀分布的伪随机数，上述累积概率为概率密度分布函数的积分。

比较单元360用于将追踪的光子个数与预设的光子个数进行比较；其中，当光子的新权重小于预设阈值或者光子未发生碰撞时，比较单元将追踪的光 子个数与预设的光子个数进行比较；其中，当追踪的光子个数小于预设的光子个数时，初始化单元重新初始化光子权重、光子在大气顶层的随机位置和发射方向。

统计单元370用于进行光子误差统计，获取大气透过率和大气顶层平均反射率；其中，当光子新权重获取单元获取的光子的新权重比预设阈值小，且追踪的光子个数达到预设的光子个数时，统计单元进行光子误差统计。

根据上述，本发明提供的基于蒙特卡洛算法的气溶胶多次散射模拟方法及系统是将光子传输过程分解为发射、吸收、散射和反射等一系列子过程，通过建立每个子过程的概率模型，将其转化为随机问题来求解，因此，本发明不仅能够处理任意单次散射反照率和各向异性很强的散射相函数，还能够处理任意几何形状下的大气辐射传输，以及对其他大气辐射传输算法的计算结果进行验证，进而准确地对气溶胶的多次散射过程进行模拟。

如上参照附图以示例的方式描述了根据本发明的基于蒙特卡洛算法的气溶胶多次散射模拟方法及系统。但是，本领域技术人员应当理解，对于上述本发明所提出的基于蒙特卡洛算法的气溶胶多次散射模拟方法及系统，还可以在不脱离本发明内容的基础上做出各种改进。因此，本发明的保护范围应当由所附的权利要求书的内容确定。