碳卫星高光谱数据CO2柱浓度反演初步研究

王鼎益，刘冬冬

（1.纽布朗什维克大学 物理系，费德尔顿 E3B5A3；2.西安交通大学 理学院，西安 710049）

大气温室效应引起全球变暖，全球极端天气频发，威胁着人类的生存和发展，已成为世界各国关注的焦点问题。大型水库会对CO2释放和吸收产生影响。三峡水库这类大型水利工程对CO2释放和吸收的影响越来越被重视。三峡库区CO2浓度分布及其时空变化规律是目前的一个研究热点[1-2]。为了准确监测大气中的CO2浓度变化及其对气候的影响，有必要准确地量化全球二氧化碳源和汇的分布和变化。然而，由于温室气体地面观测站点稀疏，对全球碳循环的认识在时间和空间上都非常有限。大气、陆地生态系统和海洋的碳通量测量不准确，特别是在极地、沙漠、海洋和人烟稀少的高原等地区缺乏观测资料。随着技术发展，遥感卫星观测突破了地基单点观测的缺陷，使区域或全球范围探测成为可能。卫星遥感观测正在成为监测全球CO2浓度时空分布与变化最有效的手段[3]。

CO2的人为排放活动主要发生在大气底层。获取大气边界层CO2浓度是分析其产生与消耗（源与汇）机理、监测其时空分布与变化的核心问题[4-5]。研究表明，大气中CO2在近红外波段的2.0 μm和1.6 μm吸收带对近地面层的CO2浓度变化比较敏感[6]。而CO2平均干空气柱浓度摩尔分数（XCO2）能够反映地表附近的CO2浓度的微小变化[7]。为了提高二氧化碳源和汇估计的准确性，需要在区域和全球范围内高精度地解决和量化XCO2中0.25%～2.00%的变化[8]。

目前国际上具有代表性的近红外波段CO2探测仪器包括欧洲的大气制图扫描成像吸收分光计SCIAMACHY （SCanning Imaging Absorption spec⁃troMeter for Atmospheric CHartographY）[9-10]，日本的热和近红外碳观测传感器-傅立叶变换光谱仪TANSO-FTS（Thermal and Near Infrared Sensor for Carbon Observation-Fourier Transform Spectrome⁃ter）[11]和美国的轨道碳观测卫星OCO-2（Orbiting Carbon Observatory）[12-13]。为了提高在全球碳排放和气候变化研究上的国际话语权，中国也开展了卫星探测大气CO2方面的研究工作[14]。在2016年12月成功发射全球CO2监测科学实验卫星碳卫星（TanSat），在2018年5月9日成功发射搭载有大气主要温室气体监测仪的高分五号卫星。这是继日本GOSAT卫星和美国OCO-2卫星之后，世界上第三颗和第四颗用于温室气体监测的卫星[15]。碳卫星上搭载的高光谱大气CO2光栅光谱仪ACGS（At⁃mospheric CO2Grating Spectroradiometer）专门用于探测大气CO2，已经成功获取观测数据，并通过国家卫星气象中心网站面向社会公众开放共享，标志着中国继日本和美国之后成为第三个可以提供碳卫星数据的国家[16-17]。

针对中国全球CO2监测科学实验卫星碳卫星TanSat，中科院大气物理研究所提出了基于最优估计和矢量线性离散坐标传输模型VLIDORT（vec⁃tor linearized discrete ordinate radiative transfer）的CO2反演算法[15，18]，并对碳卫星观测数据进行了反演研究[19]。本文针对TanSat也开展了相应的反演试验研究，且不同于上述算法。本文采用最早针对欧洲空间局环境卫星ENVISAT搭载的大气绘图扫描成像吸收光谱分析仪SCIAMACHY（SCanning Imaging Absorption SpectroMeter for Atmospheric CHartographY）开发的SCIATRAN辐射传输模式和最优化估计算法，结合碳卫星的仪器特性和观测几何，把水汽廓线、温度廓线、地表反照率等作为非需反演的模型参数由其他独立测量或者模型提供，不考虑频率漂移，从碳卫星的观察光谱反演XCO2的信息。本文第一部分概括讨论卫星遥感XCO2的正演模型和反演原理，第二部分简要介绍碳卫星搭载的高光谱大气CO2光栅光谱仪ACGS的特性。在第三部分中展示分析基于SCIATRAN模型和碳卫星数据反演得到的XCO2信息。第四部分给出结论。

1 正演模型和反演原理

碳卫星测量大气温室气体CO2所用的方法是被动遥感测量，观测的目标是CO2分子吸收光谱。当太阳光穿过大气层时，一些光会被气体分子吸收，一些光会被气溶胶和云吸收或散射。每种不同的大气成分气体由于其化学结构不同，将使光在不同的特定波长被吸收。因此，对特定波长吸收光谱的遥感测量可以提供相应气体成分的信息，推断出其浓度[20]。

测量分子的吸收光谱需要专门的仪器（硬件）。分析测量光谱、提取大气成分的相关信息，如浓度、温度、压强以及运动速度等，则需要首先模拟仪器观察的物理过程和光线在大气中的传输过程，产生出仿真模拟的光谱（即正演模型）；然后比较模拟光谱和实际观察光谱，并不断改进模拟光谱，使其与实际观察光谱达到尽可能一致，从而获得大气状态的信息（即反演算法）。以下简单介绍正演模型和反演算法的原理，以及在CO2遥感测量中的应用。

1.1 正演模型

在本文讨论的高光谱卫星遥感中，正演模型是描述一个卫星仪器遥感观察过程的数值模拟模型。正演模型的已知输入参数包括太阳光谱、空间大气、地表的物理和化学性质及其时空变化，以及卫星仪器的机械和遥感性能等。正演的过程就是输入这些已知参数，结合模型，模拟出卫星仪器观察接收的太阳光谱，即根据需要给定系统的参数和状态，预先模拟出系统的结果。

辐射传输模型是正演模型的核心。该模型旨在模拟大气辐射场的全物理—化学过程，其中包括入射太阳光的特性，其在大气中的透射、反射、折射、散射和辐射，以及地面的热辐射和双向表面反射等。理论上讲，当建立一定的观测坐标系后，卫星仪器在大气层顶观察到的辐射强度完全可以由这些参数和边界条件所确定。考虑平行平面介质中一定体积内的能量守恒，可以导出经过散射、吸收和辐射的总Stokes矢量（直接和扩散部分的和）的矢量辐射传输方程[21]为：

式中：Jtot(τ,Ω)和Je(τ,Ω)分别为多次散射源函数和内部辐射的源函数；τ为光学厚度τ∈[0,τ0]，在平行平面介质的顶端为0，底部为τ=τ0；μ=cosϑ∈[-1,1]；Ω为变量{μ,φ}；ϑ是太阳天顶角，ϑ∈[0,π ]；φ为方位角，φ∈[0,2π]；Itot(τ,Ω)是总斯托克斯矢量。

在一定的边界条件下，方程（1）一般只有形式解[21]：

式中：I(τ,±μ,ϕ)表示向下(+μ)与向上(-μ)方向的斯托克斯矢量；J(τ′,±μ,ϕ)表示向下(+μ)与向上(-μ)方向的源函数；I(0,+μ,ϕ)代表上边界条件；I(τ0,-μ,ϕ)代表下边界条件。实际上，方程（1）～（3）是极其复杂的微分—积分方程，是通过计算机模拟计算来实行的。SCIANTRAN传输模型就是广泛使用的正演模型之一，可用于模拟计算平行平面大气和球面大气下的大气辐射传输过程，包括瑞利散射、大气成分的吸收以及地表的双向反射等；可在紫外—可见光—红外光谱范围内（175～4 000 nm），作为地面、机载仪器或星载在天底、临边和掩星观测方式下测量地球大气的散射吸收太阳辐射的正演模型[22]。它还提供了非常丰富的参数化输入接口，用户可根据自己的需要对其进行改造，以完成自己的个性化任务[23]。

1.2 反演原理

卫星遥感观测的最终目的是从卫星仪器观测的光谱或图像中分析和计算出关于大气或者目标的物理及化学性质。为了达到这一目的，就需要有相应的反演算法，正演模型是反演算法的基础。反演算法是利用正演模型产生的光谱或者图像与仪器实际观测的光谱或者图像进行对比，不断调整正演模型的输入参数（大气成分或者目标的物理和化学性质参数），使模拟预测的结果和仪器实际观测的结果尽量达到一致，由此得到一系列大气或者目标的成分和化学性质参数，并将这些输入参数当作实际观测到的大气信息[24]。实际操作中将大气CO2的初始浓度代入正演模型进行计算，得到模拟的大气辐射光谱，将模拟光谱与测量光谱进行比较，如果两者的残差满足反演收敛条件，则迭代计算结束，所输出的大气CO2的初始浓度就是所求的浓度；如果两者的残差超出了反演收敛条件，则调整正演模型的输入参数，使模拟光谱和测量光谱尽量达到一致，再重新进行迭代计算，直到满足反演收敛条件，则经不断调整后的输入参数就是所求的值。

正演模型主要包含辐射传输模型，结合卫星观测模式，考虑太阳辐射和卫星仪器模型函数[24]，可以形式上表示为：

式中：y为测量矢量；x为反演状态矢量；b为非反演状态矢量；ε为误差矢量。F（x，b）是状态矢量的函数即正演模型，描述整个测量的物理过程。

成本函数代表了通过优化测量和模拟光谱之间的差异以及状态向量与先验信息之间的差异所产生的不同最低成本。求解式可通过最优化算法使成本函数最小化[24]。

式中：Sa为先验误差方差矩阵；xa为先验状态矢量；Sε为误差协方差矩阵。

最可能的状态矢量可能与真正的解相去甚远。为了更接近真实的状态，可以使用迭代方法，每次迭代都要重新运行正演模型，但是要使用前次迭代中获得的先验信息。这种先验信息很可能更接近于真正的解，并会导致正演模型产生一个更接近真实的状态矢量。通常使用Levenberg-Mar⁃quardt迭代法求解：

式中：xi+1,xi分别表示第i+1和i次的反演状态矢量；Ki为第i次迭代时的权重函数矩阵；γ为使J(xi+1)最小化而引入的非负常数；D为比例变换矩阵。

1.3 CO2柱浓度 （XCO2）

CO2的柱平均干空气柱浓度摩尔分数（简称CO2的平均柱浓度）是将二氧化碳柱总量用同时从O2-A带反演得到的氧气柱总量归一化后得到的。因为O2分子在空气中的变化十分微小，是一种被广泛认可的、可以准确计算空气柱含量的气体[25]。所以近地面CO2平均柱浓度（干燥空气下）[11]可以表达为：

式中：XCO2表示CO2平均柱浓度（干燥空气下），单位为mg/L；CO2col表示反演的CO2的绝对柱总量，单位为mol/cm2；O2col表示反演的O2绝对柱总量，单位为mol/cm2；O2mf是转换常数，用于将O2的柱含量转化为干燥空气的柱含量，一般取值为0.209 5。CO2绝对柱总量和O2绝对柱总量是分别反演得到的。

2 碳卫星探测仪ACGS

全球CO2监测科学实验卫星碳卫星TanSat上搭载的高光谱大气CO2光栅光谱仪ACGS是由中国科学院长春光学精密机械与物理研究所研制的。采用结构相似的3套光谱仪系统分别对0.76 μm O2分子A带和CO2分子的1.61 μm和2.06 μm波段进行高光谱探测（相关参数如表1所示）[26-27]。由1块Si-CCD探测器（O2-A）和2块MCT探测器（Weak CO2和Strong CO2）分别接收3个波段的气体吸收光谱[26，28]。这3个光谱仪使用相似的光学设计，并被集成到一个共同的结构中以提高系统的刚度和热稳定性。在其焦平面阵列中，面阵横向为光谱维，纵向为空间维。在横向上，O2分子A带的光谱像素个数为1 242，2个CO2波段的光谱像素个数都为500；在纵向上，有9个幅宽采样点。1个面阵数据为一帧数据（frame），即飞行方向某一时间采样的观测数据。数据帧数（飞行方向数据长度）与观测时间相关。

表1 碳卫星大气CO2遥感探测仪的相关参数Tab.1 Relevant parameters of TanSat ACGS

TanSat卫星主要有3种观测模式，分别是天底模式、耀斑模式和目标模式。探测仪器的视线指向当地的最低点（即天底观测模式，Nadir observa⁃tion）或者是闪烁的光点（即耀斑观测模式，Glint observation），还可以瞄准选定的地球表面校准和验证点（即目标观测模式，Target observation）。在这些地方，阳光从地球表面反射出来经过大气进入探测仪[29]（如图1所示）。

Nadir观测模式提供了最佳的水平空间分辨率，并有望在部分多云地区或地形上产生更多有用的XCO2探测。Glint观测模式在黑暗、镜面表面有比较大的信噪比，预计在海洋上会产生更有用的探测结果。通常，碳卫星在Nadir观测模式和Glint观测模式之间交替进行。Target观测是在碳卫星验证点上进行的，并收集成千上万的观测数据，大量的测量减少了随机误差的影响，并提供了识别目标附近XCO2场空间变异性的信息。

目前，碳卫星已经对外共享了经过定标后的L1B光谱数据集，所有产品文件都是以层次型科学数据格式HDF-5发布。这种格式有助于创建逻辑数据结构，通过将数据产品组织到文件夹和子文件夹中，每个文件对应一个轨道连续模式的数据集。

图1 天底（a）、耀斑（b）和目标（c）观测模式[29]Fig.1 Nadir(a),glint(b)and target(c)observation mode

光谱辐射包含在L1B产品的光谱测量文件夹中，其中包括各种类型，在每个类型中又包含多种变量，每个变量可以是单值的或多值的。每个变量的元数据描述变量的属性，如维度、数据解释和单位等。碳卫星L1B产品中包含了用于反演XCO2产品的测量校准光谱，如图2所示。其中图（a）、图（b）、图（c）分别为波段1（758～778 nm）、波段2（1 594～1 624 nm）、波段3（2 041～2 081 nm）定标后的光谱图。此外，还包含了所有测量的观测几何信息（太阳天顶角、视线极轴角和视线方位角等），卫星轨道信息（轨道高度和飞行姿态等），观测点的地理位置信息（日期、经度、纬度等），以及相应的仪器和校准信息。

图2 TanSat卫星近红外测量光谱示例Fig.2 Sample spectra of TanSat near infrared measurements

3 反演试验

本节针对碳卫星TanSat探测仪ACGS的仪器特性和观察几何特征，基于SCIATRAN辐射传输模式和最优化算法，利用碳卫星L1B光谱数据，对碳卫星TanSat观测的XCO2浓度进行了反演试验研究。

在三种观测模式下分别选择了陆地（北京）、海洋（印度洋）和冰川（南极大陆）为研究对象。以TanSat卫星实际测量的L1B光谱数据作为输入光谱来反演XCO2。计算中，XCO2的先验初始浓度值设定为405 mg/L；采用SCIATRAN模型模拟大气中的辐射传输过程，并将其与TanSat仪器特性和观测场景相结合。温度廓线采用欧洲中期天气预报中心ECMWF的ERA-Interim数据库与美国标准大气模型US76的融合模式，海拔范围是0～60 km。地球表面被假定为朗伯体，地表反照率采用中分辨率成像光谱仪MODIS的MCD43C3数据，当观测点的经度和纬度位于0.1°×0.1°的数据网格中时，地表反照率是四个顶点的线性插值。全球高程数据采用GTOPO30数据。大气痕量气体廓线的初始输入数据采用不莱梅大学B2D二维大气模型作为输入数据库，并在涉及高度廓线的数据计算中会根据大气分层模式对廓线插值。HITRAN2012大气分子吸收谱数据库被用作分子光谱参数数据库[30]。气体吸收光谱由逐线积分模型计算，积分步长为0.005 nm。气溶胶数据来自LOWTRAN 7模型，该模型可以设置不同的气溶胶类型（海洋型、城市型和乡村型等），其中气溶胶层的数量和高度是可选择的，可根据类型选择不同的能见度和相对湿度（0，70%，80%和99%）。入射的太阳辐照度从Kurucz发布的太阳吸收线表获得[31]。采用离散坐标法计算权重函数矩阵。

图3所示为目标观测模式下北京地区（2017年5月27日）的TanSat卫星数据XCO2反演结果，从图中可以看出，北京地区XCO2浓度分布不均匀，城市中心区域的XCO2浓度显著高于边缘区域，且随着远离城市中心XCO2浓度是逐渐降低的。北京城北城南XCO2高浓度区域明显大于城东城西区域，而且城区西北边缘的山区形成明显的低值，这些都与北京的城市格局相一致，反映了人类活动频繁地区XCO2浓度高的现象和对环境的影响。北京城区 （39°～40.5°N，116°～117°E） 的XCO2浓度平均值为405.9 mg/L，而国家卫星气象中心观测的当天北京城区的XCO2浓度平均值为406.9 mg/L（私人通信，国家卫星气象中心杨忠东），两者仅相差1 mg/L。

图3 2017年5月27日目标观测模式下TanSat数据反演XCO2浓度分布图（北京）Fig.3 XCO2distribution of target observation retrieved from TanSat data（2017.05.27 Beijing）

图4和图5所示分别为耀斑观测模式下海洋区域（2017年5月10日）和天底观测模式下冰川区域（2017年3月12日）TanSat卫星数据XCO2反演结果。从图中可以看出，海洋区域和冰川区域XCO2浓度呈现相对均匀的分布状态，这是因为在这些地区人类活动较少，对环境的影响较小。上述反演结果与GOSAT[32]和OCO-2[33]发布的XCO2浓度全球分布在对应地区是定性一致的，其浓度分布相对均匀，浓度集中在390～410 mg/L。针对不同观测模式和不同下表面条件的研究结果证明了本文建议的反演算法应用于碳卫星观测数据的可行性。

图4 2017年5月10日耀斑观测模式下TanSat数据反演XCO2浓度分布图（印度洋）Fig.4 XCO2distribution of glint observation retrieved from TanSat data（2017.05.10 Indian Ocean）

图5 2017年3月12日天底观测模式下TanSat数据反演XCO2浓度分布图（南极冰川）Fig.5 XCO2distribution of nadir observation retrieved from TanSat data（2017.03.12 South Pole）

4 结论

本研究针对碳卫星TanSat探测仪ACGS的仪器特性和观察几何特征，基于SCIATRAN辐射传输模式和最优化算法，以碳卫星三种观测模式（天底、耀斑和目标模式）下的陆地（北京）、海洋（印度洋）和冰川（南极大陆）为研究对象，利用碳卫星L1B光谱数据，对碳卫星TanSat观测的XCO2浓度进行了初步反演试验研究。在三种观测模式和不同下表面条件下的反演结果与国家卫星气象中心、GOSAT和OCO-2的结果是基本一致的，表明本文所采用的反演算法对于处理碳卫星数据是可行的。但本研究只是针对少量数据的初步研究结果，还需要针对全球不同地区和场景的观测数据做大量和进一步的反演试验研究，并对反演算法做相应的改进，最终的研究成果将为在中国领土上对温室气体CO2浓度变化的全方位、高精度监测提供技术支持。