大气污染物垂直分布的实时反演方法、装置、计算机设备和存储介质与流程

1.本发明涉及光学测量方法领域，特别是涉及一种大气污染物垂直分布的实时反演方法、装置、计算机设备和存储介质。


背景技术：

2.各种污染源直接排放的颗粒物、二次颗粒物和污染气体大量增加，使得对流层大气气溶胶和污染气体垂直分布存在更加快速和剧烈的变化。大气气溶胶和部分痕量气体等大气污染物对人体和生态环境有重要危害。
3.其中，no2是主要的的大气污染痕量气体之一，其参与臭氧的生成的光化学反应和气溶胶的生成。no2和大气气溶胶对人体的健康均有危害。大气中的hcho是重要的致癌物，so2会形成酸雨，hono对大气氧化性有重要影响。因此，对大气气溶胶和痕量气体等大气污染物的浓度和分布进行长期的监测，对于认识大气污染物的演化过程、校准卫星观测数据和质量模型具有重要意义。
4.现有技术中，采用被动光源的超光谱遥感技术可以对大气气溶胶和痕量气体的垂直廓线进行长时间、高时间分辨率的监测。采用差分光谱技术可以观测在大气中分布接近恒定的氧二聚体o4，以此来反映大气气溶胶对氧二聚体o4吸收峰光路的影响，从而反演气溶胶粒子属性，也可以直接测定大气中痕量气体(如no2、so2、hcho、chocho、hono、o3等)的浓度垂直分布，常用的大气污染物垂直浓度分布反演方法为最优化迭代估计方法和查找表方法。
5.其中，查找表方法计算速度快，但反演的大气污染物浓度垂直分布不稳定，易出反常规结果。最优化迭代估计方法反演的大气污染物浓度垂直分布较为稳定，但通常使用数值差分方式计算迭代所需的权重函数，具体过程为：需要使用已知先验廓线模拟一次光路总浓度，视为标准光路总浓度，紧接着需要将廓线某一层的浓度增加一个浓度增量，再进行光路总浓度模拟，视为波动光路总浓度。将波动光路总浓度与标准光路总浓度做差并除以浓度增量，即可获得对应该层的光路总浓度权重函数。将所有层都增加浓度增量模拟光路总浓度即可获得完整的光路总浓度权重函数。因为每一层的改变求取光路总浓度均需一次模拟，再加上模拟标准光路总浓度需要求解一次辐射传输模型，差分方法获得光路总浓度权重函数总共需要廓线层数 1次辐射传输模型模拟光路总浓度，这极大地降低了最优化迭代估计方法的时间效率，无法实现实时输出大气污染物垂直分布。
6.因此，最优化估计迭代反演方法和查找表方法存在的问题限制了超光谱遥感技术在大气污染物实时监测和溯源方面的应用。因此发展一种稳定且实时的超光谱遥感反演技术对实现大气污染物的实时监测和溯源具有重要意义。


技术实现要素：

7.鉴于上述技术问题，本发明的目的在于提供大气污染物垂直分布的实时反演方
法、装置、计算机设备和存储介质，降低反演计算开销，实现对大气污染物垂直分布的在线实时反演监测，同时保证反演监测的稳定性。
8.为实现上述发明目的，本发明提供以下技术方案：
9.第一方面，实施例提供了一种大气污染物垂直分布的实时反演方法，包括以下步骤：
10.步骤1，获取利用超光谱遥感技术实时采集的多仰角观测的高分辨率光谱，根据高分辨率光谱拟合得到大气污染物多仰角观测的实测光路总浓度；
11.步骤2，采用辐射传输模型模拟地面光强得到模拟光强，并结合大气污染物特征吸收截面将模拟光强转换为模拟光路总浓度；
12.步骤3，计算模拟光强对分层的大气污染物浓度分布的偏微分以得到模拟光强的雅克比矩阵，利用模拟光强的雅克比矩阵并结合大气污染物特征吸收截面计算模拟光路总浓度的雅克比矩阵，即迭代权重函数；
13.步骤4，以最小化带迭代权重函数的实测光路总浓度和模拟光路总浓度为目标函数，利用采用基于最小二乘法的最优化迭代估计方法迭代求解以得到大气污染物垂直分布。
14.在一个实施例中，步骤1中，在获得多仰角观测的高分辨率光谱后，从高分辨光谱中过滤掉随波长做慢变化的光谱结构，将剩余光谱与大气污染物的标准参考光谱进行最小二乘法拟合，以得到每种大气污染物每个角度观测的实测光路总浓度。
15.在一个实施例中，步骤2中，通过采用辐射传输模型模拟地面光强得到模拟光强包括目标污染气体存在于大气中时的模拟光强i1和目标污染气体不存在于大气中时的模拟光强i0，利用模拟光强i1和模拟光强i0，并结合大气污染物特征吸收截面σ计算模拟光路总浓度dscd
model
，其表达式为：
[0016][0017]
其中，模拟光强i的上标表示观测仰角，其中，90表示90度观测仰角，β为除90度之外的其它观测仰角，和分别表示目标污染气体不存在于大气中时以β度观测仰角和90度观测仰角观测时的模拟光强，和分别表示目标污染气体存在于大气中时以β度观测仰角和90度观测仰角观测时的模拟光强。
[0018]
在一个实施例中，步骤3中，采用以下公式计算计算模拟光强对分层的大气污染物浓度分布的偏微分以得到模拟光强的雅克比矩阵j：
[0019][0020]
其中，模拟光强i包括目标污染气体存在于大气中时的模拟光强i1和目标污染气体不存在于大气中时的模拟光强i0，c表示每层目标污染气体的浓度。
[0021]
在一个实施例中，步骤3中，利用模拟光强的雅克比矩阵并结合大气污染物特征吸收截面计算模拟光路总浓度的雅克比矩阵，即迭代权重函数 wf表示为：
[0022][0023]
其中，上标表示观测仰角，其中，90表示90度观测仰角，β为除90 度之外的其它观测仰角，下标表示目前污染物是否存在于大气中，其中， 0表示目前污染物是不存在于大气中，1表示目前污染物是存在于大气中，和分别表示目标污染气体不存在于大气中时以β度观测仰角和90度观测仰角观测时的模拟光强，和分别表示目标污染气体存在于大气中时以β度观测仰角和90度观测仰角观测时的模拟光强，分别表示模拟光强的雅克比矩阵，σ表示大气污染物特征吸收截面。
[0024]
在一个实施例中，步骤4中，利用迭代权重函数并采用基于最小二乘法的最优化迭代估计方法迭代求解时，获得的目标大气污染气体垂直浓度分布迭代解为：
[0025][0026]
其中，c表示污染气体的浓度或消光系数等，下标n表示迭代步数，下标a表示反演所需先验信息，sa表示先验误差矩阵，表示迭代权重函数，表示实测dscd误差矩阵，dscd
meas
表示实测光路总浓度， dscd
model
表示模拟光路总浓度。γn是迭代因子，当γn大时迭代接近于梯度下降法，当γn小时迭代接近于高斯牛顿法。
[0027]
在一个实施例中，步骤4中，迭代求解过程为：
[0028]
采用初始先验廓线和初始γ计算初始目标函数，然后采用初始目标函数使用迭代式(4)计算迭代下一步的状态矢量c
next
，该状态矢量c
next
为污染气体的浓度或消光系数，并计算状态矢量c
next
对应的新目标函数；
[0029]
当新目标函数较前一步的目标函数减小，且减小值大于预设值，且没有达到最大迭代次数时，减小γ，根据状态矢量c
next
采用迭代式(4)进行迭代计算，继续循环；若减小值小于等于预设值，则以状态矢量c
next
为最优解并输出；
[0030]
当新目标函数较前一步的目标函数没减小时，若γ大于预设γ值，停止循环，此时无法找到最优解，若γ小于等于预设γ值，增大γ值后，使用新γ采用迭代式(4)重新计算c
next
，继续循环。
[0031]
第二方面，实施例提供了一种大气污染物垂直分布的实时反演装置，包括：
[0032]
实测光路总浓度计算模块，用于获取利用超光谱遥感技术实时采集的多仰角观测的高分辨率光谱，根据高分辨率光谱拟合得到大气污染物多仰角观测的实测光路总浓度；
[0033]
模拟光路总浓度计算模块，用于采用辐射传输模型模拟地面光强得到模拟光强，并结合大气污染物特征吸收截面将模拟光强转换为模拟光路总浓度；
[0034]
迭代权重函数计算模块，用于计算模拟光强对分层的大气污染物浓度分布的偏微分以得到模拟光强的雅克比矩阵，利用模拟光强的雅克比矩阵并结合大气污染物特征吸收
截面计算模拟光路总浓度的雅克比矩阵，即迭代权重函数；
[0035]
迭代计算模块，以最小化带迭代权重函数的实测光路总浓度和模拟光路总浓度为目标函数，利用采用基于最小二乘法的最优化迭代估计方法迭代求解以得到大气污染物垂直分布。
[0036]
第三方面，实施例提供了一种计算机设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上执行的计算机程序，所述处理器执行所述计算机程序时实现第一方面所述的大气污染物垂直分布的实时反演方法的步骤。
[0037]
第四方面，实施例提供了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理执行时实现第一方面所述的大气污染物垂直分布的实时反演方法的步骤。
[0038]
与现有技术相比，实施例提供的大气污染物垂直分布的实时反演方法、装置、计算机设备和存储介质具有的有益效果至少包括：
[0039]
通过步骤2和步骤3两步骤，只需要使用两次辐射传输模型的模拟光强i即可以求出模拟光路总浓度及迭代权重函数，这大大减少了求解迭代权重函数的计算开销和消耗时长，提高了计算效率，保证了反演的实时性，在此基础上，通过步骤4，采用基于最小二乘法的最优化迭代估计方法迭代求解大气污染物垂直分布，保证了反演污染物浓度垂直分布稳定性。
附图说明
[0040]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动前提下，还可以根据这些附图获得其他附图。
[0041]
图1是一实施例提供的大气污染物垂直分布的实时反演方法的流程图；
[0042]
图2是一实施例提供的迭代求解过程的流程图；
[0043]
图3是一实施例提供的气溶胶消光系数廓线，其中，虚线为气溶胶消光系数的先验廓线，实线表示采用本发明提供的实时反演方法得到的反演廓线；
[0044]
图4是一实施例提供的o4的模拟光路总浓度的对比图，其中， 号代表观测值，*号代表模拟值；
[0045]
图5是一实施例提供的采用不同方法对aeronet气溶胶总量反演结果的对比验证图，其中，十字代表基于差分求解权重函数反演结果，方框代表本发明提供的实时反演方法的反演结果；
[0046]
图6是一实施例提供的大气污染物垂直分布的实时反演装置的结构示意图。
具体实施方式
[0047]
为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。
[0048]
稳定且实时的超光谱遥感反演技术可以实现大气污染物的实时监测和溯源，有利
于研究污染物的来源、途径及污染物在大气中发生的物理化学变化，能够为污染物的防控提供科学指导。
[0049]
图1是一实施例提供的大气污染物垂直分布的实时反演方法的流程图。如图1所示，实施例提供的大气污染物垂直分布的实时反演方法，包括以下步骤：
[0050]
步骤1，获取利用超光谱遥感技术实时采集的多仰角观测的高分辨率光谱，根据高分辨率光谱拟合得到大气污染物多仰角观测的实测光路总浓度。
[0051]
在进行高分辨率光谱采集时，采用高分辨率光谱仪进行多仰角观测角度的天空散射光谱采集，以得到多仰角观测的高分辨率光谱。实施例中，可以设置仰角序列为1、2、3、4、5、6、8、10、15、30、90度。限定每个仰角的高分辨率光谱的采集时间为1分钟，则针对所有仰角的高分辨率光谱的一轮总采集时间即为x分钟，x为采集仰角的个数。为了实现大气污染物垂直分布的实时反演，要求依据高分辨率光谱反演得到大气污染物垂直分布的时间要小于每轮总采集时间。
[0052]
实施例采用差分的思想，进行实测光路总浓度的拟合。具体地，将大气污染物的吸收截面变为随波长作慢变化的部分σb(λ)和快变化部分σ
′
(λ), 那么，差分光学厚度为：
[0053][0054]
其中，i0′
是未经过气体吸收的光强，i是经过气体吸收的光强，∑σi′
(λ)dscdi就是差分光学厚度。
[0055]
在此基础上，从高分辨光谱中过滤掉随波长做慢变化的光谱结构，将剩余光谱与大气污染物的标准参考光谱(也即是特征吸收光谱)进行最小二乘法拟合，以得到每种大气污染物每个角度观测的实测光路总浓度。
[0056]
步骤2，采用辐射传输模型模拟地面光强得到模拟光强，并结合大气污染物特征吸收截面将模拟光强转换为模拟光路总浓度。
[0057]
实施例中，通过采用辐射传输模型模拟地面光强得到模拟光强包括目标污染气体存在于大气中时的模拟光强i1和目标污染气体不存在于大气中时的模拟光强i0，利用模拟光强i1和模拟光强i0，并结合大气污染物特征吸收截面σ计算模拟光路总浓度dscd
model
，其表达式为：
[0058][0059]
其中，模拟光强i的上标表示观测仰角，其中，90表示90度观测仰角，β为除90度之外的其它观测仰角，和分别表示目标污染气体不存在于大气中时以β度观测仰角和90度观测仰角观测时的模拟光强，和分别表示目标污染气体存在于大气中时以β度观测仰角和90度观测仰角观测时的模拟光强。
[0060]
步骤3，计算模拟光强对分层的大气污染物浓度分布的偏微分以得到模拟光强的雅克比矩阵，利用模拟光强的雅克比矩阵并结合大气污染物特征吸收截面计算模拟光路总浓度的雅克比矩阵，即迭代权重函数。
[0061]
实施例中，采用以下公式计算计算模拟光强对分层的大气污染物浓度分布的偏微
分以得到模拟光强的雅克比矩阵j：
[0062][0063]
其中，模拟光强i包括目标污染气体存在于大气中时的模拟光强i1和目标污染气体不存在于大气中时的模拟光强i0，c表示每层目标污染气体的浓度。
[0064]
随着模拟光强i的上下标改变时，雅克比矩阵j同理也有相同上下标。基于此，利用模拟光强的雅克比矩阵并结合大气污染物特征吸收截面计算模拟光路总浓度的雅克比矩阵，即迭代权重函数wf表示为：
[0065][0066]
其中，上标表示观测仰角，其中，90表示90度观测仰角，β为除90 度之外的其它观测仰角，下标表示目前污染物是否存在于大气中，其中， 0表示目前污染物是不存在于大气中，1表示目前污染物是存在于大气中，和分别表示目标污染气体不存在于大气中时以β度观测仰角和90度观测仰角观测时的模拟光强，和分别表示目标污染气体存在于大气中时以β度观测仰角和90度观测仰角观测时的模拟光强，分别表示模拟光强的雅克比矩阵，σ表示大气污染物特征吸收截面。
[0067]
步骤4，以最小化带迭代权重函数的实测光路总浓度和模拟光路总浓度为目标函数，采用基于最小二乘法的最优化迭代估计方法迭代求解以得到大气污染物垂直分布。
[0068]
实施例中，在获得迭代权重函数wf后，构建目标函数为最小化带迭代权重函数的实测光路总浓度dscd
meas
和模拟光路总浓度dscd
model
，然后，采用基于最小二乘法的最优化迭代估计方法迭代获取实测光路总浓度 dscd
meas
和模拟光路总浓度dscd
model
的最佳你和状态，通过求解以下公式以得到大气污染物垂直分布。
[0069][0070]
其中，c表示污染气体的浓度或消光系数等，下标n表示迭代步数，γ表示迭代参数，下标a表示反演所需先验信息，sa表示先验误差矩阵，表示迭代权重函数，表示实测dscd误差矩阵，dscd
meas
表示实测光路总浓度，dscd
model
表示模拟光路总浓度。γn是迭代因子，当γn大时迭代接近于梯度下降法，当γn小时迭代接近于高斯牛顿法。
[0071]
图2是一实施例提供的迭代求解过程的流程图。如图2所示，非线性的最小二乘法迭代求解过程为：
[0072]
采用初始先验廓线(先验廓线通过污染物总量和标高计算得到)和初始迭代因子γ计算初始目标函数。然后采用初始目标函数使用迭代式(4) 计算迭代下一步的状态矢量c
next
，该状态矢量c
next
为污染气体的浓度或消光系数，并计算状态矢量c
next
对应的新目标函数；
[0073]
当新目标函数较前一步的目标函数减小，且减小值大于预设值，且没有达到最大迭代次数时，减小γ，根据状态矢量c
next
采用迭代式(4)进行迭代计算，继续循环；若减小值小于等于预设值，则以状态矢量c
next
为最优解并输出；
[0074]
当新目标函数较前一步的目标函数没减小时，若γ大于预设γ值，停止循环，此时无法找到最优解，若γ小于等于预设γ值，增大γ值后，使用新γ采用迭代式(4)重新计算c
next
，继续循环。
[0075]
在有效观测高度为4km时，表1给出利用上述大气污染物垂直分布的实时反演方法(也即是雅克比方法)和差分方法进行反演的时间效率对比结果，分析表1可得，除4000米之前，雅克比方法耗时与差分耗时的比值均小于1，表明雅克比矩阵方法耗时短于差分耗时，且随着垂直分辨高度的递减，雅克比矩阵方法耗时越来越短于差分耗时，计算效率逐渐增强。
[0076]
表1
[0077][0078]
图3是一实施例提供的气溶胶消光系数廓线，其中，虚线为气溶胶消光系数的先验廓线，实线表示采用本发明提供的实时反演方法得到的反演廓线。分析图3可得，实际消光系数廓线相较先验廓线形态类似但值偏大，实际消光廓线的垂直分布为e指数型，符合常规的对消光廓线垂直分布的总体特征。
[0079]
图4是一实施例提供的o4的光路总浓度的对比图，其中， 号代表观测值，*号代表模拟值。分析图4可得，o4的模拟和实测的光路总浓度可以较好的符合，说明在图3的气溶胶消光廓线情况下，实测与预期的o4的光路总浓度相近，表明了采用本发明提供的反演方法反演的气溶胶消光廓线与实际消光廓线相近。
[0080]
图5是一实施例提供的采用不同方法对aeronet气溶胶总量反演结果的对比验证图，其中，十字代表基于差分求解权重函数反演结果，方框代表本发明提供的实时反演方法的反演结果。分析图5可得，aeronet 实测的气溶胶总量与高光谱反演所得的气溶胶总量具有很高的正相关性，体现了本发明方法反演气溶胶浓度的可靠性。
[0081]
实施例提供的大气污染物垂直分布的实时反演方法，通过步骤2和步骤3两步骤，只需要使用两次辐射传输模型的模拟光强i即可以求出模拟光路总浓度及迭代权重函数，这大大减少了求解迭代权重函数的计算开销和消耗时长，提高了计算效率，保证了反演的实时性，在此基础上，通过步骤4，采用基于最小二乘法的最优化迭代估计方法迭代求解大气污染物垂直分布，保证了反演污染物浓度垂直分布稳定性。
[0082]
图6是一实施例提供的大气污染物垂直分布的实时反演装置的结构示意图。如图6所示，实施例提供的大气污染物垂直分布的实时反演装置，包括：
[0083]
实测光路总浓度计算模块，用于获取利用超光谱遥感技术实时采集的多仰角观测的高分辨率光谱，根据高分辨率光谱拟合得到大气污染物多仰角观测的实测光路总浓度；
[0084]
模拟光路总浓度计算模块，用于采用辐射传输模型模拟地面光强得到模拟光强，
并结合大气污染物特征吸收截面将模拟光强转换为模拟光路总浓度；
[0085]
迭代权重函数计算模块，用于计算模拟光强对分层的大气污染物浓度分布的偏微分以得到模拟光强的雅克比矩阵，利用模拟光强的雅克比矩阵并结合大气污染物特征吸收截面计算模拟光路总浓度的雅克比矩阵，即迭代权重函数；
[0086]
迭代计算模块，以最小化带迭代权重函数的实测光路总浓度和模拟光路总浓度为目标函数，采用基于最小二乘法的最优化迭代估计方法迭代求解以得到大气污染物垂直分布。
[0087]
需要说明的是，需要说明的是，上述实施例提供的大气污染物垂直分布的实时反演装置在进行大气污染物垂直分布的实时反演时，应以上述各功能模块的划分进行举例说明，可以根据需要将上述功能分配由不同的功能模块完成，即在终端或服务器的内部结构划分成不同的功能模块，以完成以上描述的全部或者部分功能。另外，上述实施例提供的大气污染物垂直分布的实时反演装置与大气污染物垂直分布的实时反演方法实施例属于同一构思，其具体实现过程详见大气污染物垂直分布的实时反演方法实施例，这里不再赘述。
[0088]
实施例还提供一种计算机设备，包括存储器、处理器以及存储在存储器中并可在处理器上执行的计算机程序，所述处理器执行计算机程序时实现上述大气污染物垂直分布的实时反演方法，即包括以下步骤：
[0089]
步骤1，获取利用超光谱遥感技术实时采集的多仰角观测的高分辨率光谱，根据高分辨率光谱拟合得到大气污染物多仰角观测的实测光路总浓度；
[0090]
步骤2，采用辐射传输模型模拟地面光强得到模拟光强，并结合大气污染物特征吸收截面将模拟光强转换为模拟光路总浓度；
[0091]
步骤3，计算模拟光强对分层的大气污染物浓度分布的偏微分以得到模拟光强的雅克比矩阵，利用模拟光强的雅克比矩阵并结合大气污染物特征吸收截面计算模拟光路总浓度的雅克比矩阵，即迭代权重函数；
[0092]
步骤4，以最小化带迭代权重函数的实测光路总浓度和模拟光路总浓度为目标函数，采用基于最小二乘法的最优化迭代估计方法迭代求解以得到大气污染物垂直分布。
[0093]
实际应用中，存储器可以为在近端的易失性存储器，如ram，还可以是非易失性存储器，如rom，flash，软盘，机械硬盘等，还可以是远端的存储云。处理器可以为中央处理器(cpu)、微处理器(mpu)、数字信号处理器(dsp)、或现场可编程门阵列(fpga)，即可以通过这些处理器实现大气污染物垂直分布的实时反演步骤。
[0094]
实施例还提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理执行时实现上述大气污染物垂直分布的实时反演方法，即包括以下步骤：
[0095]
步骤1，获取利用超光谱遥感技术实时采集的多仰角观测的高分辨率光谱，根据高分辨率光谱拟合得到大气污染物多仰角观测的实测光路总浓度；
[0096]
步骤2，采用辐射传输模型模拟地面光强得到模拟光强，并结合大气污染物特征吸收截面将模拟光强转换为模拟光路总浓度；
[0097]
步骤3，计算模拟光强对分层的大气污染物浓度分布的偏微分以得到模拟光强的雅克比矩阵，利用模拟光强的雅克比矩阵并结合大气污染物特征吸收截面计算模拟光路总浓度的雅克比矩阵，即迭代权重函数；
[0098]
步骤4，以最小化带迭代权重函数的实测光路总浓度和模拟光路总浓度为目标函
数，采用基于最小二乘法的最优化迭代估计方法迭代求解以得到大气污染物垂直分布。
[0099]
其中，计算机可读存储介质可以是rom、随机存取存储器(ram)、 cd-rom、磁带、软盘和光数据存储设备等。
[0100]
以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。