一种基于MAX-DOAS的50米分辨率痕量气体廓线反演方法

一种基于max
‑
doas的50米分辨率痕量气体廓线反演方法
技术领域
1.本发明涉及环境监测技术领域，特别涉及一种基于max
‑
doas的50米分辨率痕量气体廓线反演方法。


背景技术：

2.随着城市化进程，大气污染形势日益严峻，大气污染物的实时监测在污染物的治理过程中显得尤为重要。max
‑
doas(多轴差分吸收光谱)作为地基遥测的一种重要手段，可以连续自动工作并提供污染物(气溶胶、二氧化硫、二氧化氮、甲醛、气态亚硝酸等)的立体分布廓线，被广泛使用。
3.但是目前的基于max
‑
doas痕量廓线反演方法仍然具有一些明显的缺陷：
4.1、现有max
‑
doas反演算法的廓线立体分辨率均为200米，意味着对于气体的高值立体分布状况并不能清晰定位，尤其是对于一些近地面分布的气体(例如气态亚硝酸等)并不能显示真正的立体分布情况。同时由于立体分辨率太高，反演误差中平滑误差占比通常超过70％，极大降低了算法的精确度。
5.2、现有的max
‑
doas气体廓线反演算法都是为了满足价值函数最小而设计的，无一例外的采用了基于高斯牛顿方法迭代求解。但是迭代求解的方法会有三点重要缺陷：
6.1)效率低；常用的迭代求解算法(例如priam等)使用最大迭代20次求解廓线，这极大增加了矩阵运算量，使得算法的求解效率很低，消耗太多求解时长。
7.2)出现奇异值；算法的迭代过程采用统一的迭代步长，使得高层(1km以上)的反演通常出现奇异值，这限制了迭代算法在高分辨率上的应用。
8.3)精度差；迭代过程采用阈值判定法，阈值的采用有两个：最大迭代次数阈值和最小价值函数阈值，这意味着迭代求解的廓线是局部最优解，并不能获得最精确的全局最优解；同时，最大迭代次数也导致了奇异值的出现。


技术实现要素：

9.为解决上述问题，本发明旨在提出一种基于max
‑
doas的50米分辨率痕量气体廓线反演方法，通过蒙特卡罗方法第一步反演气溶胶消光系数廓线，第二步通过气溶胶消光系数廓线反演痕量气体垂直廓线，通过蒙特卡罗方法进行两步反演法极大程度的降低了反演的误差，提高了算法的可靠性和精度。
10.为达到上述目的，本发明的技术方案是这样实现的：
11.一种基于max
‑
doas的50米分辨率痕量气体廓线反演方法，所述反演方法包括以下步骤：
12.s1、输入原始光谱，通过光谱拟合获得o4斜柱浓度以及气体斜柱浓度；
13.s2、所述o4斜柱浓度结合气溶胶先验廓线通过蒙特卡罗采样算法第一次反演计算得到气溶胶廓线；
14.s3、所述气溶胶廓线结合气体斜柱浓度、气体先验廓线通过蒙特卡罗采样算法第
二次反演计算得到痕量气体廓线。
15.进一步的，所述蒙特卡罗采样算法的主体流程为：
16.输入斜柱浓度和先验廓线到辐射传输模型中，产生气体的垂直柱浓度与各层的权重函数；
17.接下来将权重函数和先验廓线进行归一化处理，输入到蒙特卡罗随机采样计算流程中，计算得到最佳线型；
18.最佳线型与垂直柱浓度点乘，获得立体廓线。
19.进一步的，所述斜柱浓度为o4斜柱浓度或气体斜柱浓度。
20.进一步的，所述先验廓线为气溶胶先验廓线或气体先验廓线。
21.进一步的，所述蒙特卡罗随机采样计算流程为：
22.首先按照输入的最大采样次数产生0和1之间的随机采样序列，该序列为权重因子；
23.结合权重因子表征不同的随机采样线型，记录该线型下的价值函数，并储存下来；
24.当采样过程结束后，采用排序方法搜索采样过程中的最小价值函数，输出对应的线形和权重因子，即最佳线型。
25.进一步的，所述蒙特卡罗采样算法的两次反演计算过程为：
26.均使用最小化价值函数：
[0027][0028]
其中，y为max
‑
doas观测值，f为辐射传输模型模拟值，ε
y
表示观测值y的误差，x
p
表示大气成分的先验廓线，ε
x
表示先验廓线的误差项，其中m，n分别表示网格化模型的高度格点层数；
[0029]
上述反演需要使用到辐射传输模型来模拟气体在大气中的真实立体分布状态，大气辐射传输模型模拟值f受到观测、气象、观测误差多种影响，计算时简化为如下模型：
[0030]
f(y)＝k(y，b)+ε
ꢀꢀꢀꢀ
(2)；
[0031]
上式中k表示各仰角下的权重函数，y，b，ε分别表示观测、气象、观测误差三项影响因素；权重函数k具体形式可表示模型对于各分量的偏导数为：
[0032][0033]
其中i表示仰角为i时的模拟值，j表示高度格点层数；
[0034]
采用蒙特卡罗随机抽样模型求解(1)式，采用权重因子α衡量模型和先验之间的关系，解的表现形式为：
[0035][0036]
[0037][0038]
其中k和vcd为辐射传输模型模拟的输出，这两项综合了观测与模拟的结果，x
p
为先验廓线，box表示盒子模型的高度格点数量，ea表示观测仰角数量，||k||和||x
p
||分别表示权重函数与先验廓线的一阶范数，通过归一化的除式构造可以获得对应的线型，对于相同的价值函数处的值x
α
与x，满足：
[0039]
||x
α
‑
x||
→0ꢀꢀꢀꢀꢀ
(7)；
[0040]
该算法的误差可分为先验误差s
p
、测量误差s
m
和模型误差s
k
，总误差s
total
服从误差传递公式，其中先验误差s
p
与模型误差s
k
，通过权重因子α计算：
[0041][0042][0043]
总误差s
total
服从误差传递公式，可表示为：
[0044][0045]
有益效果：本发明首次实现了50m分辨率大气痕量气体廓线反演，将200m分辨率提升至50m，极大程度降低了反演的误差。本发明采用了全新的更高效率的求解算法，通过采用蒙特卡罗采样算法求解，避免了传统的牛顿迭代法造成的多点问题，同时极大增加了算法求解的效率，提高求解的精度与可靠性。
附图说明
[0046]
构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
[0047]
图1为本发明实施例所述的基于max
‑
doas的50米分辨率痕量气体廓线反演方法的流程图；
[0048]
图2为本发明实施例所述的基于max
‑
doas的50米分辨率痕量气体廓线反演方法的蒙特卡罗采样算法的主体流程图；
[0049]
图3为本发明实施例所述的基于max
‑
doas的50米分辨率痕量气体廓线反演方法的蒙特卡罗随机采样计算流程图。
具体实施方式
[0050]
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
[0051]
下面将参考附图并结合实施例来详细说明本发明。
[0052]
实施例1
[0053]
参见图1：本实施例的一种基于max
‑
doas的50米分辨率痕量气体廓线反演方法，所述反演方法包括以下步骤：
[0054]
s1、输入原始光谱，通过光谱拟合获得o4斜柱浓度以及气体斜柱浓度；
[0055]
s2、所述o4斜柱浓度结合气溶胶先验廓线通过蒙特卡罗采样算法第一次反演计算得到气溶胶廓线；
[0056]
本实施例的第一次反演气溶胶的状态采用分辨率为200m。
[0057]
s3、所述气溶胶廓线结合气体斜柱浓度、气体先验廓线通过蒙特卡罗采样算法第二次反演计算得到痕量气体廓线。
[0058]
本实施例的第二次用于反演痕量气体的垂直廓线，采用50m分辨率。
[0059]
综上所述，本实施例实现了痕量气体廓线分辨率由200m提升为50m，这极大降低了反演算法的平滑误差，提升了立体分辨率。另外，本实施例采用蒙特卡罗方法替代高斯牛顿迭代求解，蒙特卡罗求解过程通过调整采样次数，使得算法的解无限接近全局最优解，极大提高了算法的可靠性、精度和求解速率，提升了算法的效率。
[0060]
参见图2：在一具体的实例中，所述蒙特卡罗采样算法的主体流程为：
[0061]
输入斜柱浓度和先验廓线到辐射传输模型中，产生气体的垂直柱浓度与各层的权重函数；
[0062]
接下来将权重函数和先验廓线进行归一化处理，输入到蒙特卡罗随机采样计算流程中，计算得到最佳线型；
[0063]
最佳线型与垂直柱浓度点乘，获得立体廓线。
[0064]
需要说明的是，各种max
‑
doas廓线反演算法都是采用迭代法来逼近最优解，本实施例创新地提出解的表现形式，直接表现出解的两大影响因素：模型与先验廓线，其中模型对垂直柱浓度和权重产生影响，先验廓线仅仅对线型产生影响。本实施例引进线型的概念与表征，通过线型这一概念，用来表述气体立体分布的状态，将垂直柱浓度与线型分开处理，使得求解过程更加清晰简洁。
[0065]
另外，本实施例的归一化处理方法，采用一阶范数归一化处理先验廓线与权重函数，有效地防止了数字大小的误差对于求解过程的影响，解决了先验廓线不精确时反演出现奇异值的情况。
[0066]
在一具体的实例中，所述斜柱浓度为o4斜柱浓度或气体斜柱浓度，所述先验廓线为气溶胶先验廓线或气体先验廓线。
[0067]
参见图3：在一具体的实例中，所述蒙特卡罗随机采样计算流程为：
[0068]
首先按照输入的最大采样次数产生0和1之间的随机采样序列，该序列为权重因子；
[0069]
结合权重因子表征不同的随机采样线型，记录该线型下的价值函数，并储存下来；
[0070]
当采样过程结束后，采用排序方法搜索采样过程中的最小价值函数，输出对应的线形和权重因子，即最佳线型。
[0071]
可以理解的是，本实施例采用蒙塔卡罗方法求解最优权重因子来表现线型，蒙塔卡罗方法是采用随机抽样的方法来求解数学模型，这一采样过程与计算机的计算过程高度相关，所以计算机的求解速度大为提升，对cpu的利用率也有效提升，使得廓线的求解速度得到很大提高。
[0072]
本实施例的权重因子采样区间的设置，将权重因子限制在0和1之间，不仅使得解的形式具有实际物理意义，更有效避免了权重因子数值溢出对于反演的影响，使反演算法
必定收敛，不存在奇异值的情况。
[0073]
在一具体的实例中，所述蒙特卡罗采样算法的两次反演计算过程为：
[0074]
均使用最小化价值函数：
[0075][0076]
其中，y为max
‑
doas观测值，f为辐射传输模型模拟值，ε
y
表示观测值y的误差，x
p
表示大气成分的先验廓线，ε
x
表示先验廓线的误差项，其中m，n分别表示网格化模型的高度格点层数；
[0077]
上述反演需要使用到辐射传输模型来模拟气体在大气中的真实立体分布状态，大气辐射传输模型模拟值f受到观测、气象、观测误差多种影响，计算时简化为如下模型：
[0078]
f(y)＝k(y，b)+ε(2)；
[0079]
上式中k表示各仰角下的权重函数，y，b，ε分别表示观测、气象、观测误差三项影响因素；权重函数k具体形式可表示模型对于各分量的偏导数为：
[0080][0081]
其中i表示仰角为i时的模拟值，j表示高度格点层数；
[0082]
采用蒙特卡罗随机抽样模型求解(1)式，采用权重因子α衡量模型和先验之间的关系，解的表现形式为：
[0083][0084][0085][0086]
其中k和vcd为辐射传输模型模拟的输出，这两项综合了观测与模拟的结果，x
p
为先验廓线，box表示盒子模型的高度格点数量，ea表示观测仰角数量，||k||和||x
p
||分别表示权重函数与先验廓线的一阶范数，通过归一化的除式构造可以获得对应的线型，对于相同的价值函数处的值x
a
与x，满足：
[0087]
||x
α
‑
x||
→0ꢀꢀꢀꢀ
(7)；
[0088]
该算法的误差可分为先验误差s
p
、测量误差s
m
和模型误差s
k
，总误差s
total
服从误差传递公式，其中先验误差s
p
与模型误差s
k
，通过权重因子α计算：
[0089][0090][0091]
总误差s
total
服从误差传递公式，可表示为：
[0092][0093]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。