基于EMI超光谱卫星载荷的自适应迭代二氧化硫反演方法

基于emi超光谱卫星载荷的自适应迭代二氧化硫反演方法
技术领域
1.本发属于空气质量监测技术领域，具体涉及全球大气中二氧化硫反演算法。


背景技术：

2.大气污染精准监测是污染问题管控和治理的重要前提。大气中二氧化硫(so2)主要有两种来源：一是火山喷发等自然源，二是化石燃烧等人为源。自工业革命以来，以化石燃料燃烧为代表的人类活动成为大气so2浓度增长的主要因素，并且引发酸雨和硫酸盐气溶胶等一系列环境问题。卫星观测通常能够实现大范围全局观测，有着地面监测无法企及的优势。但是针对so2的反演，由于so2的强吸收波段一般位于观测波段边缘，o3的强吸收干扰，宇宙射线辐射、氧原子和紫外暴露等复杂空间条件，往往给so2的反演带来更多的困难。gf-5(高分五号卫星)是我国首颗用于痕量气体监测的卫星,其搭载的emi超光谱载荷可实现多种气体的每日全球化观测。我国有望依托于emi实现全球化so2的天基监测业务，实现该任务的前提是自主开发出适用于emi仪器特性的大气so2反演方法。


技术实现要素：

3.本发明的目的在于，针对so2反演存在的困难，提供一种基于国产emi仪器特性和光谱特征的，能够实现光谱信噪比低和稳定性差情形下的高精度的自适应迭代二氧化硫反演方法。
4.本发明提供的基于国产emi超光谱的自适应迭代二氧化硫反演方法，采用自适应迭代技术，通过对包含so2柱总量在内的多种参数同时迭代拟合，结合大气辐射传输模型和大气消光定律，使观测与模拟的大气反照率数据之间达到最优拟合。具体步骤如下。
5.(一)计算模拟大气反照率：
6.按照大气辐射传输模型计算包含臭氧吸收，氧气和氮气对太阳光的散射和吸收，气溶胶和云对太阳光的散射吸收等因素光谱信息；根据大气消光定律，添加其它痕量气体吸收信息，计算公式为：
[0007][0008]
其中，i
simu
为模拟大气反照率，i
model
为按照大气传输模型计算的大气反照率，w为波长，crosec为气体吸收特征，isrf为高分五号卫星的仪器函数，v为气体的柱总量，a为大气质量因子，igas为计算中添加的痕量吸收气体；表示卷积。
[0009]
(二)获取观测光谱数据波长信息：
[0010]
针对emi未对紫外-可见甚高光谱数据(包括对地球观测光谱数据和对太阳观测光谱数据)直接提供对应波长信息(1473个波长维信息)的情况，对观测数据进行观测波长计算，计算公式为：
[0011]wi,j
＝a
1,j
+(a
2,j
×
(1437-i))+(a
2,j
×
(1437-i)2)，
ꢀꢀ
(2)；
[0012]
其中，a为观测数据中提供的波长因子，w
i,j
为对于第j个列像元第i个波长像元处
的波长。
[0013]
(三)对观测光谱进行预处理：
[0014]
在复杂的太空环境中，仪器性能会发生变化；因此需要对观测光谱进行预处理，预处理包括对在轨仪器函数和波长校正、太阳光谱重构、在轨辐射校正；具体流程为：
[0015]
(1)对在轨仪器函数和波长重定标：
[0016]
针对emi仪器未提供发射前仪器函数数据，选择在轨观测光谱计算在轨仪器函数，并对观测的波长进行重定标：选取相对干净的太平洋区域的对地面观测光谱数据；需要选择晴朗的像元，降低云和气溶胶的吸收散射，对于so2的吸收波段，该地区除了o3以外的其他气体对太阳光的吸收强度低，可以忽略不计；由于卫星在复杂的太空中性能会发生变化，观测的波长与真实波长之间存在偏移，因此需要对观测光谱的波长进行重定标：
[0017][0018]
其中，w
scaling
为波长拉伸量，w
baseline
为波长偏移量，w为波长，为反演波段的波长均值，w
update
为重定标波长；为了避免波长重定标时出现过拟合现象，在拟合的过程选择十个波峰点进行粗定标，需要先通过粗定标判断波长重定标波动的范围；
[0019]
将高分辨率太阳参考谱与仪器响应函数卷积作为模拟的大气光谱，其中需要基于比尔朗伯定律考虑o3的吸收，用低阶多项式考虑o2和n2的散射吸收，得到模拟观测光谱；其中的每日在轨仪器响应函数为对称超高斯函数：
[0020][0021]
其中，a
scal
，ai为拟合参数；i0为高分辨太阳参考谱，crosec
ozone
为臭氧特征吸收结构，s
ozone
为臭氧的斜柱总量，rad
simu
为模拟光谱；
[0022]
重复迭代波长重定标值和仪器响应函数，至模拟与观测光谱数据之间的残差稳定且最小，达到最优拟合，得到每日的在轨仪器函数和波长重定标值。
[0023]
(2)重构太阳光谱：
[0024]
每日的对太阳观测光谱主要用于计算大气反照率，计算仪器函数和波长重定标；其中，观测大气反照率计算式为：
[0025][0026]
其中，i
meas
为观测大气反照率，pi为圆周率，rad和irrad为emi观测的对地球和对太阳光谱数据，solzenang为观测的太阳天定角数据，w为波长；
[0027]
由于emi未提供在轨每日太阳光谱数据，因此在流程(1)中选择干净太平洋区域的光谱对仪器函数和波长偏移进行计算；但是仍然需要在轨每日太阳光谱来计算大气反照率；通过由算式(4)计算的仪器函数和波长重定标值，重构每日太阳观测光谱，具体如下：
[0028]
通过流程(1)中得到的仪器响应函数和波长重定标结果，将高分辨太阳参考谱卷积在轨当日仪器狭缝函数，同时考虑当日的日地距离和波长偏移，重构得到每日的太阳参考谱：
[0029][0030]
其中，irrad
simu
为重构的每日太阳光谱，i0为高分辨率太阳参考谱，r为当日的日地距离。
[0031]
(3)在轨辐射校正：
[0032]
由于so2的强吸收波段主要位于光谱观测边缘，观测光谱存在明显的观测强度偏差；因此需要对观测大气反照率进行在轨辐射定标；并按如下方式计算得到辐射强度校正系数：
[0033]
通过算式(1)对干净太平洋地区的大气反照率进行单次拟合，计算模拟大气反照率：选择在云量和气溶胶较低的像元，在so2的强吸收波段，仅仅在辐射传输模型中输入o3的廓线信息，其中的o3廓线由两部分组成：1)当压强大于215百帕，选择来自mcpeter的气候学臭氧廓线和2)当压强小于215百帕时，选择已经过验证的每日mls臭氧廓线(ml2o3.004)；通过算式(5)，利用流程(2)中计算得到的重构每日太阳参考谱计算观测大气反照率；
[0034]
通过对比单次模拟和观测大气反照率，得到与波长和观测列像元相关的比例系数；该比例系数作为辐射校正函数在后续对观测大气反照率进行辐射校正。
[0035]
(四)对so2进行反演，具体流程为：
[0036]
(1)反演需要考虑的吸收：
[0037]
emi的紫外观测甚高光谱范围为306-403nm；so2的强吸收波段为305-330nm；此波段包含有hcho、bro、o3的吸收干扰，因此在反演参数中需要加入hcho、o3和bro的柱总量在算式(1)重的v
igas
中；同时，云高、云量、气溶胶等参数都会影响光的吸收散射，考虑在大气辐射传输模型中i
model
。
[0038]
(2)获取反演参数初始值：
[0039]
o3的初始值和先验值相同，由两部分组成：1)当压强大于215百帕，选择来自mcpeter的气候学臭氧廓线；2)当压强小于215百帕时，选择已经过验证的每日mls臭氧廓线(ml2o3.004)；其他气体hcho、bro和so2的柱总量先验值来自于geos-chem模型模拟的月平均结果，初始值均为0。对于卫星痕量气体的计算，云高，云量的初始值和先验值一般来自于该卫星的云产品。由于目前未提供有效的emi云产品，因此采用空间分辨率更高的哨兵五号云产品数据进行重新格点，得到与emi像元对应的云产品作为云高和云量的初始值。气溶胶的初始值和先验值来自于geos-chem模型模拟的气溶胶吸收强度。
[0040]
(3)选择气体吸收截面：
[0041]
其中，hcho的吸收截面选择的是298k下meller and noortgat观测数据，bro的吸收截面选择的是wilmouth等人观测的228k下的数据。so2的吸收截面为由203,223,243,273,293k5个温度下bogumil等人观测的so2吸收截面考虑温度廓线组合而成。o3的吸收截面为由218,228,243,295k4个温度下brion等人观测的o3吸收截面组合而成。
[0042]
(4)获取每日的气象数据：
[0043]
包括：每日地表压强数据，每日地标温度数据，每日温度廓线，每日对流层顶压强数据。数据来自于每日ncep在分析数据。
[0044]
(5)筛选so2的反演窗口：
[0045]
emi的紫外观测甚高光谱范围为306～403nm。so2的强吸收波段为305～330nm。考虑到不同仪器的观测光谱特性不同，需要针对emi的观测光谱对so2的反演波段进行筛选。选择的起始波段为307～315nm，终止波段为320～330nm。步长为0.5nm。其中选择的反演区域为干净的太平洋地区，在无火山喷发的情况下，该区域的so2均值接近与0du。针对不同的筛选反演窗口，主要分析其中的so2的反演柱总量接近与0，同时考虑到反演不确定度和反
演拟合残差尽量小，反演成功率尽量高，反演rms尽量接近与1。最后得到的最优反演窗口为310.5～320.0nm。
[0046]
(6)自适应多次迭代反演so2柱总量：
[0047]
将流程(1)-(5)中的参数设置输入到算式(1)中，得到模拟的大气反照率，和算式(5)计算的观测的大气反照率进行对比，通过多次自适应多次迭代算式(7)直至算式(8)观测与模拟的大气反照率和先验与反演的反演参数的区别达到最小且稳定。
[0048][0049][0050]
其中，retri表示第i次迭代下的反演结果，其中包含了so2的柱总量。i
meas
和i
simu
分别表示观测和模拟的大气反照率，retr
apri
为反演参数的先验值，err
apri
为反演参数先验值的误差，来自于先验模型输出。err
meas
为观测光谱的误差。avgker为反演参数的权重函数，
[0051]
(五)二级产品的校正：
[0052]
针对emi搭载的二维ccd性能参数在行像元上存在不一致性，对反演的so2进行背景校正，包括：沿着行像元的去条带校正和沿着纬度方向的校正。得到校正以后的不规则so2垂直柱总量。
[0053]
(六)格点得到三级产品：
[0054]
该方法还包括：通过对像元18km周围作平均，将不规则的二级so2垂直柱总量进行重新格点化，得到全球区域的规则网格化三级so2柱总量产品。
[0055]
本方法针对我国首颗用于痕量气体监测的星载光谱仪emi在紫外波段信噪比较低，辐射定标偏差严重，每日对太阳观测光谱缺失的情况下，对光谱进行预处理，通过对反演参数设置的筛选，提出了一种大气so2柱总量的反演方法，与未进行校正的方法相比，反演结果残差大大降低，反演结果更符合我国so2排放逐年降低的现象。
附图说明
[0056]
图1为emi紫外波段仪器函数半高宽。
[0057]
图2为不同气体的吸收截面。
[0058]
图3为emi so2波段测试结果对比。
[0059]
图4为emi 2020年3月月平均so2柱总量和拟合残差的分布。
具体实施方式
[0060]
下面通过具体实施例进一步介绍本发明方法。
[0061]
实施例：
[0062]
计算2019年8月so2柱总量：
[0063]
获取观测光谱的波长信息：由于观测光谱中未提供波长信息，因此需要对每个波长像元对应的波长进行计算(1,2,3
…
1473个波长像元)。利用在一级光谱文件中提供了波长因子(3*191)和列像元位置(1，2，3
…
191)计算得到观测光强对应的波长信息。
[0064]
观测光谱预处理：
[0065]
(1)由于gf-5未提供仪器响应函数，因此需要计算在轨仪器响应函数。由于仪器响应函数随着波长变化，因此在此算法中选取用于计算仪器响应函数的波长涵盖so2的强吸收波段(305～340nm)。计算步骤为：选取每日的太平洋地区的光谱平均得到观测光谱，以降低光谱的观测误差。同时考虑到仪器在复杂太空中可能存在的性能变化，光谱的波长定标可能存在飘移，在拟合的同时对波长进行重定标：
[0066][0067]wscaling
为波长拉伸量，w
baseline
为波长偏移量，w为波长，为反演波段的波长均值，w
update
为重定标波长。其中，为了避免波长重定标时出现过拟合现象，在拟合的过程选择十个波峰点进行粗定标，需要先通过粗定标判断波长重定标波动的范围。
[0068]
将高分辨率的太阳谱与仪器函数进行卷积得到低分辨率的光谱，同时考虑o3吸收和其他气体的散射吸收，得到模拟光谱：
[0069][0070]
其中a
scal
，ai为拟合参数。i0为高分辨太阳参考谱，crosec
ozone
为臭氧特征吸收结构，s
ozone
为臭氧的斜柱总量，rad
simu
为模拟光谱,表示卷积。
[0071]
多次迭代拟合观测和模拟光谱至拟合残差最小，得到在轨计算的仪器函数(图1)。
[0072]
(2)针对观测光谱中未提供每日的对太阳观测光谱的问题，需要重构每日对太阳观测光谱。通过以上计算出来的仪器响应函数与高分辨率的太阳谱卷积，同时考虑到波长重定标波动值，重构当日观测太阳光谱：
[0073][0074]
其中，irrad
simu
为重构的每日太阳光谱，i0为高分辨率太阳参考谱，r为当日的日地距离。
[0075]
(3)针对仪器观测边缘的光谱存在变形等光强偏差等问题，对观测大气反照率进行在轨辐射定标。对比首次模拟得到的大气反照率和观测大气反照率，得到与波长和列像元相关的辐射强度校正系数。应用于后续观测光谱中。
[0076]
模拟大气反照率：按照大气辐射传输模型计算包含臭氧吸收，氧气和氮气对太阳光的散射和吸收，气溶胶和云对太阳光的散射吸收等因素光谱信息；根据大气消光定律，添加其它痕量气体吸收信息，计算公式为：
[0077][0078]
其中，i
simu
为模拟大气反照率，i
model
为按照大气传输模型计算的大气反照率，w为波长，crosec为气体吸收特征，isrf为emi的仪器函数，v为气体的柱总量，a为大气质量因子，igas为计算中添加的痕量吸收气体；表示卷积。
[0079]
对于首次模拟，仅仅考虑o3的吸收，其中o3的廓线信息由两部分组成：1)当压强大于215百帕，选择来自mcpeter的气候学臭氧廓线；2)当压强小于215百帕时，选择已经过验证的每日mls臭氧廓线(ml2o3.004)。
[0080]
利用(2)中计算得到的重构每日太阳谱和观测对地光谱计算观测大气反照率：
[0081][0082]
其中，i
meas
为观测大气反照率，pi为圆周率，rad和irrad为emi观测的对地球和对太阳光谱数据，solzenang为观测的太阳天顶角数据，w为波长。
[0083]
考虑so2强吸收波段可能存在的各种吸收：对比在so2强吸收波段的各种气体吸收截面，其中hcho，bro，o3在该波段有较强的吸收(图2)，因此在反演参数中需要添加这些气体的吸收截面。这些气体的吸收截面选择为：hcho的吸收截面选择的是298k下meller and noortgat观测数据，bro的吸收截面选择的是wilmouth等人观测的228k下的数据。so2和o3的吸收截面考虑了温度的影响。将多个温度的吸收截面用多项式拟合得到随温度变化的吸收截面函数。其中，so2的吸收截面为由203,223,243,273,293k5个温度下bogumil等人观测的so2吸收截面考虑温度廓线组合而成。o3的吸收截面为由218,228,243,295k4个温度下brion等人观测的o3吸收截面组合而成。
[0084]
反演设置的初始值来自于模型或者观测数据：o3的初始值和先验值相同，由两部分组成：1)当压强大于215百帕，选择来自mcpeter的气候学臭氧廓线；2)当压强小于215百帕时，选择已经过验证的每日mls臭氧廓线(ml2o3.004)；其他气体hcho、bro和so2的柱总量先验值来自于geos-chem模型模拟的月平均结果，初始值均为0。对于卫星痕量气体的计算，云高，云量的初始值和先验值一般来自于该卫星的云产品。由于目前未提供有效的emi云产品，因此采用空间分辨率更高的哨兵五号云产品数据进行重新格点，得到与emi像元对应的云产品作为云高和云量的初始值。气溶胶的初始值和先验值来自于geos-chem模型模拟的气溶胶吸收强度。
[0085]
每日的气象数据包括：每日地表压强数据，每日地标温度数据，每日温度廓线，每日对流层顶压强数据。从ncep在分析数据网站上下载获取：
[0086]
(https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html)。
[0087]
对于so2和o3的吸收截面，根据每日温度廓线进行吸收截面的计算。
[0088]
so2反演波段的选择：由于so2在大气中浓度低，并且该波段的o3吸收明显。因此对so2反演窗口针对不同的仪器性能存在区别，并且对so2的反演至关重要。通过权力要求11的方法，最终计算的最优反演窗口为310.5～320.0nm(图3)。
[0089]
将以上参数设置输入到算式(4)中，得到模拟的大气反照率，和算式(5)计算的观测的大气反照率进行对比，通过多次自适应多次迭代算式(6)直至算式(7)观测与模拟的大气反照率和先验与反演的反演参数的区别达到最小且稳定。
[0090]
[0091][0092]
其中，retri表示第i次迭代下的反演结果，其中包含了so2的柱总量。i
meas
和i
simu
分别表示观测和模拟的大气反照率，retr
apri
为反演参数的先验值，err
apri
为反演参数先验值的误差，来自于先验模型输出。err
meas
为观测光谱的误差。avgker为反演参数的权重函数，
[0093]
由于gf-5的探测器为二维仪器，不同列像元之间存在仪器区别，导致反演结果存在条带现象，因此需要对结果进行沿着行像元的去条带校正和沿着纬度方向的校正。选取干净的太平洋区域，该区域的理论so2柱总量为0，计算该区域的观测so2柱总量与理论值沿着行像元和纬度向的波动作为校正值。校正以后得到每日观测的二级不规则so2柱总量产品。
[0094]
在将二级产品进行重格点为规则像元。首先对所选择的区域进行格点化，选择的分辨率为0.1度*0.1度，以规则格点中心为圆心，半径18km以内的所有二级产品像元进行平均，作为该规则格点处的so2柱总量，得到规则的三级产品(图4(c)和(d))。相对于没有经过预处理的反演结果，反演的拟合残差大幅度降低(图4(a)和(c))，反演结果更符合我国so2排放逐年降低的事实(图4(b)和(d))。