气溶胶细粒子比偏振遥感最优化反演方法研究

郑逢勋， 李正强， 侯伟真， 董晓刚， 周志远

1. 河南科技大学 农业装备工程学院, 洛阳 471003;

2. 中国科学院空天信息创新研究院 国家环境保护卫星遥感重点实验室, 北京 100094;

3. 北京环境特性研究所 光学辐射重点实验室, 北京 100854

1 引 言

人为活动向大气中排放的细颗粒物会对人体健康造成严重影响（Hoek 等，2013；Cohen 等，2017；Song 等，2017）。描述气溶胶细粒子占比的参数FMF（Fine Mode Fraction）是近地面PM2.5纯物理遥感估算模型中的关键参数之一（Li 等，2016）。FMF 遥感研究对于理解人为因素对气候和环境的影响和大气颗粒物对人类健康的影响具有重要意义（Yan 等，2017a）。

目前FMF 数据产品主要有：（1）AERONET（AErosol RObotic NETwork）地基FMF 产品，利用光谱退卷积算法SDA （Spectral Deconvolution Algorithm），由气溶胶光学厚度AOD （Aerosol Optical Depth）直接计算得到（O’Neill 等，2003）。地基FMF 产品具有较高精度，然而覆盖度低，主要用于卫星反演产品的地面验证（Li 等，2020；Che 等，2018）。（2）MODIS （Moderate Resolution Imaging Spectroradiometer）FMF 产品，基于暗目标方法，反演不确定性较大（Levy 等，2010），精度上较难满足PM2.5估算需求（Zhang 和Li，2015）。Zhao 等（2017）尝试融合MODIS 和AERONET 的FMF 产 品 以支撑PM2.5遥感估 算。（3） POLDER（POLarization and Directionality of the Earth’s Reflectances） FMF 产 品， 采 用GRASP 算 法（Román 等，2018；Torres 等，2017）反演得到。中国区域验证结果表明该产品与地基观测具有较好的一致性（Wei 等，2020），是目前精度较高的FMF 产品。此外，Yan 等（2017b，2019）将SDA方法应用于MODIS 的AOD 产品，反演的FMF 在精度上比官方产品有所提升；Zhang 等（2016a）基于POLDER 偏振和强度数据分别反演细粒子AOD和总AOD，进而得到FMF。总体上，现有FMF 数据产品和反演研究较少，且受到传感器观测信息量的限制，精度难以满足实际应用需求。

继POLDER 传感器后，具备偏振能力的载荷成为气溶胶遥感领域的发展方向（Dubovik 等，2019）。中国高分5 号卫星搭载的DPC（Directional Polarimetric Camera）（Gu 和Tong，2015）、TanSat卫星搭载的CAPI （Cloud and Aerosol Polarimetric Imager）（Chen 等，2017a） 和HJ-2 卫星搭载的PSAC （Polarized Scanning Atmospheric Corrector）（Hou 等，2019），都具备偏振测量能力。偏振探测相比传统强度测量能够提供更多的大气信息，尤其是长波近红外波段如1610 nm 的偏振测量，对于粗模态气溶胶反演有重要提升价值（Qie 等，2015），为精确反演FMF 提供了可能。现有反演算法中，SDA 算法属于半经验方法，MODIS 的暗目标 法 利 用 查 找 表LUT （Look Up Table） 反 演，GRASP 算法基于最优估计OE（Optimal Estimation）理论（Rodgers，2000），由AERONET 反演算法（Dubovik 和King，2000）发展而来，属于纯物理的全反演方法。LUT 反演效率高，反演时需要进行插值，精度受到一定限制（Kim 等，2018），OE方法具有比LUT 方法更高的理论反演精度，缺点是计算复杂，效率较低（Diner 等，2012；Jeong等，2016）。由于OE 算法的复杂性，国内相关研究仍然以理论分析为主（Chen 等，2017b；Zheng等，2019），基于OE 反演框架和长波近红外偏振测量的FMF反演研究有待深入开展。

本文利用长波近红外偏振测量包含的气溶胶信息量优势，基于现有OE反演框架（Li 等，2018a），提出了一种联合强度和线偏振度DOLP（Degree of Linear Polarization）测量的FMF 参数反演方法。由于目前近红外偏振波段反演尚缺乏实际数据支撑，因此本研究使用矢量化的辐射传输模型针对地基天顶观测模式开展偏振观测仿真，分析偏振数据对FMF参数的波段敏感性，通过对仿真数据叠加观测误差对算法进行测试，验证反演算法的可行性和一致性，为后续实际数据反演研究提供参考。

2 模型与方法

2.1 前向建模

线偏振度（DOLP）是重要的偏振物理量。相对于强度测量，DOLP 作为“相对测量”量，在观测方面天然具有精度高等优势（Yan 等，2018）。目前，强度测量的不确定性在5%左右，而在440—1640 nm 波长范围内，DOLP 能够达到0.005 的绝对误差（Li 等，2018b）。作为描述大气偏振特性的重要参数，DOLP 对粒子尺寸和形状非常敏感，能够提高气溶胶微物理参数反演精度（Li 等，2009；Chen 等，2018）。DOLP 可以通过Stokes 矢量的前3个分量计算：

式中，I为强度分量，Q和U为偏振分量。

对于线偏振度测量的前向建模，本研究采用线性化的矢量辐射传输模式UNL-VRTM（UNified Linearized Vector Radiative Transfer Model）进行辐射传输计算（Wang 等，2014）。UNL-VRTM 能够模拟给定观测几何和光谱的辐射和偏振量（Stokes矢量），并同时输出Stokes 矢量对于气溶胶参数的敏感度矩阵（Jacobian matrix），非常有利于OE 反演算法的实现（Hou 等，2016 和2017），且成功应用于气溶胶遥感研究（Hou 等，2020；Xu 等，2015）。

对于下行辐射建模，大气散射的变化主要依赖于气溶胶粒子的物理和光学特性，地表强度反射率的不确定性影响较小（Dubovik 和King，2000）。研究表明，DOLP 为0.2 以及单次散射反照率为0.85 情况下，地表对下行偏振基本无影响（Li 等，2006）。因此，在前向模拟中本文采用了朗伯地表假设，并忽略地表偏振对测量观测的贡献，对气溶胶的垂直分布，本研究采用高斯分布模型，标高设置为2 km。

本研究采用双模态对数函数分布形式描述气溶胶粒子谱分布情况：

式中，V0、rv和σg分别为气溶胶的体积柱浓度、体积中值半径和标准方差，上角标i 表示与粒子尺寸相关的气溶胶模态；进一步，气溶胶光学厚度τ可以通过上述物理量计算得到（Xu 等，2013）：

式中，Qext为气溶胶消光效率因子，有效半径reff、有效方差νeff与体积参数之间通过下式进行转换：

以f 和c 分别表征细模态（i=1）和粗模态（i=2），则气溶胶光学厚度τ和细粒子比FMF 可以表示为：

2.2 最优估计反演方法

根据OE理论，复杂观测系统可以简单描述为：

式中，y、x分别表示观测向量、状态向量；F为大气辐射传输过程；ϵ是经验误差项。本研究中，y由强度和线偏振度测量构成，x由AOD 和FMF 参数构成。反演本质上是对上式进行反向求解，在不考虑模型误差情况下，反演结果对应非线性代价函数J（x）的最小值。J（x）的形式为

式中，Sy为观测误差协方差矩阵；xa和Sa分别为x的先验估计和误差协方差矩阵，描述了测量前对状态量的预估及其误差范围；上角标“-1”和“T”分别代表矩阵求逆和转置操作；γ为正则化系数（Xu等，2015）。上式右侧分别描述了观测误差和先验误差对反演的影响。代价函数J(x)是一个非线性函数，其梯度向量形式为：

式中，K表示前向模型对状态量的敏感度矩阵（Jacobians Matrix），其中Ki，j=∂Fi/∂xj。本研究采用拟牛顿算法程序包L-BFGS-B 对J(x)的最小值进行求解。通过将L-BFGS-B程序模块、前向辐射传输模拟模块和代价函数的计算模块相结合，进而建立最优估计迭代反演的程序框架。在程序框架下，经过多次迭代计算寻找到状态向量的最优估计，迭代结束的收敛条件为：

式中，‖ ‖2为L2范数，δ为判断迭代结束的收敛阈值。最优估计反演流程如图1所示。

图1 最优估计反演流程Fig.1 The flowchart of the optimal estimation inversion

2.3 模型及算法参数设置

2.3.1 气溶胶模型参数

表1 给出了研究采用的气溶胶模型参数。其中，气溶胶粒子谱分布参数（有效半径reff和有效方差veff）来自北京地区AERONET站点多年观测数据的统计。复折射指数实部mr和虚部mi来自于Zhang 等（2016b）建立的双模态气溶胶模型数据。由于FMF 表征了不同粒子混合的情况，因此气溶胶类型可以通过调节FMF 实现。例如FMF 大于0.8表示细模态主导类型，如城市污染型气溶胶，FMF 小于0.5 表征了粗粒子主导气溶胶类型，如沙尘型气溶胶。

表1 研究采用的气溶胶模型参数Table 1 Parameter settings of aerosol model

2.3.2 观测几何设置

表2 给出了前向模拟时采用的观测几何设置，太阳天顶角为60°，观测天顶角为0°（天顶观测模式），传感器位于太阳主平面内。建立前向仿真数据时考虑了不同的气溶胶负载和不同的细粒子占比情况。AOD 从0.1 到3.0 变化，表征了不同气溶胶含量的情况，0.1 到1 之间AOD 间隔0.1，1 到2之间间隔0.2，0.2 到3 之间间隔0.4。FMF 参数从0.05到0.99变化，步长0.05。

表2 前向模拟观测几何及气溶胶设置Table 2 Observation geometry and aerosol settings in the forward simulation

2.3.3 观测量和反演量设置

观测向量y由总辐射强度和DOLP 两种具有不同测量精度的测量组成，其中辐射强度测量包含了5 个波段，DOLP 测量包含4 个波段，如表3 所示。在误差设置方面，对于强度测量，观测不确定性εI设为5%，对于星载偏振测量，偏振度0.4情况下，绝对测量精度可保障在0.005（杨洪春等，2019），因此，本研究DOLP 相对误差εDOLP设为1%。根据式（7）可知，在最优估计反演过程中，先验估计的不确定性会影响到反演结果。考虑到气溶胶模型获取的不确定性，先验估计误差设定为100%（Hou 等，2017）。观测量、反演量和误差协方差矩阵的组成如表3 所示，其中，diag(I)表示构建矩阵操作，以向量I的元素作为矩阵的对角元素，矩阵非对角元素填充零。

表3 观测向量、状态向量及相应的误差定义Table 3 Definition of observation vector， state vector and prior estimation error

3 结果与分析

3.1 仿真结果分析

基于表1 和表2 中的参数设置，使用UNLVRTM 模式进行地基天顶观测模式的仿真。图2 给出了天顶观测模式下不同气溶胶光学厚度和不同的粗、细气溶胶粒子混合比下天空光辐射测量（图2（a））和偏振测量（图2（b））的模拟结果。横坐标轴为AOD，纵坐标轴为FMF。由图2可知：

（1）对于强度测量：天空光辐射强度随着气溶胶光学厚度AOD 的增加而增强，在AOD 为2.5时达到最大值，然后逐渐减小。在相同AOD 情况下，天空光的强度随着细粒子含量的增加而增强。低AOD 情况下，强度测量受FMF 变化的影响较小。

（2）对于线偏振度测量：由图2（b）可知，在490 nm波段，DOLP随气溶胶负载的增加而逐渐减小，但不同FMF取值下，DOLP 变化不明显，这表明该波段观测量的变化主要由AOD的变化引起，而对FMF 参数变化不敏感；在670 nm、870 nm 和1610 nm 波段，相同AOD 情况下，DOLP 随FMF 的增加出现先减小后增大的情况。DOLP 的最小值与FMF 大小有关，670 nm、870 nm 和1610 nm 这3 个波段的DOLP 的最小值分别出现在FMF 为0.1、0.2和0.7附近。表明随着波长的增大，DOLP最小值所对应的FMF 逐渐增大。这是由于在近红外870 nm和长波近红外1610 nm 波段，粗粒子气溶胶的偏振贡献相对于可见光（490 nm 和670 nm）有所增强，即近红外波段的观测对粗模态气溶胶更为敏感；在1610 nm 波段，DOLP 的变化主要由FMF 的变化引起，与490 nm 存在显著差异，表明该波段对FMF参数的变化比对AOD的变化更敏感。

图2 天空光辐射强度和线偏振度测量仿真（天顶观测模式，太阳天顶角60°）Fig. 2 Simulation of intensity and degree of linear polarization （DOLP） for sky light measurements （vertical upward observation with the sun zenith angle equaling 60°）

（3）由强度和偏振结果对比可知：辐射强度随波长的增加而减小，而DOLP 的波段间差异并没有这一趋势。尤其对于细模态主导的气溶胶类型（FMF>0.8），随着气溶胶光学厚度的增加，近红外波段的DOLP显著高于可见光波段。

3.2 反演结果分析

基于模拟观测数据，对上述OE 反演算法进行测试。研究采用了3 种不同的情形，测试说明如表4 所示。S1 为未引入误差情况下的强度和偏振联合反演，用于测试算法的自我一致性；S2 和S3采用服从高斯分布的随机噪声对模拟观测值施加误差扰动，用以模拟真实观测，其中，S2 为纯强度反演情形，S3 为强度偏振联合反演情形。S2 和S3的差异用于对比分析偏振测量对反演的贡献。

表4 反演测试说明Table 4 Data scenarios adopted for inversion test

图3 给出了模拟数据为S3 条件下，AOD 和FMF 均为0.6 时，最优化反演过程中辐射强度和线偏振度随迭代次数增加的变化情况。由图3 可知，仿真模拟值随着迭代次数增加逐渐逼近观测数据。研究采用反演残差描述反演结果对应的模拟值与观测值的差异，反演残差定义为

式中，yfit表示最后一次迭代时的模拟值，y为观测值。图3 中第15 次迭代时的反演残差仅为0.6%，S3 条件下全体模拟数据的平均反演残差为5.2%，表明最优估计反演算法能够根据参数敏感度矩阵和最小值求解过程，实现对反演参数的动态调整，并最终在前向光辐射模拟上与观测值达成很好的一致性。

图3 最优估计反演的迭代过程Fig. 3 The iterative process of the optimal estimation inversion

图4 给出了550 nm 波段下不同气溶胶光学厚度和细粒子比条件下的反演结果和真实值之间的散点图。其中图4（a）和（d）为没有观测噪声引入情况下的反演结果，反演结果和真值的相关系数大于0.99，由算法本身引起的AOD 和FMF 反演误差分别为0.085%和0.014%。由此可知，本算法具有很好的闭合性。图4（b）和（e）为仅采用强度观测数据时的反演结果，图4（c）和（f）为强度和偏振联合反演的结果，通过对比可知，在测量误差存在的情况下，引入DOLP 测量对于气溶胶反演精度的提高非常显著。表5给出了反演误差统计结果，增加DOLP 数据参与反演的情况下，AOD反演误差从1%下降到了0.3%，FMF参数反演误差从1.4%下降到了0.18%，反演精度相对于强度反演有明显的提升。

表5 3种情况下AOD和FMF的平均反演误差Table 5 Average inversion error of AOD and FMF/%

4 讨 论

根据代价函数和梯度向量的表达式可知，OE反演受到观测和先验估计的约束，反演结果的精度不仅取决于观测数据的不确定性，同时受到先验知识的影响（Xu 和Wang，2015）。获取反演参数的先验知识的是最优化反演算法的重要组成，基于先验知识的平滑约束有助于参数反演的求解（Fedarenka 等，2016）。先验知识来源于对反演参数的实验室测量、气象数据统计以及其他途径的遥感观测。合理的先验估计能够约束反演结果的合理性，减少算法的迭代次数（Dubovik 等，2011）。

最优估计反演充分利用强度和偏振数据中的有效信息，通过在连续的向量空间中搜索最优解，理论上具有比传统查找表方法更高的反演精度（Jeong 等，2016）。然而，由于反演过程中，需要进行多次迭代和实时前向辐射传输计算，因此算法具有较高的耗时。在反演测试中，平均迭代次数在10次左右。未来将OE方法应用于大数据量的反演需要进一步优化算法的计算效率。

目前缺乏长波近红外偏振数据支持，因此本研究基于模拟数据开展了算法测试。采用模拟数据和真实数据反演的主要差异在于前者不存在强度和偏振数据自洽性问题。尽管研究对模拟数据叠加了噪声，但实测数据的不确定性和误差要比模拟数据更复杂，后续研究需要针对实测数据进行参数优化。此外，测试针对单角度观测的仿真数据开展，而多角度观测能够显著降低参数的反演误差（Wu 等，2015）。因此，要进一步降低FMF 参数反演不确定性，采用多角度的偏振测量数据进行反演是一个有效途径和方法。

5 结 论

本研究基于OE 反演框架，联合强度和线偏振度的多光谱测量，开展了气溶胶FMF 参数的OE 反演算法研究。基于矢量化的UNL-VRTM 模式进行了前向模拟仿真，分析了线偏振度对FMF 参数的波段敏感性，利用仿真数据对最优化估计反演算法进行了反演一致性测试。研究结果表明：

（1）偏振测量的波段间差异显著，粗粒子对偏振的贡献在长波近红外要高于可见光波段。AOD 不变的情况下，线偏振度作为FMF 的函数存在极小值，且极小值对应的FMF 值随着波长的增加而增加。在670 nm、870 nm 和1610 nm 波段，线偏振度极小值对应的FMF 分别在0.1、0.2 和0.7附近。表明在长波近红外波段，线偏振度对FMF的敏感性高于可见光波段。

（2）OE 反演框架本身具有很好的闭合性。在考虑强度和偏振的测量误差的情况下，基于模拟数据的平均反演残差仅为5.2%，由算法本身引起FMF 反演误差仅为0.014%。给定气溶胶模态参数情况下，基于单角度多光谱强度和线偏振度，OE反演算法能够较好的约束气溶胶的粗、细模态AOD反演，从而得到FMF参数。

（3）线偏振度测量为粗模态气溶胶反演提供了更多的信息量支持，在地基天顶观测模式下，强度联合线偏振度反演，FMF 反演误差从1.4%下降到了0.18%。可见，引入线偏振度测量参与反演能够极大的降低FMF 参数的反演误差，改善反演结果的不确定性。

志 谢前向模拟基于国产卫星载荷同步大气校正仪，传感器波段及带宽等参数设置得到了中国科学院安徽光学精密机械研究所孙晓兵、李梦凡等人的支持，在此表示衷心的感谢！