CN111881569B - 二氧化碳柱浓度的反演方法、装置、存储介质及电子设备 - Google Patents

二氧化碳柱浓度的反演方法、装置、存储介质及电子设备 Download PDF

Info

Publication number
CN111881569B
CN111881569B CN202010721395.XA CN202010721395A CN111881569B CN 111881569 B CN111881569 B CN 111881569B CN 202010721395 A CN202010721395 A CN 202010721395A CN 111881569 B CN111881569 B CN 111881569B
Authority
CN
China
Prior art keywords
sample
carbon dioxide
state vector
cost function
determining
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010721395.XA
Other languages
English (en)
Other versions
CN111881569A (zh
Inventor
田向军
金哲
张洪芹
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Atmospheric Physics of CAS
Original Assignee
Institute of Atmospheric Physics of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Atmospheric Physics of CAS filed Critical Institute of Atmospheric Physics of CAS
Priority to CN202010721395.XA priority Critical patent/CN111881569B/zh
Publication of CN111881569A publication Critical patent/CN111881569A/zh
Application granted granted Critical
Publication of CN111881569B publication Critical patent/CN111881569B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Abstract

本发明提供了一种二氧化碳柱浓度的反演方法、装置、存储介质及电子设备，其中，该方法包括：获取目标大气数据和历史大气数据，确定先验状态向量、观测光谱、样本状态向量和模拟光谱；确定反演模型的代价函数，基于非线性最小二乘四维变分同化方法极小化代价函数；在达到收敛标准时，确定此时最优分析场中的二氧化碳廓线分量，并确定二氧化碳柱浓度。通过本发明实施例提供的方案，将传统的数据反演代价函数改写为增量形式的代价函数；基于非线性最小二乘四维变分同化方法迭代求解以极小化代价函数，不需要计算切线性模型与伴随模型，从而极大减少了计算量，且正演模型更新时也不需要重新构建切线性模型和伴随模型，维护、更新成本低。

Description

二氧化碳柱浓度的反演方法、装置、存储介质及电子设备
技术领域
本发明涉及大气监测技术领域，具体而言，涉及一种二氧化碳柱浓度的反演方法、装置、存储介质及电子设备。
背景技术
全球用来监测大气CO2浓度的卫星主要有GOSAT(Greenhouse gases ObservingSATellite)，OCO-2(Orbiting Carbon Observatory-2)以及TanSat。日本于2009年发射的GOSAT卫星是第一颗专门用于温室气体监测的卫星。美国于2014年发射的OCO-2卫星可以获得区域尺度上的碳源汇分布。中国于2016年发射了首颗碳卫星TanSat。目前，从卫星观测的光谱数据得到大气CO2浓度需要通过反演系统实现。反演系统主要包括两个部分，分别为正演模型和反演模型。正演模型由输入数据得到对应的模拟光谱；反演模型通过最优化方法，结合先验信息以及卫星观测数据，得到对当前状态的最优估计，从而推断出大气CO2柱浓度XCO2。
目前针对卫星观测已经开发了多套CO2反演系统，包括应用于GOSAT观测的NIES(National Institute for Environment Studies)系统，同时应用于GOSAT观测以及OCO-2观测的ACOS(Atmospheric CO2 Observations from Space)系统和RemoTeC系统，同时应用于GOSAT观测以及TanSat观测的UoL-FP系统等。需要指出的是，以上系统的反演模型均采用了Levenberg-Marquardt method这一最优化算法，或对其做了微小修改。这些反演模型在运算时需要迭代求解，且在求解过程中无法避免地需要计算雅各比矩阵(正演模型对状态变量的梯度，即切线性模型)与雅各比矩阵的转置(伴随模型)。切线性模型与伴随模型的编码难度与计算复杂度都很高，再加上卫星观测的数据量巨大，在卫星观测的CO2反演过程中这一部分所占的计算量十分庞大。此外，这些系统的切线性模型与伴随模型是针对特定的正演模型构建的，而正演模型又在不断地发展，这导致系统的可移植性较差，正演模型的每一次改进都要重新构建切线性模型与伴随模型，系统的维护、更新成本大大提高。
发明内容
为解决上述问题，本发明实施例的目的在于提供一种二氧化碳柱浓度的反演方法、装置、存储介质及电子设备。
第一方面，本发明实施例提供了一种二氧化碳柱浓度的反演方法，包括：
获取目标大气数据和历史大气数据，根据所述目标大气数据确定先验状态向量和观测光谱，根据所述历史大气数据确定样本状态向量，并根据正演模型确定相应的模拟光谱；
确定反演模型的代价函数，将所述先验状态向量、所述观测光谱、所述样本状态向量和所述模拟光谱作为已知信息，基于非线性最小二乘四维变分同化方法极小化所述代价函数；所述代价函数为基于非线性最小二乘四维变分同化方法对数据反演代价函数改写后所确定的函数；
在所述代价函数迭代后达到收敛标准时，确定此时最优分析场中的二氧化碳廓线分量，并根据所述二氧化碳廓线分量确定二氧化碳柱浓度。
第二方面，本发明实施例还提供了一种二氧化碳柱浓度的反演装置，包括：
获取模块，用于获取目标大气数据和历史大气数据，根据所述目标大气数据确定先验状态向量和观测光谱，根据所述历史大气数据确定样本状态向量，并根据正演模型确定相应的模拟光谱；
处理模块，用于确定反演模型的代价函数，将所述先验状态向量、所述观测光谱、所述样本状态向量和所述模拟光谱作为已知信息，基于非线性最小二乘四维变分同化方法极小化所述代价函数；所述代价函数为基于非线性最小二乘四维变分同化方法对数据反演代价函数改写后所确定的函数；
柱浓度确定模块，用于在所述代价函数迭代后达到收敛标准时，确定此时最优分析场中的二氧化碳廓线分量，并根据所述二氧化碳廓线分量确定二氧化碳柱浓度。
第三方面，本发明实施例还提供了一种计算机存储介质，所述计算机存储介质存储有计算机可执行指令，所述计算机可执行指令用于上述任意一项所述的二氧化碳柱浓度的反演方法。
第四方面，本发明实施例还提供了一种电子设备，包括：
至少一个处理器；以及，
与所述至少一个处理器通信连接的存储器；其中，
所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器能够执行上述任意一项所述的二氧化碳柱浓度的反演方法。
本发明实施例上述第一方面提供的方案中，通过将数据同化中的观测时刻限定为初始时刻、数据同化中的预报模型设置为恒等算子、数据同化中的观测算子取为正演模型等方式，从而将非线性最小二乘四维变分同化方法应用到反演系统中，将传统的数据反演代价函数改写为增量形式的代价函数，并且在极小化代价函数的过程中，基于非线性最小二乘四维变分同化方法进行迭代求解，不需要计算切线性模型与伴随模型，从而极大减少了计算量。即使正演模型发生变化改进，也不需要重新构建新的与正演模型对应的切线性模型与伴随模型，维护、更新成本低。
为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1示出了本发明实施例所提供的一种二氧化碳柱浓度的反演方法的流程图；
图2示出了验证本发明实施例所提供的二氧化碳柱浓度的反演方法中，观测系统模拟实验选取的四个观测点的位置示意图；
图3示出了本发明实施例提供的NEW反演系统，在观测系统模拟实验中经过一次、两个、三次迭代得到的XCO2与先验XCO2以及真实XCO2的对比示意图；
图4示出了本发明实施例提供的NEW反演系统，在观测系统模拟实验中经过一次、两次、三次迭代得到的CO2廓线与先验CO2廓线以及真实CO2廓线的对比示意图；
图5示出了本发明实施例提供的NEW反演系统，在真实实验中XCO2中位数与OCO-2相同观测点以及TCCON对应观测中位数的对比示意图；
图6示出了本发明实施例提供的NEW反演系统，在真实实验中XCO2与TCCON观测的线性拟合；
图7示出了本发明实施例所提供的一种二氧化碳柱浓度的反演装置的结构示意图；
图8示出了本发明实施例所提供的用于执行二氧化碳柱浓度的反演方法的电子设备的结构示意图。
具体实施方式
在本发明的描述中，需要理解的是，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。
本发明实施例提供的一种二氧化碳柱浓度的反演方法，用于确定大气中二氧化碳(CO2)的柱浓度XCO2，参见图1所示，该方法包括：
步骤101：获取目标大气数据和历史大气数据，根据目标大气数据确定先验状态向量和观测光谱，根据历史大气数据确定样本状态向量，并根据正演模型确定相应的模拟光谱。
本发明实施例中，当需要确定某个观测点的二氧化碳柱浓度XCO2时，需要获取到该观测点处的大气数据，即目标大气数据，之后基于该目标大气数据即可确定相应观测光谱(卫星观测光谱)。同时，基于该目标大气数据也可以确定相应的状态向量，即先验状态向量；该状态向量中包含二氧化碳廓线(即CO2廓线，或CO2浓度廓线)以及一些卫星观测光谱对其敏感的参数，如表面气压、地表反射率参数等，根据其中的CO2廓线即可计算得到二氧化碳柱浓度XCO2。其中，该大气数据具体包括卫星提供的Level 1B数据和气象数据，该Level1B数据包含三个光谱波段校准的太阳辐射，指向以及地理位置信息，根据卫星接收到的太阳辐射即可组成观测光谱。该气象数据具体可以插值到该观测点，气象数据包括表面气压、水汽廓线、温度廓线、水平风速等。可选地，可以对大气数据进行预筛选，过滤掉低信噪比以及有云或者多气溶胶的观测，将预筛选后的大气数据作为目标大气数据。
正演模型包含辐射传输模型和辅助模型，具体地，该辐射传输模型为LIDORT模型(Linearized Discrete Ordinate Radiative Transfer)，同时采用了ACOS(AtmosphericCO2 Observations from Space)系统中的辅助模型，辅助模型包括太阳模型、大气模型、地表模型以及仪器模型等。其中，观测光谱y与真实状态向量xt之间的关系可用下式表示：
y＝F(xt,b)+ε；
其中，F为正演模型，b是一系列固定的输入参数，ε包含仪器误差以及正演模型误差。
本实施例中，根据目标大气数据中的太阳辐射可以组成观测光谱y，具体可以指卫星接收到的经过地表、大气反射后的太阳辐射光谱；同时，也可以确定相应的状态向量，即先验状态向量xa；且基于正演模型可以确定相应的模拟光谱，如与先验状态向量xa相对应的模拟光谱F(xa)。同时，本实施例中还确定历史大气数据，该历史大气数据与目标大气数据相似，都为一种大气数据；根据该历史大气数据可以确定另一组状态向量，即样本状态向量；且基于正演模型也可以确定与该样本状态向量对应的模拟光谱，例如，样本状态向量xj的模拟光谱F(xj)。通过该先验状态向量和样本状态向量之间的差值即可确定样本扰动，方便后续基于增量形式的代价函数进行处理。
步骤102：确定反演模型的代价函数，将先验状态向量、观测光谱、样本状态向量和模拟光谱作为已知信息，基于非线性最小二乘四维变分同化方法极小化代价函数；代价函数为基于非线性最小二乘四维变分同化方法对数据反演代价函数改写后所确定的函数。
本发明实施例中所采用的反演模型的代价函数与传统反演模型的Levenberg-Marquardt method不同，本实施例中基于非线性最小二乘四维变分同化方法(NonlinearLeast Squares four-Dimensional Variational data assimilation method，NLS-4DVar)，将数据反演代价函数改写为增量形式的代价函数，之后根据该先验状态向量、观测光谱、样本状态向量和模拟光谱对代价函数进行极小化处理，通过多次迭代使得代价函数收敛。其中，本实施例中也基于非线性最小二乘四维变分同化方法进行迭代求解，以极小化该代价函数；该迭代求解方式不需要切线性和伴随模型，从而极大减少了计算量，能够保证计算简单，且计算精度高。此外，即使正演模型发生变化改进，本实施例中也不需要重新构建与正演模型相对应的新的切线性模型和伴随模型，即本实施例中只需要更新正演模型即可，维护、更新成本低。
其中，NLS-4DVar本身是一种数据同化的方法，与数据反演属于两个不同的领域；这两个不同的领域对应的名词符号含义有所差异，最为明显的是数据反演中的正演模型与数据同化中的观测算子。从物理意义上分析，以辐射传输模型为核心的正演模型是由输入的大气、仪器等参数得到对应的经过地表、大气反射后的太阳辐射，而数据同化中的观测算子一般为将变量(如温度、气压、湿度、风速等)从模式格点插值到观测点的插值算子，两者的作用完全不同，即数据反演领域中的正演模型与数据同化领域中的观测算子不同，很难将二者联系起来，因此一般难以想到将数据同化方法应用于数据反演，尤其是将集合数据同化方法(比如NLS-4DVar)应用于卫星观测二氧化碳柱浓度的反演。申请人在实现本发明创造的过程中，通过数学分析发现，数据同化与数据反演本质上是相通的，两者存在着类似“父与子”的关系。
具体地，数据同化代价函数为：
Figure BDA0002600151510000071
数据反演代价函数为：
Figure BDA0002600151510000072
上述的数据同化代价函数和数据反演代价函数为本领域公知的公式，此处不做详述。本实施例中，将数据反演的代价函数改写为增量形式，以方便适用数据同化方法。其中，将数据反演代价函数改写为增量形式以确定代价函数时，将数据同化中的观测时刻限定为初始时刻，数据同化中的预报模型设置为恒等算子，且数据同化中的观测算子取为正演模型。即观测时刻S＝0，预报模型M＝I(I表示单位矩阵)，观测算子H取为CO2反演中的正演模型F。经过上述适应性修改，从而可以将数据同化方法应用于卫星二氧化碳柱浓度XCO2的反演。具体地，增量形式的代价函数即为本实施例中的代价函数，其为：
Figure BDA0002600151510000081
其中，xa表示状态向量的先验值，即先验状态向量，
Figure BDA0002600151510000082
nx为状态向量的维数，先验状态向量与样本状态向量的维数均为nx；x′为先验状态向量的扰动，且x′＝x-xa；y为根据目标大气数据所确定的观测光谱，
Figure BDA0002600151510000083
ny为观测光谱的维数；Sa为先验协方差矩阵，F(x)表示状态向量x经过正演模型F()模拟得到的模拟光谱，Sε表示观测误差协方差矩阵，T表示矩阵的转置。在迭代过程中，状态向量x是变化的，而观测光谱y固定不变，通过多次迭代以确定状态向量x的最优值。
此外，由于卫星提供的观测误差必定存在偏差，因此在本实施例中通过一个修正系数γ来更准确地刻画观测误差协方差矩阵Sε，得到γSε，这里的修正系数γ具体可通过敏感性试验确定，修正后所确定的代价函数为：
Figure BDA0002600151510000084
步骤103：在代价函数迭代后达到收敛标准时，确定此时最优分析场中的二氧化碳廓线分量，并根据二氧化碳廓线分量确定二氧化碳柱浓度。
本发明实施例中，在极小化代价函数的过程中，每轮过程均可以确定一个分析场，该分析场表示状态向量的分析值，即当前轮迭代所确定的最优的状态向量。当达到收敛标准时，此时所确定的分析场即为最优分析场。其中，该收敛标准为预先设置的，例如，该收敛标准可以为迭代轮数，一般经过三次迭代即可认为达到收敛。该最优分析场中包含二氧化碳廓线分量，根据该二氧化碳廓线分量即可计算得到二氧化碳柱浓度XCO2。
本发明实施例提供的一种二氧化碳柱浓度的反演方法，通过将数据同化中的观测时刻限定为初始时刻、数据同化中的预报模型设置为恒等算子、数据同化中的观测算子取为正演模型等方式，从而将非线性最小二乘四维变分同化方法应用到反演系统中，且对观测误差协方差矩阵添加修正因子，将传统的数据反演代价函数改写为增量形式的代价函数；在极小化代价函数的过程中，基于非线性最小二乘四维变分同化方法进行迭代求解，从而不需要计算切线性模型与伴随模型，极大减少了计算量，能够保证计算简单，且计算精度高。此外，即使正演模型发生变化改进，该方法也不需要重新构建与正演模型相对应的新的切线性模型和伴随模型，维护、更新成本低。
在上述实施例的基础上，由于本实施例在最优化过程中采用的是一种集合的方法，故需要构造一组状态向量的样本，即多个样本状态向量。同时，由于状态向量中二氧化碳廓线分量具有明确的物理意义，且为矢量，而其余参数或者物理意义不够直接，或者为单一标量，因而本实施例中对状态向量中的二氧化碳廓线分量和其余参数采用不同的样本生成方法。具体地，上述步骤101“根据历史大气数据确定样本状态向量”包括：
步骤A1：确定历史大气数据中同一个模式网格点在不同时刻的二氧化碳廓线，并确定日平均二氧化碳廓线。
步骤A2：通过多次样条插值的方式将日平均二氧化碳廓线插值到目标气压层上，目标气压层为先验状态向量中二氧化碳廓线所在的气压层。
本发明实施例中，历史大气数据包含多个模式网格点的大气数据，其中包含每个模式网格点在不同时刻的二氧化碳廓线，基于多个二氧化碳廓线的平均值即可确定该模式网格点的日平均二氧化碳廓线。例如，历史大气数据可以采用Geos-chem模式(v11-01)来模拟得到二氧化碳廓线，并采用由Tan-Tracker(v1)同化OCO-2观测后的二氧化碳浓度和碳通量驱动该模式。其中，Geos-chem模式的水平分辨率为2.5°×2°，垂直气压层为47层，一天内在8个时刻输出结果，分别为03时，06时，09时，12时，15时，18时，21时和24时，即每个模式网格点具有8个二氧化碳廓线，对一天内8个时刻的CO2廓线做平均，得到日平均CO2廓线。
同时，在构建历史大气数据作为样本时，样本个数不宜过多，本实施例中从一个或多个纬度带中，间隔选取多个模式网格点的大气数据形成历史大气数据；同时，由于CO2具有纬度变化特征，所选取的纬度带为与目标观测点相邻的纬度带，该目标观测点为目标大气数据所对应的观测点，即选取目标观测点上下两个纬度带的模式模拟结果来体现纬度变化特征。例如，经度方向每4个模式格点抽取一个结果，即Geos-chem的水平分辨率取为10°×2°，然后将日平均二氧化碳廓线通过多次(如三次)样条插值方法(自然边界条件)插值到先验状态向量中的二氧化碳廓线所在气压层上，即插值到目标气压层上，该目标气压层为20层气压层。由于二氧化碳廓线存在季节变化与纬度变化特征，根据观测点测量的月份和所在的纬度，选取测量月份一整个月的，观测点所在纬度上下两个纬度带的Geos-chem模式模拟结果作为原始的大样本，大样本个数一般在2000-2300个。以2016年1月30日位于36.641N°，97.441W°的观测点为例，将Geos-chem的模拟数据中2016年1月36.0N°，38.0N°两个纬度带上共36×2×31＝2232个日平均二氧化碳廓线作为大样本。本实施例中，基于上述步骤A1-A2即可确定每个模式网格点在目标气压层上的日平均二氧化碳廓线。
步骤A3：将一个或多个纬度带中的多个模式网格点所对应的插值后的日平均二氧化碳廓线作为大样本，通过奇异值分解和随机正交矩阵从大样本中选取出多个小样本。
本发明实施例中，由于大样本数量较多，影响计算效率，本实施例从大样本中选取出部分作为小样本，从而可以减少样本数量；其中，每个小样本即为一个二氧化碳廓线。其中，小样本的样本数量可通过敏感性试验测试确定。本实施例中，通过敏感性试验测试了样本个数为30，50，80，100个的反演结果，最终选定样本个数取50个，即涉及50个样本二氧化碳廓线。具体地，可以通过奇异值分解和随机正交矩阵从大样本中选取出多个小样本，其过程为现有的成熟方案，此处不做赘述。其中，“大样本”和“小样本”为相对的概念，其中的“大”和“小”仅用于表示二者之间的大小关系，即大样本所包含的样本数量大于小样本所包含的样本数量。上述样本库的生成方式，计算经济有效且反演精度高
步骤A4：在其余参数的先验值上添加正态分布的随机扰动确定其余参数的样本值，先验值为先验状态向量中与其余参数相关的值。
步骤A5：根据样本二氧化碳廓线和其余参数的样本值确定样本状态向量。
本发明实施例中，如上所述，状态向量中包含二氧化碳廓线分量，还包含其余参数的分量，该其余参数指的是除二氧化碳廓线之外的其他参数，具体包括：温度廓线订正因子、表面气压、水汽廓线缩放因子、气溶胶参数、地表反照率参数、仪器订正参数、光谱订正参数、太阳诱导的荧光参数等。为了方便确定其余参数的数值，本实施例中利用先验状态向量中相应的先验值，通过添加正态分布的随机扰动来确定。其中，根据先验状态向量中该其余参数的先验值，并将该先验值作为正态分布的均值，并预先确定每个其余参数所对应的正态分布的标准差，进而可以添加正态分布的随机扰动，并确定其余参数的数值，即样本值。最后将样本二氧化碳廓线和其余参数的样本值合并到一起，即可生成状态向量的样本，即样本状态向量。该样本状态向量用于在后续处理过程中确定初始样本扰动。
其中，状态向量中部分其余参数的标准差可参见下表1所示：
表1
Figure BDA0002600151510000111
Figure BDA0002600151510000121
此外，考虑到部分参数的物理性质，某些参数在正态分布生成样本的过程中还需要对其施加额外的约束，即上述步骤A4“在其余参数的先验值上添加正态分布的随机扰动确定其余参数的样本值”包括：对正态分布的随机扰动施加约束条件，生成其余参数的样本值。其中，约束条件包括：水汽廓线缩放因子约束条件、第一气溶胶约束条件、第二气溶胶约束条件、地表反射参数约束条件、风速约束条件中的一种或多种。即，当添加正态分布的随机扰动所确定的样本值符合该约束条件时，说明该样本值异常，需要重新确定，此时可以基于均匀分布添加随机扰动。
本实施例中，水汽廓线缩放因子约束条件为：若水汽廓线缩放因子的样本值小于0，则舍弃相应的样本值，重新赋值为0到水汽廓线缩放因子的先验值之间均匀分布的随机数。
第一气溶胶约束条件为：两种对流层气溶胶在预设波长(如755nm)处的气溶胶光学厚度的自然对数若大于第一临界厚度(如-1.6095)，则舍弃相应的样本值，重新赋值为第一气溶胶光学厚度的自然对数的先验值到第一临界厚度之间均匀分布的随机数。
第二气溶胶约束条件为：水云、冰云或平流层气溶胶在预设波长(如755nm)处的气溶胶光学厚度的自然对数若大于等于第二临界厚度(如0)，则舍弃相应的样本值，重新赋值为第二气溶胶光学厚度的自然对数的先验值到第二临界厚度之间均匀分布的随机数。
地表反射参数约束条件为：朗伯地表反射参数若小于等于预设阈值(如0)，则舍弃相应的样本值，重新赋值为预设阈值到朗伯地表反射参数的先验值之间均匀分布的随机数。
风速约束条件为：海洋上的风速若小于等于0，则舍弃相应的样本值，重新赋值为0到海洋风速的先验值之间均匀分布的随机数。
本发明实施例中，样本状态向量中的二氧化碳廓线和其余参数分别确定，通过上述小样本的生成方法可以确定样本二氧化碳廓线，计算经济有效且反演精度高；在先验值的基础上增加正态分布的随机扰动，可以方便快速地确定其余参数的样本值，进而生成样本状态向量。
在上述实施例基础上，上述步骤102“基于非线性最小二乘四维变分同化方法极小化代价函数”包括：
步骤B1：根据样本状态向量和先验状态向量确定初始样本扰动Px，且Px＝(x′1,x′2,...,x′N)；xj′表示第j个初始样本扰动，x′j＝xj-xa(j＝1,2,...,N)，xj表示第j个样本状态向量，N为样本状态向量个数。
步骤B2：将分析增量x′，*转换为初始样本扰动Px的线性组合，即x′，*＝Pxβ，其中β＝(β1,β2,...,βN)T为权重系数向量。
步骤B3：用样本协方差矩阵Be代替先验协方差矩阵Sa，且
Figure BDA0002600151510000131
步骤B4：将代价函数转换为以权重系数向量β为自变量的非线性最小二乘形式：
Figure BDA0002600151510000132
其中，
Figure BDA0002600151510000133
且
Figure BDA0002600151510000134
Sε,+表示经过平方根法分解(Cholesky分解)得到的矩阵，F′(Pxβ)＝F(xa+Pxβ)-F(xa)，y′＝y-F(xa)。
步骤B5：基于高斯牛顿迭代法极小化非线性最小二乘形式的代价函数，权重系数更新公式为：
Figure BDA0002600151510000141
Figure BDA0002600151510000142
xi,*＝xi-1,*+δxi-1,*；
其中，xi,*为当前迭代的分析场，xi-1,*为上一迭代的分析场，x′,i-1,*＝xi-1,*-xa为上一迭代的分析增量；δβi-1为当前迭代样本扰动
Figure BDA00026001515100001410
的权重系数向量，δxi-1,*表示相较于上一次迭代的分析场的增量，y′,i-1＝y-F(xi-1,*)；
其中，
Figure BDA0002600151510000143
Figure BDA0002600151510000144
I表示单位矩阵，
Figure BDA0002600151510000145
Figure BDA0002600151510000146
Figure BDA0002600151510000147
本发明实施例中，如上所述，代价函数为：
Figure BDA0002600151510000148
同时，定义F′(x′)＝F(xa+x′)-F(xa)，y′＝y-F(xa)；其中，y为根据目标大气数据所确定的观测光谱，
Figure BDA0002600151510000149
ny为观测光谱的维数。为了避免对先验协方差矩阵Sa的直接求逆，本实施例提供的NLS-4DVar方法中，将分析增量x′，*表示为初始样本扰动Px的线性组合，即x′，*＝Pxβ，其中，Px＝(x′1,x′2,...,x′N)为初始样本扰动，x′j＝xj-xa表示第j个样本的初始扰动，j＝1,2,...,N，β＝(β1,β2,...,βN)T为权重系数向量。其中，该分析增量x′，*为上述先验状态向量的扰动x′的一种特殊情况，即分析增量也为一种扰动，且为一种最优的扰动，可以用分析增量x′，*代替扰动x′以进行后续处理。
进一步地，假设先验协方差矩阵Sa可以用样本协方差矩阵Be(ensemblecovariance matrix)近似代替，即用样本协方差矩阵Be代替先验协方差矩阵Sa，且
Figure BDA0002600151510000151
在接下来的计算过程中，由于样本扰动矩阵发生变化，使得Sa随着迭代的进行不断变化，这也使得它对真实误差的刻画更加准确，而传统方法中Sa是静态的不变的。
将上式样本协方差矩阵Be带入到代价函数中即可得到新的代价函数：
Figure BDA0002600151510000152
为了避免计算切线性模型与伴随模型，NLS-4DVar方法将上式新的代价函数转换为如下非线性最小二乘形式：
Figure BDA0002600151510000153
其中，
Figure BDA0002600151510000154
F′(Pxβ)表示状态向量分析值xa+Pxβ的模拟光谱F(xa+Pxβ)与先验值xa的模拟光谱F(xa)之间的差值。使用高斯-牛顿迭代法极小化该代价函数，并经过简单的数学推导，很好地避免了切线性模型与伴随模型的计算，并最终得到能够更新权重系数向量的系数更新公式为：
Figure BDA0002600151510000155
Figure BDA0002600151510000156
xi,*＝xi-1,*+δxi-1,*； (3)
其中，xi,*为当前迭代的分析场，xi-1,*为上一迭代的分析场；x′,i-1,*＝xi-1,*-xa为上一迭代的分析增量；δβi-1整体为一个参数，其表示当前迭代样本扰动
Figure BDA00026001515100001610
的权重系数向量，δxi-1,*整体也为一个参数，其表示相较于上一次迭代的分析场xi-1,*的增量。y′,i-1＝y-F(xi -1,*)，其为观测光谱y与上一次迭代分析场模拟光谱F(xi-1,*)的差值；且：
Figure BDA0002600151510000161
Figure BDA0002600151510000162
其中，
Figure BDA0002600151510000163
Figure BDA0002600151510000164
Figure BDA0002600151510000165
Figure BDA0002600151510000166
本发明实施例中，该代价函数一般在三次迭代后达到极小化收敛标准，此时即可得到的分析场xi,*中包含二氧化碳廓线分量，根据该二氧化碳廓线可计算得到二氧化碳柱浓度XCO2。
具体地，在第一轮迭代中，i＝1，上一迭代的分析场xi-1,*为所述先验状态向量xa，即xi-1,*＝x0,*＝xa；当前迭代样本扰动
Figure BDA0002600151510000167
为所述初始样本扰动Px，即
Figure BDA0002600151510000168
上一迭代的分析增量x′,i-1,*为0，即x′,i-1,*＝x′,0,*＝x0,*-xa＝0。且此时F(xi-1,*)＝F(xa)，即先验状态向量的模拟光谱；F(xa+x′j)＝F(xj)，即样本状态向量的模拟光谱。因此，在第一轮迭代中，先验状态向量、观测光谱、样本状态向量和模拟光谱均为已知信息，根据上述的公式(6)、(7)等可以确定此时的
Figure BDA0002600151510000169
进而基于公式(4)确定
Figure BDA0002600151510000171
由于第一轮迭代x′,i-1,*＝0，故此时基于公式(1)即可确定第一轮迭代后更新的系数δβ0，最后基于公式(2)、(3)最终确定第一轮的分析场x1,*。之后即可进行下一轮迭代处理，直至达到收敛标准。
可选地，在步骤103“根据二氧化碳廓线分量确定二氧化碳柱浓度”之后，该方法还包括：基于订正方式对二氧化碳柱浓度进行偏差订正，确定最终二氧化碳柱浓度；订正方式包括：依赖于观测足迹的订正、和其他反演变量之间的虚假相关引起的偏差订正、和地基观测对比的缩放因子订正中的一种或多种。本发明实施例中，由于天基观测的XCO2均存在系统误差，主要来源于观测信息的不全面，仪器误差以及模型误差，因而需要对XCO2进行偏差订正。
本发明实施例提供的一种二氧化碳柱浓度的反演方法，通过将数据同化中的观测时刻限定为初始时刻、数据同化中的预报模型设置为恒等算子、数据同化中的观测算子取为正演模型等方式，从而将非线性最小二乘四维变分同化方法应用到反演系统中，将传统的数据反演代价函数改写为增量形式的代价函数；在极小化代价函数的过程中，基于非线性最小二乘四维变分同化方法进行迭代求解，从而不需要计算切线性模型与伴随模型，从而极大减少了计算量，能够保证计算简单，且计算精度高。即使正演模型发生变化改进，也不需要重新构建与正演模型对应的新的切线性模型和伴随模型，维护、更新成本低。样本状态向量中的二氧化碳廓线和其余参数分别确定，通过上述小样本的生成方法可以确定样本二氧化碳廓线，计算经济有效且反演精度高；在先验值的基础上增加正态分布的随机扰动，可以方便快速地确定其余参数的样本值，进而生成样本状态向量。用样本协方差矩阵Be代替先验协方差矩阵Sa，使得Sa随着迭代的进行不断变化，这也使得它对真实误差的刻画更加准确。
下面通过观测系统模拟实验以及真实实验对上述方法进行验证，以证明本方法的有效性，并证明其一定程度上具有优越性。本实施例中，基于该二氧化碳柱浓度的反演方法生成相应的反演系统，本实施例中简称为NEW反演系统。
在观测系统模拟实验中，在全球共选择了四个观测点做反演，其中三个陆地观测点选在了TCCON(Total Carbon Column Observing Network)观测站上空，分别为Lamont站，Bremen站以及Wollongong站，还有一个点选在了北太平洋洋面上。将这四个观测点分别命名为la点，br点，wo点以及oc点，这个四个点的信息如表2所示，具体位置参见图2。四个点的观测日期分别为2016年1月3日，2016年3月17日，2016年2月10日以及2016年7月1日。在选取观测点时同时考虑了不同下垫面属性，不同纬度，海陆分布对反演结果的影响，以保证实验的全面性。
表2：选取的四个观测点的信息
Figure BDA0002600151510000181
实验中用到的卫星数据(即目标大气数据)为version 8r版本的OCO-2数据，包括Level 1B数据OCO2_L1B_Calibration 8r,气象数据OCO2_L2_Met 8r,预筛选数据OCO2_L2_ABand 8r和OCO2_IMAPDOAS 8r,以及分析数据OCO2_L2_Diagnostic 8r；其中，预筛选数据是剔除卫星观测中不好的观测点。三个陆地上的观测点选择目标观测模式的数据，海洋上的点选取耀斑观测模式的数据。经敏感性实验测试后样本个数取为50。
观测系统模拟实验中，状态向量的先验值采用OCO-2根据输入数据计算得到的先验值，状态向量的真值采用该观测点或附近点的OCO-2反演结果。将构造的真值输入正演模型，得到真值模拟的光谱，在该光谱上加上OCO-2的观测误差，作为观测。从XCO2和CO2廓线两方面对结果进行分析，XCO2是目前卫星提供的产品数据，准确的CO2廓线则是反演系统努力的目标。
本实施例中，图3与图4分别展示了先验值、真值以及本实施例提供的方法经过一次、两次、三次迭代得到的XCO2和CO2廓线。图3和图4分别展示了(a)la点,(b)br点,(c)wo点以及(d)oc点的结果，
对于XCO2的结果。设计实验时，先验的XCO2与真实XCO2的差异构造的是比较大的，四个地点的差值分别为1.88ppm,1.70ppm,-3.15ppm,2.27ppm。从图3可以看出，在本实施例提供的方法的作用下，随着迭代的进行反演结果逐渐向真值靠近，并在第三次迭代后达到最优，此时反演的XCO2与真值的差异分别是0.33ppm,0.11ppm,0.11ppm,-1.13ppm，较先验值有了大幅改进，这说明本实施例提供的方法在反演XCO2方面是有效的。
对于CO2廓线的反演情况。从图4中可以看出，CO2廓线存在南北半球差异，北半球的CO2廓线垂直梯度较大，南半球的垂直梯度较小。先看la点和br点这两个北半球点的结果。OCO-2计算得到的先验CO2廓线比较平稳，垂直梯度较小。在本实施例提供的方法的作用下，随着迭代的进行反演的CO2廓线逐步靠近真值，特别是在对流层下层，近地面附近，反演结果与真值十分接近。这一点尤为重要，因为用户关心的碳源汇都集中在近地面附近。再来看wo点和oc点这两个南半球点的反演结果。可以看出，这两个点反演的CO2廓线能够大致把握住真值的特点，和先验值相比有了较大改进，但结果不如北半球的理想。这可能是由于先验的CO2廓线垂直梯度很小，在先验值基础上分解得到的CO2廓线样本的垂直梯度也较小，而分析值是样本的线性组合，因而无法反演出大的垂直梯度的情况。
在真实实验中，考虑到南北半球差异以及海陆性质差异，选取的两个目标观测地点为北美洲平原的Lamont站以及大洋洲沿岸的Darwin站。这两个地点的信息如表3所示。在两个地点各选取了约一年的观测数据，数据的观测日期见表4。
表3：验证的目标观测点，包括该点的经纬度以及海拔(km)
Figure BDA0002600151510000201
表4：真实反演选取站点与日期
Figure BDA0002600151510000202
受计算量的制约，Lamont位置每个日期在OCO-2通过预筛选的观测点中随机选取10个点做反演，Darwin位置从陆地和海洋观测中各随机选取10个点做反演。此时认为在OCO-2目标观测这段时间±30分钟内的TCCON观测和OCO-2的观测是同时发生的，如果这段时间内TCCON观测少于5个，这段时间将被拓宽至±120分钟。然后将选取的10个点反演结果的中位数与TCCON观测的中位数进行对比。真实实验中样本个数取为50个，迭代次数为3次。反演结束后用OCO-2的偏差订正方法对结果进行了校正。
图5展示了NEW反演系统在Lamont站以及Darwin站反演的XCO2的中位数，相同的观测点OCO-2通过ACOS系统反演的XCO2中位数以及对应的TCCON观测的中位数，图5中包括(a)Lamont站结果(b)Darwin站陆地点结果(c)Darwin站海洋点结果。为了更好地评估反演结果，定义偏差(b)为反演的XCO2的中位数与TCCON观测中位数的差值的平均值，观测精度(σ)为差值的标准差。Lamont站点通过NEW反演系统(图5中以NEW表示)得到的b为-1.93ppm，σ为0.91ppm；OCO-2反演得到的b为-2.95ppm，σ为1.49ppm。Darwin站点NEW反演系统计算得到的b为-2.60ppm，σ为1.30ppm；OCO-2得到的b为-3.10ppm，σ为1.31ppm。可以看出，NEW反演系统的反演结果普遍要比OCO-2的结果更加接近TCCON观测，这也说明了NEW反演系统的有效性，以及和OCO-2的结果是可比的。
图6展示了NEW反演系统得到的XCO2与TCCON观测线性拟合的结果，其中(a)图为偏差订正前的结果，(b)图为偏差订正后的结果。偏差订正前NEW反演系统的XCO2与TCCON观测的相关系数为0.89，偏差订正后为0.90，均通过了显著性水平为5％的显著性检验。这从一定程度上反映了反演结果是可信的，偏差订正是有效的。
以上详细介绍了大气柱浓度的反演方法流程，该方法也可以通过相应的装置实现，下面详细介绍该装置的结构和功能。
基于同样的发明构思，本发明实施例还提供一种大气柱浓度的反演装置，参见图7所示，该装置包括：
获取模块71，用于获取目标大气数据和历史大气数据，根据所述目标大气数据确定先验状态向量和观测光谱，根据所述历史大气数据确定样本状态向量，并根据正演模型确定相应的模拟光谱；
处理模块72，用于确定反演模型的代价函数，将所述先验状态向量、所述观测光谱、所述样本状态向量和所述模拟光谱作为已知信息，基于非线性最小二乘四维变分同化方法极小化所述代价函数；所述代价函数为基于非线性最小二乘四维变分同化方法对数据反演代价函数改写后所确定的函数；
柱浓度确定模块73，用于在所述代价函数迭代后达到收敛标准时，确定此时最优分析场中的二氧化碳廓线分量，并根据所述二氧化碳廓线分量确定二氧化碳柱浓度。
在上述实施例的基础上，所述获取模块71根据所述历史大气数据确定样本状态向量包括：
确定历史大气数据中同一个模式网格点在不同时刻的二氧化碳廓线，并确定日平均二氧化碳廓线；
通过多次样条插值的方式将所述日平均二氧化碳廓线插值到目标气压层上，所述目标气压层为所述先验状态向量中二氧化碳廓线所在的气压层；
将一个或多个纬度带中的多个模式网格点所对应的插值后的日平均二氧化碳廓线作为大样本，通过奇异值分解和随机正交矩阵从所述大样本中选取出多个小样本；
在其余参数的先验值上添加正态分布的随机扰动确定所述其余参数的样本值，所述先验值为所述先验状态向量中与所述其余参数相关的值；
根据所述样本二氧化碳廓线和其余参数的所述样本值确定样本状态向量。
在上述实施例的基础上，所述获取模块71在其余参数的先验值上添加正态分布的随机扰动确定所述其余参数的样本值，包括：
预先给定其余参数的标准差，将所述其余参数的先验值作为均值，对正态分布的随机扰动施加约束条件，生成所述其余参数的样本值；所述约束条件包括：水汽廓线缩放因子约束条件、第一气溶胶约束条件、第二气溶胶约束条件、地表反射参数约束条件、风速约束条件中的一种或多种；
其中，所述水汽廓线缩放因子约束条件为：若水汽廓线缩放因子的样本值小于0，则舍弃相应的样本值，重新赋值为0到水汽廓线缩放因子的先验值之间均匀分布的随机数；
所述第一气溶胶约束条件为：两种对流层气溶胶在预设波长处的气溶胶光学厚度的自然对数若大于第一临界厚度，则舍弃相应的样本值，重新赋值为第一气溶胶光学厚度的自然对数的先验值到所述第一临界厚度之间均匀分布的随机数；
所述第二气溶胶约束条件为：水云、冰云或平流层气溶胶在预设波长处的气溶胶光学厚度的自然对数若大于等于第二临界厚度，则舍弃相应的样本值，重新赋值为第二气溶胶光学厚度的自然对数的先验值到所述第二临界厚度之间均匀分布的随机数；
所述地表反射参数约束条件为：朗伯地表反射参数若小于等于预设阈值，则舍弃相应的样本值，重新赋值为所述预设阈值到朗伯地表反射参数的先验值之间均匀分布的随机数；
所述风速约束条件为：海洋上的风速若小于等于0，则舍弃相应的样本值，重新赋值为0到海洋风速的先验值之间均匀分布的随机数。
在上述实施例的基础上，所述代价函数为：
Figure BDA0002600151510000231
其中，xa表示先验状态向量，
Figure BDA0002600151510000232
nx为状态向量的维数；x′为先验状态向量的扰动，且x′＝x-xa；y为根据目标大气数据所确定的观测光谱，
Figure BDA0002600151510000233
ny为观测光谱的维数；Sa为先验协方差矩阵，F(x)表示状态向量x经过正演模型F()模拟得到的模拟光谱，Sε表示观测误差协方差矩阵，γ为修正系数，T表示矩阵的转置。
在上述实施例的基础上，所述处理模块72根据所述先验状态向量、所述观测光谱、所述模拟光谱和所述样本状态向量极小化所述代价函数，包括：
根据所述样本状态向量和所述先验状态向量确定初始样本扰动Px，且Px＝(x′1,x′2,...,x′N)；xj′表示第j个初始样本扰动，x′j＝xj-xa(j＝1,2,...,N)，xj表示第j个样本状态向量，N为样本状态向量个数；
将分析增量x′，*转换为初始样本扰动Px的线性组合，即x′，*＝Pxβ，其中β＝(β1,β2,...,βN)T为权重系数向量；
用样本协方差矩阵Be代替先验协方差矩阵Sa，且
Figure BDA0002600151510000241
将代价函数转换为以权重系数向量β为自变量的非线性最小二乘形式：
Figure BDA0002600151510000242
其中，
Figure BDA0002600151510000243
且
Figure BDA0002600151510000244
Sε,+表示经过平方根法分解得到的矩阵，F′(Pxβ)＝F(xa+Pxβ)-F(xa)，y′＝y-F(xa)；
基于高斯牛顿迭代法极小化非线性最小二乘形式的代价函数，且权重系数更新公式为：
Figure BDA0002600151510000245
Figure BDA0002600151510000246
xi,*＝xi-1,*+δxi-1,*；
其中，xi,*为当前迭代的分析场，xi-1,*为上一迭代的分析场，x′,i-1,*＝xi-1,*-xa为上一迭代的分析增量；δβi-1为当前迭代样本扰动
Figure BDA0002600151510000247
的权重系数向量，δxi-1,*表示相较于上一次迭代的分析场的增量，y′,i-1＝y-F(xi-1,*)；
且其中，
Figure BDA0002600151510000251
Figure BDA0002600151510000252
I为单位矩阵，
Figure BDA0002600151510000253
Figure BDA0002600151510000254
Figure BDA0002600151510000255
在上述实施例的基础上，在第一轮迭代中，上一迭代的分析场xi-1,*为所述先验状态向量xa，当前迭代样本扰动
Figure BDA0002600151510000256
为所述初始样本扰动Px，上一迭代的分析增量x′,i-1,*为0。
在上述实施例的基础上，该装置还包括订正模块；
在所述柱浓度确定模块73根据所述二氧化碳廓线分量确定二氧化碳柱浓度之后，所述订正模块用于：基于订正方式对所述二氧化碳柱浓度进行偏差订正，确定最终二氧化碳柱浓度；所述订正方式包括：依赖于观测足迹的订正、和其他反演变量之间的虚假相关引起的偏差订正、和地基观测对比的缩放因子订正中的一种或多种。
本发明实施例还提供了一种计算机存储介质，所述计算机存储介质存储有计算机可执行指令，其包含用于执行上述的大气柱浓度的反演方法的程序，该计算机可执行指令可执行上述任意方法实施例中的方法。
其中，所述计算机存储介质可以是计算机能够存取的任何可用介质或数据存储设备，包括但不限于磁性存储器(例如软盘、硬盘、磁带、磁光盘(MO)等)、光学存储器(例如CD、DVD、BD、HVD等)、以及半导体存储器(例如ROM、EPROM、EEPROM、非易失性存储器(NANDFLASH)、固态硬盘(SSD))等。
图8示出了本发明的另一个实施例的一种电子设备的结构框图。所述电子设备1100可以是具备计算能力的主机服务器、个人计算机PC、或者可携带的便携式计算机或终端等。本发明具体实施例并不对电子设备的具体实现做限定。
该电子设备1100包括至少一个处理器(processor)1110、通信接口(Communications Interface)1120、存储器(memory array)1130和总线1140。其中，处理器1110、通信接口1120、以及存储器1130通过总线1140完成相互间的通信。
通信接口1120用于与网元通信，其中网元包括例如虚拟机管理中心、共享存储等。
处理器1110用于执行程序。处理器1110可能是一个中央处理器CPU，或者是专用集成电路ASIC(Application Specific Integrated Circuit)，或者是被配置成实施本发明实施例的一个或多个集成电路。
存储器1130用于可执行的指令。存储器1130可能包含高速RAM存储器，也可能还包括非易失性存储器(non-volatile memory)，例如至少一个磁盘存储器。存储器1130也可以是存储器阵列。存储器1130还可能被分块，并且所述块可按一定的规则组合成虚拟卷。存储器1130存储的指令可被处理器1110执行，以使处理器1110能够执行上述任意方法实施例中的大气柱浓度的反演方法。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (8)

1.一种二氧化碳柱浓度的反演方法，其特征在于，包括：
获取目标大气数据和历史大气数据，根据所述目标大气数据确定先验状态向量和观测光谱，根据所述历史大气数据确定样本状态向量，并根据正演模型确定相应的模拟光谱；
确定反演模型的代价函数，将所述先验状态向量、所述观测光谱、所述样本状态向量和所述模拟光谱作为已知信息，基于非线性最小二乘四维变分同化方法极小化所述代价函数；所述代价函数为基于非线性最小二乘四维变分同化方法对数据反演代价函数改写后所确定的函数；所述代价函数为：
Figure FDA0002860665260000011
其中，xa表示先验状态向量，
Figure FDA0002860665260000012
nx为状态向量的维数；x′为先验状态向量的扰动，且x′＝x-xa；y为根据目标大气数据所确定的观测光谱，
Figure FDA0002860665260000013
ny为观测光谱的维数；Sa为先验协方差矩阵，F(x)表示状态向量x经过正演模型F()模拟得到的模拟光谱，Sε表示观测误差协方差矩阵，γ为修正系数，T表示矩阵的转置；
在所述代价函数迭代后达到收敛标准时，确定此时最优分析场中的二氧化碳廓线分量，并根据所述二氧化碳廓线分量确定二氧化碳柱浓度；
其中，所述基于非线性最小二乘四维变分同化方法极小化所述代价函数，包括：
根据所述样本状态向量和所述先验状态向量确定初始样本扰动Px，且Px＝(x′1,x′2,...,x′N)；xj′表示第j个初始样本扰动，x′j＝xj-xa(j＝1,2,...,N)，xj表示第j个样本状态向量，N为样本状态向量个数；
将分析增量x′，*转换为初始样本扰动Px的线性组合，即x′，*＝Pxβ，其中β＝(β1,β2,...,βN)T为权重系数向量；
用样本协方差矩阵Be代替先验协方差矩阵Sa，且
Figure FDA0002860665260000021
将代价函数转换为以权重系数向量β为自变量的非线性最小二乘形式：
Figure FDA0002860665260000022
其中，
Figure FDA0002860665260000023
且
Figure FDA0002860665260000024
Sε,+表示经过平方根法分解得到的矩阵，F′(Pxβ)＝F(xa+Pxβ)-F(xa)，y′＝y-F(xa)；
基于高斯牛顿迭代法极小化非线性最小二乘形式的代价函数，且权重系数更新公式为：
Figure FDA0002860665260000025
Figure FDA0002860665260000026
xi,*＝xi-1,*+δxi-1,*；
其中，xi,*为当前迭代的分析场，xi-1,*为上一迭代的分析场，x′,i-1,*＝xi-1,*-xa为上一迭代的分析增量；δβi-1为当前迭代样本扰动
Figure FDA0002860665260000027
的权重系数向量，δxi-1,*表示相较于上一次迭代的分析场的增量，y′,i-1＝y-F(xi-1,*)；
且其中，
Figure FDA0002860665260000028
Figure FDA0002860665260000029
I为单位矩阵，
Figure FDA00028606652600000210
Figure FDA00028606652600000211
Figure FDA00028606652600000212
Figure FDA0002860665260000031
2.根据权利要求1所述的方法，其特征在于，所述根据所述历史大气数据确定样本状态向量包括：
确定历史大气数据中同一个模式网格点在不同时刻的二氧化碳廓线，并确定日平均二氧化碳廓线；
通过多次样条插值的方式将所述日平均二氧化碳廓线插值到目标气压层上，所述目标气压层为所述先验状态向量中二氧化碳廓线所在的气压层；
将一个或多个纬度带中的多个模式网格点所对应的插值后的日平均二氧化碳廓线作为大样本，通过奇异值分解和随机正交矩阵从所述大样本中选取出多个小样本；
在其余参数的先验值上添加正态分布的随机扰动确定所述其余参数的样本值，所述先验值为所述先验状态向量中与所述其余参数相关的值；
根据所述样本二氧化碳廓线和其余参数的所述样本值确定样本状态向量。
3.根据权利要求2所述的方法，其特征在于，所述在其余参数的先验值上添加正态分布的随机扰动确定所述其余参数的样本值，包括：
预先给定其余参数的标准差，将所述其余参数的先验值作为均值，对正态分布的随机扰动施加约束条件，生成所述其余参数的样本值；所述约束条件包括：水汽廓线缩放因子约束条件、第一气溶胶约束条件、第二气溶胶约束条件、地表反射参数约束条件、风速约束条件中的一种或多种；
其中，所述水汽廓线缩放因子约束条件为：若水汽廓线缩放因子的样本值小于0，则舍弃相应的样本值，重新赋值为0到水汽廓线缩放因子的先验值之间均匀分布的随机数；
所述第一气溶胶约束条件为：两种对流层气溶胶在预设波长处的气溶胶光学厚度的自然对数若大于第一临界厚度，则舍弃相应的样本值，重新赋值为第一气溶胶光学厚度的自然对数的先验值到所述第一临界厚度之间均匀分布的随机数；
所述第二气溶胶约束条件为：水云、冰云或平流层气溶胶在预设波长处的气溶胶光学厚度的自然对数若大于等于第二临界厚度，则舍弃相应的样本值，重新赋值为第二气溶胶光学厚度的自然对数的先验值到所述第二临界厚度之间均匀分布的随机数；
所述地表反射参数约束条件为：朗伯地表反射参数若小于等于预设阈值，则舍弃相应的样本值，重新赋值为所述预设阈值到朗伯地表反射参数的先验值之间均匀分布的随机数；
所述风速约束条件为：海洋上的风速若小于等于0，则舍弃相应的样本值，重新赋值为0到海洋风速的先验值之间均匀分布的随机数。
4.根据权利要求1所述的方法，其特征在于，在第一轮迭代中，上一迭代的分析场xi-1,*为所述先验状态向量xa，当前迭代样本扰动
Figure FDA0002860665260000041
为所述初始样本扰动Px，上一迭代的分析增量x′,i-1,*为0。
5.根据权利要求1-4任意一项所述的方法，其特征在于，在所述根据所述二氧化碳廓线分量确定二氧化碳柱浓度之后，还包括：
基于订正方式对所述二氧化碳柱浓度进行偏差订正，确定最终二氧化碳柱浓度；所述订正方式包括：依赖于观测足迹的订正、和地基观测对比的缩放因子订正中的一种或多种。
6.一种二氧化碳柱浓度的反演装置，其特征在于，包括：
获取模块，用于获取目标大气数据和历史大气数据，根据所述目标大气数据确定先验状态向量和观测光谱，根据所述历史大气数据确定样本状态向量，并根据正演模型确定相应的模拟光谱；
处理模块，用于确定反演模型的代价函数，将所述先验状态向量、所述观测光谱、所述样本状态向量和所述模拟光谱作为已知信息，基于非线性最小二乘四维变分同化方法极小化所述代价函数；所述代价函数为基于非线性最小二乘四维变分同化方法对数据反演代价函数改写后所确定的函数；所述代价函数为：
Figure FDA0002860665260000051
其中，xa表示先验状态向量，
Figure FDA0002860665260000052
nx为状态向量的维数；x′为先验状态向量的扰动，且x′＝x-xa；y为根据目标大气数据所确定的观测光谱，
Figure FDA0002860665260000053
ny为观测光谱的维数；Sa为先验协方差矩阵，F(x)表示状态向量x经过正演模型F()模拟得到的模拟光谱，Sε表示观测误差协方差矩阵，γ为修正系数，T表示矩阵的转置；
柱浓度确定模块，用于在所述代价函数迭代后达到收敛标准时，确定此时最优分析场中的二氧化碳廓线分量，并根据所述二氧化碳廓线分量确定二氧化碳柱浓度；
其中，所述处理模块基于非线性最小二乘四维变分同化方法极小化所述代价函数，包括：
根据所述样本状态向量和所述先验状态向量确定初始样本扰动Px，且Px＝(x′1,x′2,...,x′N)；xj′表示第j个初始样本扰动，x′j＝xj-xa(j＝1,2,...,N)，xj表示第j个样本状态向量，N为样本状态向量个数；
将分析增量x′，*转换为初始样本扰动Px的线性组合，即x′，*＝Pxβ，其中β＝(β1,β2,...,βN)T为权重系数向量；
用样本协方差矩阵Be代替先验协方差矩阵Sa，且
Figure FDA0002860665260000054
将代价函数转换为以权重系数向量β为自变量的非线性最小二乘形式：
Figure FDA0002860665260000055
其中，
Figure FDA0002860665260000056
且
Figure FDA0002860665260000057
Sε,+表示经过平方根法分解得到的矩阵，F′(Pxβ)＝F(xa+Pxβ)-F(xa)，y′＝y-F(xa)；
基于高斯牛顿迭代法极小化非线性最小二乘形式的代价函数，且权重系数更新公式为：
Figure FDA0002860665260000061
Figure FDA0002860665260000062
xi,*＝xi-1,*+δxi-1,*；
其中，xi,*为当前迭代的分析场，xi-1,*为上一迭代的分析场，x′,i-1,*＝xi-1,*-xa为上一迭代的分析增量；δβi-1为当前迭代样本扰动
Figure FDA0002860665260000063
的权重系数向量，δxi-1,*表示相较于上一次迭代的分析场的增量，y′,i-1＝y-F(xi-1,*)；
且其中，
Figure FDA0002860665260000064
Figure FDA0002860665260000065
I为单位矩阵，
Figure FDA0002860665260000066
Figure FDA0002860665260000067
Figure FDA0002860665260000068
7.一种计算机存储介质，其特征在于，所述计算机存储介质存储有计算机可执行指令，所述计算机可执行指令用于执行权利要求1-5任意一项所述的二氧化碳柱浓度的反演方法。
8.一种电子设备，其特征在于，包括：
至少一个处理器；以及，
与所述至少一个处理器通信连接的存储器；其中，
所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器能够执行权利要求1-5任意一项所述的二氧化碳柱浓度的反演方法。
CN202010721395.XA 2020-07-24 2020-07-24 二氧化碳柱浓度的反演方法、装置、存储介质及电子设备 Active CN111881569B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010721395.XA CN111881569B (zh) 2020-07-24 2020-07-24 二氧化碳柱浓度的反演方法、装置、存储介质及电子设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010721395.XA CN111881569B (zh) 2020-07-24 2020-07-24 二氧化碳柱浓度的反演方法、装置、存储介质及电子设备

Publications (2)

Publication Number Publication Date
CN111881569A CN111881569A (zh) 2020-11-03
CN111881569B true CN111881569B (zh) 2021-02-05

Family

ID=73200226

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010721395.XA Active CN111881569B (zh) 2020-07-24 2020-07-24 二氧化碳柱浓度的反演方法、装置、存储介质及电子设备

Country Status (1)

Country Link
CN (1) CN111881569B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113204543A (zh) * 2021-04-26 2021-08-03 武汉大学 基于机器学习的二氧化碳柱浓度时空序列调整方法
CN113533241B (zh) * 2021-07-19 2022-12-20 中国科学技术大学 基于卫星红外超光谱的大气二氧化碳浓度高精度反演系统
CN113834902A (zh) * 2021-08-16 2021-12-24 中国人民解放军国防科技大学 一种基于四维变分同化的二氧化硫排放源反演方法
CN114547553B (zh) * 2022-04-27 2022-08-02 河北先河环保科技股份有限公司 二氧化碳排放量的反演方法、装置、设备及存储介质
CN114970184B (zh) * 2022-06-07 2024-04-02 中国科学院地理科学与资源研究所 同步反演高分辨率人为co2排放与自然co2通量的同化方法及系统
CN116153425B (zh) * 2023-02-24 2024-04-02 北京和利时工业软件有限公司 一种氨碳比软测量方法、装置、设备及介质
CN116738232A (zh) * 2023-06-16 2023-09-12 中国科学院空天信息创新研究院 一种基于ftir光谱的城市大气碳排放分布检测方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103983599B (zh) * 2014-02-13 2016-08-17 中国科学院合肥物质科学研究院 直射太阳光谱反演环境大气中二氧化碳垂直柱浓度的方法
US20180066191A1 (en) * 2016-09-06 2018-03-08 Biomass Energy Enhancements Llc Oxygen-Deficient Thermmally Produced Processed Biogas from Beneficiated Organic-Carbon-Containing Feedstock
CN108228538A (zh) * 2017-12-29 2018-06-29 航天科工智慧产业发展有限公司 一种城市中so2浓度值的预测方法
CN108875254B (zh) * 2018-07-03 2023-05-23 南京信息工程大学 一种大气温湿廓线的一维变分反演方法
CN111257241B (zh) * 2020-01-20 2020-11-24 中国科学院地理科学与资源研究所 一种deei的基于卫星观测的大气二氧化碳浓度反演算法

Also Published As

Publication number Publication date
CN111881569A (zh) 2020-11-03

Similar Documents

Publication Publication Date Title
CN111881569B (zh) 二氧化碳柱浓度的反演方法、装置、存储介质及电子设备
Yu et al. Generic atmospheric correction model for interferometric synthetic aperture radar observations
Pail et al. First GOCE gravity field models derived by three different approaches
Scipal et al. A possible solution for the problem of estimating the error structure of global soil moisture data sets
Shi et al. Parameter estimation of a physically based land surface hydrologic model using the ensemble Kalman filter: A synthetic experiment
Bowman et al. Capturing time and vertical variability of tropospheric ozone: A study using TES nadir retrievals
Harnisch et al. Scaling of GNSS radio occultation impact with observation number using an ensemble of data assimilations
Chevallier et al. On the accuracy of the CO2 surface fluxes to be estimated from the GOSAT observations
Scherllin‐Pirscher et al. The power of vertical geolocation of atmospheric profiles from GNSS radio occultation
Poterjoy et al. Efficient assimilation of simulated observations in a high-dimensional geophysical system using a localized particle filter
CN107561557B (zh) 一种掩星探测仪大气成分反演方法
Vogelzang et al. Quadruple Collocation Analysis of In‐Situ, Scatterometer, and NWP Winds
CN111551964B (zh) 一种基于导航卫星系统的水汽观测方法
Lieb et al. Combination of various observation techniques for regional modeling of the gravity field
Zhong et al. A self‐calibration variance‐component model for spatial downscaling of GRACE observations using land surface model outputs
Liu et al. Improvements to a GPS radio occultation ray‐tracing model and their impacts on assimilation of bending angle
Schwarz et al. Integrating uncertainty propagation in GNSS radio occultation retrieval: From bending angle to dry‐air atmospheric profiles
Zhu et al. Estimating climate feedbacks using a neural network
McDonald et al. A new method to evaluate reanalyses using synoptic patterns: An example application in the Ross Sea/Ross Ice Shelf Region
Kostsov Retrieving cloudy atmosphere parameters from RPG-HATPRO radiometer data
Elvidge et al. Multi‐model ensembles for upper atmosphere models
CN107491588B (zh) 一种温压廓线的反演方法和系统
Coopmann et al. Assimilation of IASI ozone‐sensitive channels in preparation for an enhanced coupling between numerical weather prediction and chemistry transport models
Ao et al. Evaluation of CMIP5 upper troposphere and lower stratosphere geopotential height with GPS radio occultation observations
Lomidze et al. Magnetic meridional winds in the thermosphere obtained from Global Assimilation of Ionospheric Measurements (GAIM) model

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant