具体实施方式

为获得小尺度(0.1°或10km)的城市CO₂格网化通量，本发明提出利用拉格朗日粒子扩散模型计算指定三维时空位置(日期、时间、经度、纬度和高度)的逆向足迹，进而依靠已有的人为CO₂排放通量和生物CO₂通量(NEE)构建先验CO₂格网化通量场及其时空协方差，依靠上述两种数据集建立通量场与浓度场之间的定量关系，并利用贝叶斯反演获得经优化的后验CO₂格网化通量产品。

参见图1，本发明实施例提供一种基于星载激光雷达获取城市CO2格网化通量的方法，包括以下步骤：

步骤S1,卫星星下点位置的精准结算：根据需要获得CO₂通量的空间范围和时间段，利用高精度卫星轨道模拟器计算卫星星下点位置，包括经纬度和过境时间。

实施例中，本步骤使用基于python语言的ephem库和北美防空司令部所定义的二行卫星轨道参数计算卫星轨道信息，并利用自定义的时空范围信息获得指定地区和时间范围内的卫星星下点坐标信息。为便于实施参考起见，提供具体实现过程如下：

首先，利用ephem库的readtle函数和搭载了CO₂-IPDA LIDAR的卫星二行轨道参数，完成初始化。

其次，利用python的datetime数据类型设置起始时间，同时利用经纬度限定待计算的空间范围。

第三，依靠起始时间和持续时长建立的时间向量驱动ephem库的compute函数，得到指定时间范围内的全部星下点经纬度，并判定其是否位于指定的空间范围内。

第四，输出符合要求的卫星星下点时空信息，其具体形式为属性包括日期、时间、经度和纬度的数据表格。

最后，根据经度和纬度信息，利用高精度的数字高程模型为每一条数据进一步加上高程属性。至此，最终得到指定时空范围内的卫星星下点信息。

步骤S2,分层足迹场获得。

本步骤优选使用WRF-STILT模型模拟在气象场控制下指定位置(经度、维度、高度)的粒子的后向随机轨迹，并得到该点的足迹场。WRF-STILT模型为拉格日朗大气传播模型。实施例将卫星星下点的时空信息作为输入，利用WRF驱动的STILT模型计算每个星下点处0-15km的分层足迹。

本发明的目标在于获得0.1°(～10km)分辨的通量信息，因此足迹场的分辨率需要设置为0.1°。因此，气象场分辨率需要优于0.1°才行。WRF的设置需要利用三层嵌套，其具体方案为(27km，9km，3km)或者(32km，8km，2km)。WRF模拟的时间范围需要根据星下点的时间进行倒推，推荐的倒推时长为7天，该时长最短不得低于3天，最长不宜超过15天。WRF其他设置(大气方案、辐射方案、地形方案等)则可根据指定地区的情况和时间进行确定。将WRF输出的气象场进行格式转换后输入STILT模型，计算指定位置的足迹场。

STILT模型是一种拉格朗日随机游走理论的传输模型，它把观测点上游的源(汇)通量与观测点的浓度变化用足迹权重联系起来。STILT在城市地区能够精确模拟超高分辨率(<1km)大气传输过程的能力，因此该模型完全满足城市CO₂通量反演对于0.1°分辨足迹的需求。其具体原理就是通过向后释放大量的空气粒子来模拟气体在湍流和平均风向驱动下的后向运动轨迹，通过计算上游某区域边界层某高度内的所有粒子数量和每个粒子所停留时间来定量计算足迹权重的值f()，其计算公式如下：

其中，x_r为观测点所在的位置和高度，t_r为对应的观测时间，m_air为空气的摩尔重量，h为下垫面影响层的高度，ρ()为所有粒子的平均密度，t_m代表一个离散的时间步长，(x_i,y_j)代表某一区域的下垫面，p代表一个粒子，Δt_p,i,j代表每一个粒子对应在某一区域(x_i,y_j)的下垫面影响层所停留的时间。N_tot为释放的粒子总数。

由于卫星观测对象是地表至大气层顶的整层CO₂柱浓度，理论上对于各个观测点进行计算其不同高层的分层足迹。这一过程对计算资料消耗极大，本步骤特此给出一种经过多次实验获得的优选分层方案。具体为：

(1)0～0.8km范围内以0.2km进行分层，计算4个分层的足迹场；

(2)0.8～2.0km范围内以0.4km进行分层，计算3个分层的足迹场；

(3)2.0～5.0km范围内以1km进行分层，计算3个分层的足迹场；

(4)5～12km范围内计算10km和12km两处的足迹场。

(5)12km以上不再计算，最终对于每个观测点得到12层的分层足迹场。

步骤S3,先验通量场及协方差构建。

实施例利用现有的人为CO₂排放清单和自然CO₂通量(NEE)产品，在0.1°(～10km)的空间尺度上构建CO₂先验通量场，并设定先验协方差。

先验通量场由人为和自然CO₂通量场组成而生成。非美国地区的人为CO₂通量场的构建推荐利用EDGAR5.0以上版本的排放清单，美国地区推荐使用VULCAN。自然CO₂通量场推荐使用基于SIF(太阳诱导荧光)的NEE产品集，建议使用开源的SMUrF代码自行计算指定地区的高分辨NEE，该模型具备在城市地区较高的精度，因此适用于本发明所需。

先验协方差矩阵B(m×m维)可以用下式表示:

其中，D(d×d维)是时间协方差矩阵和E(e×e维)是空间协方差矩阵，是Kronecker乘积，所以有n＝d×e。

时间协方差矩阵和空间协方差矩阵可以用下式表示：

D＝V_t ^1/2M_tV_t ^1/2 3

E＝V_s ^1/2M_sV_s ^1/2 4

其中，M_t和M_s分别是时间相关矩阵和空间相关矩阵，V_t和V_s分别是时间方差和空间方差的对角矩阵。

步骤S4,权重函数计算。

实施例利用气象资料提供的大气温度、湿度和压力垂直廓线和HITRAN数据库提供的CO₂分子吸收参数，采取Voigt线型计算每个卫星CO₂浓度产品的权重函数。

公式5给出了利用IPDA(路径积分差分吸收)反演大气CO₂柱浓度的理论方程，其中P_λ和E_λ为波长为λ时的回波信号强度和发射能量，ω(p)是CO₂浓度廓线。分母部分即大气权重函数(见公式6)。公式6中σ为分子吸收截面积，M是分子质量，ρ是密度，p是大气压力，g是重力加速度。IPDA LIDAR的2个工作波长通常记为λ_on、λ_off，因此λ＝λ_on或λ_off。

从公式6可以看出，WF中的变量包括了σ和ρ_H2O。同时需要指出的是，尽管对于大气探测一般以σ(p,λ)的形式表达，但是在指定特定气体分子后，σ()实际上是一个由温度、压力和波长共同决定的变量。具体表达式由公式7-11给出，其中， ν是波数(波长的倒数)，ν_c是吸收谱的中心波数，S(T)是吸收线强度，t是积分变量，γ_L是洛仑兹半宽，γ_D是多普勒半宽。p是某一大气压力，p₀和T₀分别为标准大气压和标准温度，为101325Pa和296.15K。γ₀是标准线宽度，n是压力增宽系数的温度依赖系数，k_B是玻尔兹曼常量，E”是低能态，h是普朗克常量，c是光速。

ν_c＝ν₀+p×(δ⁰(T₀)+δ'(T-T₀)) 10

其中，ν₀是大气压力p＝0Pa,温度T＝T₀时的吸收峰波数，δ⁰(T₀)是T＝T₀时的压力平移系数，δ'(T-T₀)是压力平移系数的温度依赖系数，T₀为标准温度；ν_c是吸收谱的中心波数；

M是分子质量，c是光速，k_B是玻尔兹曼常量，γ_D是多普勒半宽；

γ_L是洛仑兹半宽，γ₀是标准线宽度，p₀为标准大气压，n是压力增宽系数的温度依赖系数；

S(T)是吸收线强度，S₀是温度为296.13K，压力为101325Pa时的吸收线强度，E”是低能态，h是普朗克常量，exp()是指数函数；

σ(p,T,ν)是温度T、压力p和波数ν下的分子吸收截面积， ν是波数(波长的倒数)；

WF(p)是权重函数，其描述了不同压力处特定大气分子对光的吸收能力，σ(p,λ_on)是压力p、波长λ_on下的分子吸收截面积，σ(p,λ_off)是压力p、波长λ_off下的分子吸收截面积，是水分子质量，M_air是空气分子量，是水汽密度，g是是重力加速度；

XCO₂是柱均CO₂干空气混合比，Psurf是地表压强，ω(p)是CO₂浓度廓线，P_on是波长为λ时的回波信号强度，P_off是波长为λ时的回波信号强度，E_on是波长为λ_on时的回波信号强度，E_off是波长为λ_off时的回波信号强度。

图2给出了由S4计算出的WF分布，该WF是依据目前大气环境星参数计算给出，可以直接使用。若将来由其他类似激光CO₂探测卫星，需要按照S4的步骤重新计算。

步骤S5,柱加权足迹场的生成。

实施例将步骤4得到的权重函数和步骤2获得的分层足迹进行向量乘法，得到柱加权足迹场。

将步骤S2所得到的分层足迹，按照对应高程乘以步骤S4所得到的WF，并对N个乘积求和，从而得到柱加权足迹场。该过程可以由公式12表达。z表示不同的分层序号，步骤2中推荐了12层(N＝12)的分层方案。除此之外，也可以根据需求设计其他分层方案。

其中，footprint(z)表示高程z处的地表通量足迹。

步骤S6,通量、足迹和浓度关系的目标方程构建。

实施例中，利用步骤3得到的先验通量信息，步骤5得到的柱加权足迹场和星载CO₂-IPDA LIDAR观测的CO₂柱浓度构建目标方程。

x^*表示CO₂通量的先验场，其与CO₂浓度的关系可以由下式给出：

y^*＝Hx^*+ε. 12

其中，y^*是XCO₂的观测值，ε是均值为μ_b且对角线协方差矩阵为R的N×1向量，其是服从ε～N(μ_b,R)正态分布的噪声。R由观测设备本身性质和观测环境共同决定，在实践中可以用指定时间范围内观测结果的波动水平来表征。H是WRF-STILT输出的足迹矩阵，其物理意义为不同位置碳源汇对指定地点CO₂浓度值的贡献。假设失配误差不相关且平均偏差为零。则可以在贝叶斯推理框架中使用这些综合观测来估计最佳CO₂通量。

步骤S7,贝叶斯反演通量场的最优解。

实施例中，使用贝叶斯反演技术优化步骤6给定的目标方程，得到经过优化的0.1°分辨CO₂通量后验解。

假设先验误差的概率分布服从高斯分布，则可给出后验CO₂通量的闭式解：

其中，x_p是先验CO₂通量，B是m×m先验误差协方差矩阵，由步骤S3给出。为本发明最终得到的高分辨格网化CO₂通量产品。

具体实施时，本发明技术方案提出的方法可由本领域技术人员采用计算机软件技术实现自动运行流程，实现方法的系统装置例如存储本发明技术方案相应计算机程序的计算机可读存储介质以及包括运行相应计算机程序的计算机设备，也应当在本发明的保护范围内。

在一些可能的实施例中，提供一种基于星载激光雷达获取城市CO2格网化通量系统，包括以下模块，

第一模块，用于根据需要获得CO₂通量的空间范围和时间段，利用高精度卫星轨道模拟器计算卫星星下点位置，包括经纬度和过境时间；

第二模块，用于将卫星星下点的时空信息作为输入，利用WRF驱动的STILT模型计算每个星下点处的分层足迹；

第三模块，用于利用现有的人为CO₂排放清单和自然CO₂通量产品，在预设空间尺度上构建CO₂先验通量场，并设定先验协方差；

第四模块，用于利用气象资料提供的大气温度、湿度和压力垂直廓线和HITRAN数据库提供的CO₂分子吸收参数，采取Voigt线型计算每个卫星CO₂浓度产品的权重函数；

第五模块，用于将第四模块得到的权重函数和第二模块获得的分层足迹进行向量乘法，得到柱加权足迹场；

第六模块，用于利用先验通量信息、柱加权足迹场和星载CO₂-IPDA LIDAR观测的CO₂柱浓度构建目标方程；

第七模块，用于使用贝叶斯反演优化第六模块给定的目标方程，得到经过优化的预设空间尺度CO₂通量后验解。

在一些可能的实施例中，提供一种基于星载激光雷达获取城市CO2格网化通量系统，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如上所述的一种基于星载激光雷达获取城市CO2格网化通量方法。

在一些可能的实施例中，提供一种基于星载激光雷达获取城市CO2格网化通量系统，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如上所述的一种基于星载激光雷达获取城市CO2格网化通量方法。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。