本发明涉及卫星遥感技术领域，具体涉及一种deei的基于卫星观测的大气二氧化碳浓度反演算法。
背景技术：
二氧化碳是重要的温室气体之一，从工业革命(1750年)开始大气二氧化碳浓度持续上升，noaa(nationaloceanicandatmosphericadministration)在美国夏威夷冒纳罗亚(maunaloa,hawaii)观象台数据显示，二氧化碳浓度于2018年已经超过410ppm(partspermillion)，近十年全球大气二氧化碳浓度年均增长量已经达到2ppm，其中2015年增长量超过了3ppm。最新的政府间气候变化专门委员会(ipcc)评估报告指出，人为排放的温室气体是大气温室效应的主要原因，人类活动很大程度上影响全球气候变化，进而导致极端气候事件增多。目前针对减少温室气体排放的议题已成为国际上密切关注与重点讨论的内容，掌握全球范围大气二氧化碳浓度分布，有助于厘清各个国家碳源汇现状，预测气候变化趋势，具有非常重要的科学意义与社会价值。
大气二氧化碳浓度可以通过卫星遥感进行观测，原理为利用探测器接收地表反射的太阳光辐射光谱信息进行xco2(thecolumn-averageddry-airmolefractionsofcarbondioxide)反演，在xco2的反演算法中，主要分为经验算法和物理算法两大类。其中经验算法是对大量的观测样本进行训练，包括统计回归法和神经网络法两种方法；全物理反演算法，是模拟辐射传输的物理过程，全面考虑大气和地表参数对光谱的影响，例如云和气溶胶的散射，并在反演过程中同时订正这些参数，以获取更高精度的模拟结果。
在反演方法中主要包括正向模型和反向求解算法，其中正向模型负责对卫星遥感观测仪器过程进行数值模拟，模拟所需的参数包括太阳光谱、大气物理化学性质、地表特征、卫星仪器性能等；反向算法可以基于正向模型利用卫星仪器的探测光谱对大气参数进行求解，迭代更新大气廓线状态参量，求出最优解即为反演结果。在反演方法中，涉及两项重要理论基础，大气辐射传输模型与反向求解算法。目前国际上已有很多成熟的大气辐射传输软件被应用于正向模拟实际问题中，比较主流的lblrtm、6s、modtran以及在gomtran 基础上研究的sciatran等，各类反演算法都是选择一个成熟稳定的大气辐射传输模型作为其正向模型，以模拟光谱在大气中辐射传输过程。反向计算一般为非线性问题，需要通过不断迭代得到最优解，目前应用最广泛的是最优估计法eom(optimalestimationmethod)，国际上近红外波段xco2反演方法大多通过最优估计法为基本思想进行求解。
目前国际上有三颗二氧化碳监测卫星，包括gosat、oco-2、tansat，gosat卫星于2009年由日本发射，其官方产品算法为nies算法；oco-2卫星于2014年由美国发射，其官方产品算法为acos算法；tansat卫星于2016年由中国发射，目前虽然有相关反演算法研究，但还未公开l2数据产品。针对不同卫星观测数据，各种算法在数据处理方面各不相同，在模型正向模拟时的参数方案也各有千秋，而不同参数的选取会直接导致模型演变结果的精度和准确性，目前来说，确定有效参数并用于tansat卫星的反演过程以使其演变结果可靠是亟需解决的问题。
技术实现要素：
本发明的目的是为了克服现有技术的不足，提供一种deei的基于卫星观测的大气二氧化碳浓度反演算法，能够至少解决现有技术中不同卫星之间反演参数不能通用、参数选取不当导致反演结果不准确等问题。
为了实现上述目的，本发明提供如下技术方案：
一种deei的基于卫星观测的大气二氧化碳浓度反演算法，所述卫星为二氧化碳监测卫星，所述方法包括以下步骤：
步骤1，提取卫星数据中的卫星参数和观测足迹参数；对云和气溶胶数据分别进行有云区域数据剔除和气溶胶高值区域剔除，得到云和气溶胶参数；
步骤2，设置大气辐射传输模型的参数，将所述卫星参数、所述观测足迹参数、所述云和气溶胶参数以及大气廓线参数输入所述大气辐射传输模型进行计算得到模拟光谱；
步骤3，对所述模拟光谱和所述卫星参数中的观测光谱进行最优估计法求解，并判断所述最优估计法求解的结果是否收敛，如果收敛则得到xco2反演结果，否则更新所述大气辐射传输模型中的参数并重复步骤2和步骤3直至所述最优估计法求解的结果收敛。
基于上述，步骤1中，所述卫星参数包括观测角度、观测时间、仪器信息、光谱信噪比和观测光谱。
基于上述，步骤1中，所述观测足迹参数包括经纬度、地表反照率、地表高程、观测质量标识以及陆地与水域标识。
基于上述，步骤1中，进行有云区域数据剔除和气溶胶高值区域剔除过程中，剔除条件分别为云检测结果中表征有云的数据和气溶胶光学厚度大于0.3的数据。
基于上述，步骤2中，所述大气辐射传输模型中的参数包括观测模式、散射模式以及痕量气体廓线。
基于上述，步骤2中，所述大气廓线参数包括co2廓线、温度廓线、位势廓线、湿度廓线以及压强廓线。
基于上述，步骤2中，计算得到模拟光谱时还得到权重函数信息，所述权重函数信息用于所述最优估计法求解。
基于上述，所述痕量气体廓线为二氧化碳廓线。
基于上述，所述卫星为tansat卫星。
与最接近的现有技术相比，本发明提供的技术方案具有如下优异效果：
本发明公开了一种deei的基于卫星观测的大气二氧化碳浓度反演算法，提取卫星数据中的卫星参数和观测参数；对云和气溶胶数据分别进行有云区域数据剔除和气溶胶高值区域剔除，得到云和气溶胶参数；设置大气辐射传输模型的反演参数，将卫星参数、观测参数、云和气溶胶参数以及大气廓线参数输入大气辐射传输模型进行计算得到模拟光谱；对模拟光谱和观测光谱进行最优估计法求解，并判断最优估计法求解的结果是否收敛，如果收敛则得到xco2反演结果，否则更新大气辐射传输模型中的参数直至最优估计法求解的结果收敛。本发明在现有技术的基础上，针对tansat特性进行各类辅助数据的选取及数据处理，通过与观测位置大气廓线及地表特征等信息进行反演，计算得到大气二氧化碳在空气中平均柱浓度，计算结果精度良好，为地表碳源汇研究提供数据支持。
附图说明
图1为本发明实例中大气二氧化碳浓度反演方法流程图；
图2为本发明实例中tansat提取参数流程图；
图3为本发明实例中北京地区tansat观测反演结果分布图；
图4为本发明实例中北京2018年3月8日tansat观测反演结果分布图；
图5为本发明实例中北京2018年4月9日tansat观测反演结果分布图；
图6为本发明实例中北京2018年4月16日tansat观测反演结果分布图；
图7为本发明实例中北京2018年5月4日tansat观测反演结果分布图；
图8为本发明实例中北京2018年5月24日tansat观测反演结果分布图；
图9为本发明实例中北京2018年5月31日tansat观测反演结果分布图；
图10为本发明实例中北京2018年6月20日tansat观测反演结果分布图；
图11为本发明实例中北京2018年8月20日tansat观测反演结果分布图；
图12为本发明实例中北京2018年11月21日tansat观测反演结果分布图；
图13为本发明实例中北京2018年12月04日tansat观测反演结果分布图；
图14为本发明实例中迭代过程示意图；
图15为本发明实例中温度廓线图；
图16为本发明实例中压强廓线图；
图17为本发明实例中湿度廓线图；
图18为本发明实例中位势廓线图；
图19为本发明实例中xco2模拟分布图。
具体实施方式
下面将参考附图并结合实施例来详细说明本发明。需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
在本发明的描述中，术语“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明而不是要求本发明必须以特定的方位构造和操作，因此不能理解为对本发明的限制。本发明中使用的术语“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接；可以是直接相连，也可以通过中间部件间接相连，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。
本实施例提供了一种deei的基于卫星观测的大气二氧化碳浓度反演算法，属于大气遥感与空间信息模拟领域，其中deei为中国科学院地理科学与资源研究所生态环境信息学研究室(departmentofecologyandenvironmentinformation，deei)。本实施例的反演方法能够从卫星观测的光谱数据，通过与观测位置大气廓线及地表特征等信息进行反演，计算得到大气二氧化碳在空气中平均柱浓度，为地表碳源汇研究提供数据支持，如图1所示，主要包括：
步骤1，数据预处理，包括从碳卫星l1b数据中提取卫星参数和观测足迹参数；对云和气溶胶数据分别进行有云区域数据剔除和气溶胶高值区域剔除，得到云和气溶胶参数；
步骤2，卫星参数和观测足迹参数作为模型观测参数；云和气溶胶参数作为模型云和气溶胶参数，大气廓线参数作为模型大气廓线参数，将模型观测参数、模型云和气溶胶参数、模型大气廓线参数输入至大气辐射传输模型进行计算得到模拟光谱与权重函数信息；
步骤3，对模拟光谱和卫星参数中的观测光谱进行最优估计法求解，并判断最优估计法求解的结果是否收敛，如果收敛则得到xco2反演结果，否则更新大气辐射传输模型中的参数，并重复步骤2和步骤3直至最优估计法求解的结果收敛。
其中卫星参数包括观测角度、观测时间、仪器信息、光谱信噪比和观测光谱等；观测足迹参数包括经纬度、地表反照率、地表高程、观测质量标识以及陆地与水域标识等；云和气溶胶参数包括云检测结果、气溶胶类型、气溶胶光学厚度等；大气廓线参数包括co2廓线、压强廓线、温度廓线、湿度廓线、位势廓线等。卫星观测光谱是实测光谱，将其与模型正向模拟结果作对比看是否结果收敛。
同时进行有云区域数据剔除和气溶胶高值区域剔除过程中，剔除条件分别为云检测结果中表征有云的数据和气溶胶光学厚度大于0.3的数据。
大气辐射传输模型中的参数包括观测模式、散射模式以及痕量气体廓线等模型相关参数，痕量气体包括许多种气体，本实施例中只针对二氧化碳廓线；本实施例中反演时更新二氧化碳廓线值，至其收敛。同时如图1所示，大气辐射传输模型中的太阳模型、散射模型、仪器模型和分子吸收谱线并不是指具体的模型，而是对大气辐射传输模型中各个部分的概括，其中太阳模型对应太阳光谱函数、散射模型对应的是大气传输过程中遇到的散射路径问题、仪器模型对应ils仪器线性模型、分子吸收谱线是指二氧化碳分子的光谱吸收特征；这些内容在大气辐射传输模型中均有相应的规定。
模型观测参数除了“观测光谱信息”，包含所有的卫星参数和观测参数，观测光谱信息作为卫星观测光谱用于与模型正向模拟结果的迭代求解二氧化碳廓线。即更新的“状态变量”就是模型中大气二氧化碳(浓度)廓线。
本实施例针对tansat卫星观测l1b光谱数据，进行大气二氧化碳空气平均柱浓度(xco2)数据反演，选取最适合tansat卫星的反演参数方案与输入数据，基于最优估计法(oem)计算获取数据产品。
本实施例中方法反演tansat卫星，云和气溶胶数据根据与tansat卫星同轨道的风云系列卫星得到，其中由风云四号(fy-4a)获取云检测数据，由风云三号(fy-3d)获取气溶胶数据。
针对tansat研究现状，本实施例方法基于最优估计法和sciatran大气辐射传输模型，通过调整优化适合tansat卫星的提取参数与控制方案，进行大气xco2反演，技术流程包含数据预处理、数据反演、数据优化三个过程。
在数据预处理过程中，利用静止轨道气象卫星，如风云四号(或称fy-4a)的云覆盖产品进行有云区域数据剔除，在此基础上，利用极轨气象卫星，如风云三号(或称fy-3d)的气溶胶产品进行气溶胶高值区域剔除，筛选阈值设置为气溶胶光学厚度大于0.3。在数据预处理方面，使用中国风云系列卫星，与tansat可以保持时空分辨率高度一致，精确的时空匹配数据为下一步反演结果的精度提供了保障。
在数据反演过程中，设置sciatran模型卫星反演参数，包括观测模式、散射模式、痕量气体廓线等，输入观测位置对应的大气压强、温度、湿度、位势廓线数据与co2廓线先验值，此处使用geos-chem大气模式模拟结果，在纬度上以每10°带将全球划分为18个纬度带进行统计，将xco2均值作为各层相同的观测位置大气co2廓线进行输入。将观测信息数据输入至反演参数设置后的sciaran模型中，将代码程序在linux环境下，多节点并行计算处理，得到初步xco2反演结果。针对观测误差，利用数理统计中测量结果符合正态分布原则，剔除观测区域中误差大于三倍标准差的观测点，优化最终反演结果。
本实施例中具体tansat参数提取过程如图2所示，针对tansatl1bv2.0版本数据，各个属性字段提取信息，包括卫星参数和观测足迹信息等，按模型输入要求提取出相关信息，经过处理作为参数输入至模型中进行模拟。图2体现了如何提取tansat数据中的信息用于二氧化碳反演，图2中右侧最终的几个输出会作为sciaran模型输入的一部分参数，例如包括仪器，光谱，信噪比等。
以本实施例方法反演tansat2018年在北京地区目标模式观测数据，反演结果显示9个传感器观测结果稳定，反演结果标准差约0.86ppm，精度良好，如图3至图13所示，分别为北京地区不同日期的tansat观测反演结果分布图，表1为北京地区不同日期的tansat观测反演结果统计表。
表1北京tansat观测反演结果统计
本实施例中涉及到的具体算法包括：
大气xco2反演算法，在模型反演中，利用大气辐射传输模型sciatran(成熟开源模型)，输入这些卫星、地表、仪器参数后，得到模拟光谱，依据选择的光谱带及通道，将模拟光谱与实测光谱代入反演算法进行迭代求解。
(1)反演算法基础
大气co2反演采用rodgers提出的最优估计法，该方法不直接比较前向模型计算辐射值和卫星实测值的差，而是将其与先验猜值廓线值和反演廓线值之间的差值一起定义代价函数，利用先验条件和牛顿迭代方式逐步逼近真解的最大似然估计。
一般反演问题可以理解为构建并求解一系列同时为线性或非线性等式的问题，假设待反演的大气状态可由如下向量形式表示：
x＝(x1,x2,…,xn)(1)
x是待反演的状态向量，下标n代表不同的大气状态参数个数。
y＝(y1,y2,…,ym)(2)
y是测量向量，为m个离散测量值的几何，则卫星观测辐亮度克表示为：
y＝f(x,b) ε(3)
式中，f即前向模型，描述大气辐射过程的探测器实际测量过程的物理特性；b为前向模型所需的其他参数：大气廓线、地表反射率、仪器的光谱响应函数、分子谱线库、太阳光谱等。ε即测量噪声，包括观测过程和前向模拟所造成的系统误差。观测误差会随观测条件的不同而变化，无法获取每次观测的准确误差值，因此无法求解公式(1)，从而转换为在预估的误差范围内寻求最优解。
根据bayesian理论，当观测光谱y已知时，x的概率表示为方程式(4)，其中p(x)表示真实情况下，x发生的概率，表明目标反演量在真实情况下的分布状态，即为先验；p(y)表示在真实条件下，y发生的概率；p(x|y)表示当y已知时，x发生的概率；p(y|x)表示x已知时，y发生的概率。
由最优估计可知，当观测光谱y已知时，概率最大的x即为反演所求的最优解，即可通过p(x|y)求导得到最优解，可用方程式(5)表示，
假设自然界中最一般的描述物理量分布特征的方式是正态分布，则p(x|y)可以表示为方程式(6)，其中sδ和sa表示观测和先验协方差，是先验状态量，c为常数。
由方程式(5)和式(6)可得方程式(7)，其中xa表示最大概率状态量，表示对x求偏导。
为了表征实测光谱与模拟计算光谱的偏差程度，定义代价函数χ(x)，由方程式(8)表示。
然而大气辐射传输模型是非线性的，因此需要通过迭代的方式一步一步寻找局地线性解，逼近最近的反演结果。如果函数f在合理的区间内是连续变化的，那么必然存在一个δx，可以使f在区间(x,x±δx)范围内近似于线性。
levenberg-marquardt方法
当先验值跟真值差距大时，levenberg-marquardt方法使迭代过程、迭代结果更稳定，可由方程式(9)表示。为了保证迭代过程的稳定性，参数λ在每次迭代过程中被调整，确保代价函数平稳的减小，参数λ初始值设为10.0。
反演结果可以分为3类，由表2表示。
表2反演结果描述
在每次更新xi之前，需要评估一下正演模型的非线性度对代价函数减小的影响，非线性度因子r为实际代价函数的变化值与线性假设下代价函数变化值的比值，可由方程式(10)表示。
其中，
(2)模型反演设置
在sciatran模型中，实现反演过程，需要在control.inp和control_ret.inp文件做相应设置。首先在control.inp文件中将反演模式“rtmmode”设置为“ret”状态，开始sciatran的反演过程。进一步在control_ret.inp中，对所需反演的气体及反演内容进行设定，包括偏移压缩量、最高迭代次数、收敛条件等。
表3反演控制文件设置
(3)状态参量初始值
模型输入参数的准确度，即太阳、大气、地表和仪器参数的准确性，直接影响大气辐射传输模型模拟真实大气的准确度，所以使用同步观测的辅助数据更能获取更准确的模拟值。在本实施例中，输入首先使用sciatran模型数据库作为基本输入参数，包括云和气溶胶目前使用其自带cloud数据库和hitran2008数据库，并假定地表为朗伯体，根据波段提取的辐照度计算地表反射率，由tansatl1b提取观测点高程信息，其他大气参数包含温湿压及co2廓线由其他来源的辅助数据提供，以提高模型反演精度，下面介绍来源ecmwf的在分析资料和geos-chem模拟的co2廓线初始值。
1)温湿压廓线数据
ecmwf(europeancentreformedium-rangeweatherforecasts)是欧洲中期天气预报中心，是世界上最大的气象数据中心之一，可提供全球数值天气预报数据等资料。按模型需要，从ecmwf官方网站下载其经过同化和订正得到的再分析数据(温度、湿度、位势、压强)，选取空间分辨率为0.5°×0.5°的大气廓线数据。再根据tansatl1b数据中的经纬度和时间信息进行提取，下面以2018年3月8日北京区域目标观测为例，提取其压强、温度、湿度、位势廓线如图15至图18所示。其中，大气分层为按照大气压强划分37层，包括1、2、3、5、7、10、20、30、50、70、100、125、150、175、200、225、250、300、350、400、450、500、550、600、650、700、750、775、800、825、850、875、900、925、950、975、1000hpa。
2)co2初始廓线数据
针对大气co2浓度廓线，利用geoschem对过去50年进行模拟，获取所反演时段大气co2浓度模拟结果，处理得到2018年4月的co2初始浓度全球分布如图19所示。将其模拟廓线按照压强权重函数(connoretal.,2008)进行插值获取xco2数据，并以10°在纬度带进行划分，以各纬度带xco2均值作为反演co2先验值进行处理，其中先验值的廓线各层相同。
压强权重函数如下：
其中，i为大气层数，通过上下层的大气进行权重求值，用于获取大气压强权重，对于上下边界层分别只取绝对值符号里面的左侧或者右侧，h表示压强权重函数，p代表压强，psurf代表地表压强。
3)云和气溶胶数据产品
在xco2反演过程中，对辐射传输过程中光学厚度影响最大的要素是云和气溶胶问题，本实施例使用相同时期fy-3d大气数据产品的气溶胶光学厚度数据和fy-4a大气数据产品的云检测数据进行筛选，剔除条件分别为云检测结果中表征有云的数据和气溶胶光学厚度大于0.3的数据。
本发明考虑到不同卫星之间、乃至不同传感器之间相同时期观测都会产生不同程度差异，所以一套反演算法无法直接应用于不同卫星观测数据处理中，本发明针对tansat卫星数据属性字段特定信息提取，同时设定tansat最优化参数及辅助数据方案，发展了一套适用于tansat卫星观测的大气二氧化碳浓度反演方法。本发明对tansatl1b数据进行信息提取，并基于风云卫星产品数据及其反演结果进行优化控制，获取xco2数据高精度时空分布。本发明在现有技术的基础上，针对tansat特性进行各类辅助数据的选取及数据处理，得到观测数据反演结果，且精度良好。
本领域技术人员在考虑说明书及实践这里记载的发明后，将容易想到本技术的其它实施方案。本申请旨在涵盖本技术的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本申请的一般性原理并包括本申请未记载的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本项发明的技术性范围并不局限于说明书上的内容，必须要根据本申请的保护范围来确定技术性范围。
应当理解的是，本申请并不局限于上面已经描述并在附图中示出的精确结构，并且可以在不脱离其范围进行各种修改和改变。本申请的范围仅由所附的代表保护范围的内容来限制。
1.一种deei的基于卫星观测的大气二氧化碳浓度反演算法，其特征在于，所述卫星为二氧化碳监测卫星，所述方法包括以下步骤：
步骤1，提取卫星数据中的卫星参数和观测足迹参数；对云和气溶胶数据分别进行有云区域数据剔除和气溶胶高值区域剔除，得到云和气溶胶参数；
步骤2，设置大气辐射传输模型的参数，将所述卫星参数、所述观测足迹参数、所述云和气溶胶参数以及大气廓线参数输入所述大气辐射传输模型进行计算得到模拟光谱；
步骤3，对所述模拟光谱和所述卫星参数中的观测光谱进行最优估计法求解，并判断所述最优估计法求解的结果是否收敛，如果收敛则得到xco2反演结果，否则更新所述大气辐射传输模型中的参数并重复步骤2和步骤3直至所述最优估计法求解的结果收敛。
2.根据权利要求1所述的deei的基于卫星观测的大气二氧化碳浓度反演算法，其特征在于，步骤1中，所述卫星参数包括观测角度、观测时间、仪器信息、光谱信噪比和观测光谱。
3.根据权利要求1所述的deei的基于卫星观测的大气二氧化碳浓度反演算法，其特征在于，步骤1中，所述观测足迹参数包括经纬度、地表反照率、地表高程、观测质量标识以及陆地与水域标识。
4.根据权利要求1所述的deei的基于卫星观测的大气二氧化碳浓度反演算法，其特征在于，步骤1中，进行有云区域数据剔除和气溶胶高值区域剔除过程中，剔除条件分别为云检测结果中表征有云的数据和气溶胶光学厚度大于0.3的数据。
5.根据权利要求1所述的deei的基于卫星观测的大气二氧化碳浓度反演算法，其特征在于，步骤2中，所述大气辐射传输模型中的参数包括观测模式、散射模式以及痕量气体廓线。
6.根据权利要求1所述的deei的基于卫星观测的大气二氧化碳浓度反演算法，其特征在于，步骤2中，所述大气廓线参数包括co2廓线、温度廓线、位势廓线、湿度廓线以及压强廓线。
7.根据权利要求1所述的deei的基于卫星观测的大气二氧化碳浓度反演算法，其特征在于，步骤2中，计算得到模拟光谱时还得到权重函数信息，所述权重函数信息用于所述最优估计法求解。
8.根据权利要求5所述的deei的基于卫星观测的大气二氧化碳浓度反演算法，其特征在于，所述痕量气体廓线为二氧化碳廓线。
9.根据权利要求1所述的deei的基于卫星观测的大气二氧化碳浓度反演算法，其特征在于，所述卫星为tansat卫星。
技术总结