卫星遥感气溶胶光学厚度（aerosol optical depth，AOD）受卫星轨道间隔、天气状况和云遮挡等因素影响，空间覆盖度偏低。为提高卫星遥感AOD数据空间覆盖度，有学者从多源AOD数据插值和时空地统计关联角度开展了系列研究，多源AOD数据插值方面的研究通过构建多源AOD数据间函数关系，基于一种AOD值估算另一种产品对应位置的AOD缺失值^[1-6]。该类方法虽然能够通过多源产品融合提升AOD空间覆盖度，但存在如下问题：①提升空间覆盖度的同时降低了数据本身的空间分辨率；②空间覆盖度优化过程中仅考虑AOD产品间的时间关联，忽略了AOD产品间的空间关联关系，时空联合关联特征考虑不全面。

根据地理学第一定律可知，任何事物间都存在时间和空间上的相互关联^[7]。鉴于此，AOD时空地统计关联方法理论上可通过考虑多源AOD产品间的时间或空间相关性，实现缺失位置AOD值填补。目前，充分考虑多源AOD数据间空间关联特征的泛克里金、空间统计融合等方法已经被广泛应用于提升卫星遥感AOD产品的空间覆盖度^[8-9]。然而，现有应用仅考虑AOD数据空间关联，未顾及时间关联特征，这将影响空间上严重缺失而时间序列上较完整区域AOD空间覆盖度提升的效果^[10-11]。也有学者试图采用贝叶斯最大熵估计方法从时空关联角度开展AOD融合实验^[12]，但该方法算法复杂且计算效率低，难以满足大范围、长时间序列应用需求。

针对上述AOD空间覆盖优化与提升方法中存在的AOD数据时空关联关系考虑不全面以及计算复杂、应用困难等问题，空间覆盖时空随机效应（spatio-temporal random effects，STRE）模型为解决上述问题提供了一种新的方案。STRE是一种将高维时空数据降到固定维时空数据的时空Kalman模型^[10-14]。然而，STRE模型在现有遥感时空数据处理中存在以下缺陷：①时空趋势去除前，通常对原始数据进行对数变换，使数据满足高斯分布，容易导致估算结果出现奇异值；②STRE建模参数估计采用Method-of-Moments方法，虽具有计算效率高的优点，但要求数据在时空上均匀分布，日均遥感AOD数据难以满足；③原始数据长时间和大范围连续缺失条件下，STRE模型估算结果存在较大不确定性。针对此，本文提出一种不确定性约束下的AOD空间覆盖STRE建模优化方法。

1 研究方法与数据 1.1 STRE模型

STRE模型是一种分层时空模型，模型将时空信息分解为全局趋势信息、局部信息、细微信息和观测噪声4层，采用多尺度空间基描述信息的空间依赖和一阶自回归模型描述信息的时间依赖关系。观测值Z(s, t)由真值Y(s, t)和观测噪声ε_t(s)两部分组成，即^[10]：

$ Z\left( {s, t} \right) = Y\left( {s, t} \right) + {\varepsilon _t}\left( s \right) $

(1)

式中，s为观测点空间位置；t为观测时刻。Y(s, t)可分解为全局时空趋势项μ_t(s)和局部信息υ(s, t)以及细微信息ξ(s, t)，即：

$ Y\left( {s, t} \right) = {\mu _t}\left( s \right) + υ\left( {s, t} \right) + \xi \left( {s, t} \right) $

(2)

对于AOD全局时空趋势，本文采用时空移动窗口去除，使其服从高斯分布、且均值为0。选取的时空窗口大小为（2M+1）×（2N+1）×（2T+1），求取窗口内AOD平均值作为中心像素所在位置趋势值。

对于局部信息υ(s, t)，一般采用空间基函数表示，即：

$ \upsilon (s, t) = {S_t}(s)a(t) $

(3)

式中，a(t)=[a₁(t), a₂(t)…a_r(t)]'为系统状态量; 空间基函数S_t(s)=[S₁, _t(s), S₂, _t(s)…S_r, _t(s)]'为随时间变化的r阶函数; υ(s, t)在时间上是相对平稳的，即在任意时刻υ(s, t)其均值为0^[15]。本文选择Bisquare空间基^[15]，具体形式如下：

$ S(s) = \left\{ \begin{array}{l} {\left\{ {1 - {{(\left\| {s - c} \right\|/g)}^2}} \right\}^2},\left\| {s - c} \right\| ＜ g\\ 0,其他 \end{array} \right. $

(4)

式中，c为空间基中心位置；g为空间基变程。为更好地模拟AOD空间变化，本文采用四叉树方法设置不同空间尺度上的空间基，直到空间基拟合残差V稳定：

$ V = {Z_t} - {\mu _t} - {S_t}({({S_t}'{S_t})^{ - 1}}{S_t}'{\rm{ }}({Z_t} - {\mu _t}){\rm{ }}) $

(5)

采用半变异函数分离式(5)得到的残差包含细微信息ξ(s, t)和观测噪声ε_t(s)，并选取AOD空间覆盖度较高的几期卫星遥感AOD数据标定上述球状半变异函数参数，作为STRE模型参数估计初始值。在此基础上，采用最大期望（expectation maximization，EM）迭代方法估计STRE模型参数$\theta = {\rm{ }}\{ \sigma _{\xi (t)}^2, \phi , Q, a_{0{\rm{ }}}^0, {K_0}\} $，其中$a_{0{\rm{ }}}^0 $与K₀分别为初始状态量a(0)的均值和协方差矩阵。

通过EM迭代参数估计和卡尔曼平滑等过程得到STRE模型参数θ及a(t)、ξ(s, t)，在此基础上，可基于式(6)估算实验区域内AOD缺失时空点位{s₀, t₀}的AOD值，同时基于式(7)估算该AOD值对应的不确定性大小：

$ \mathop {{\rm{ }}Y}\limits^\frown ({s_0},{t_0}) = {\mu _{{t_0}}}({s_0}) + S({s_0})a_{{t_0}}^{\rm{T}} + \mathop \xi \limits^\frown ({s_0},{t_0}) $

(6)

$ \begin{array}{*{20}{l}} {{\sigma ^2}{{({s_0},{t_0})}^{{\rm{FRS}}}} = S({s_0})P_{{t_0}}^{\rm{T}}S({s_0})' + \sigma _{\xi (t)}^2} \end{array} $

(7)

式中，FRS表示基于STRE模型的Kalman平滑方法。

尽管STRE模型可估算任意AOD缺失时空点位的AOD值，但是估算结果均存在一定的不确定性。特别是在AOD空间上严重缺失而时间序列上较完整的情况下，估算结果常存在较大不确定性，本文采用不确定性大小对AOD估算结果进行约束分析。具体实现过程如下：

1）给定不确定性大小初始阈值$\sigma _0^2 $，本文将所有时空点位σ²(s₀, t₀)^FRS的上四分位数作为初始阈值$\sigma _0^2 $；

2）剔除不确定性大于阈值$\sigma _i^2 $的AOD估算值，联合地基AERONET（aerosol robotic network）AOD观测数据，重新评定不确定性约束后AOD估算数据精度；

3）设定阈值下降率为α，重新设定阈值为$\sigma _i^2 = \sigma _{i - 1}^2 - \alpha $；

4）重复第2）~3）步，直到不确定性约束下AOD估算结果与地基AERONET AOD观测数据间的相关系数R>0.7。

相对已往STRE模型，本文从AOD时空趋势去除、参数估计方法及AOD估算结果不确定性约束3个方面对STRE建模过程进行了改进。STRE具体建模实现流程如图 1所示，实验过程采用MATLAB语言编程实现。

1.2 实验数据

近年来，雾霾已成为困扰我国公众健康的首要空气污染问题。城镇化水平高、人口密度大的胡焕庸线东部地区由此成为了我国大气污染治理与健康危险损害防控的重要区域。本文对此选择胡焕庸线以东区域作为研究区（以下简称东部地区），开展不确定性约束下的AOD空间覆盖STRE建模优化方法检验。原始AOD产品（http://ladsweb.nascom.nasa.gov/）为Terra卫星搭载的中分辨率成像光谱仪（Moderate-resolution Imaging Spectroradiometer，MODIS）传感器Collection 6 10 km MOD04 DT_DB_Combined数据，STRE建模优化后AOD数据空间分辨率为10 km，时间分辨率为日，空间关联窗口大小为49×49格网，时间关联窗口为3 d^[15]。AOD数据时段跨度为2013-01-19—2017-12-31，日均包括影像13~15幅，共计约28 000幅。MODIS AOD产品采用层级数据格式（hierarchical data format, HDF）文件结构存储，基于MODIS Toolkit采用IDL批量处理程序对HDF数据进行投影、拼接以及重采样等操作。为确保数据有效性，采用ArcGIS栅格空间分析工具筛选栅格质量控制标志QA=3的AOD数据作为有效AOD值。

采用地基AERONET Version 3.0 Level 2.0级AOD观测数据（https://aeronet.gsfc.nasa.gov/），对MODIS AOD产品和STRE建模优化AOD数据进行验证和精度评价。选择MODIS卫星上午过境（10:30）前后1 h内地基观测数据平均值作为AOD真实值，最终参与精度评价的AERONET站点包括Beijing、Beijing-CAMS、Xianghe、Xinglong、Taihu和Hong_Kong_Sheung 6个。利用AERONET中440 nm和870 nm波长构建440~870 nm波长指数（详细过程参见文献[16]）。

2 结果分析 2.1 AOD空间覆盖优化精度评定

图 2为原始MODIS AOD、空间覆盖建模优化后AOD（STRE AOD）与AERONET地基观测AOD（AERONET AOD）验证散点图。从图 2(a)中可知，东部地区2013-01—2017-12期间共成功匹配MODIS AOD-AERONET AOD点对568个，点对均方根误差（root mean square error，RMSE）为0.44，偏差绝对平均值Abs_basis为0.28，相关系数R为0.74。与AERONET AOD相比，MODIS AOD产品存在高估现象，其中偏差绝对值大于1.0的AOD点共计20个，主要出现在夏季，这是由于夏季降水较多、温度较高、大气湿度增加、风速较小，容易导致气溶胶的吸湿增长现象^[17]。对比图 2(a)、2(b)可知，STRE AOD-AERONET AOD匹配点对数量由原来的568个提升至942个，这说明STRE建模优化可以有效提升AOD空间覆盖度。同时，点对AOD之差的RMSE下降了0.10，偏差绝对平均值Abs_basis下降了0.04。但在匹配点对中存在部分点对的STRE AOD值小于0，这说明在STRE建模优化部分结果存在较大的不确定性。图 2（c）为不确定性约束下STRE建模优化AOD（uncertainty constraint STRE AOD，C_C_STRE AOD）与AERONET AOD结果对比。对比图 2(b)、2(c)可以发现，不确定性约束分析可将STRE建模优化估算结果中小于0的AOD值剔除，因此其匹配点对数量相较于不确定性分析前有一定的减少；但其R相较于不确定性分析前提升了0.06，RMSE和Abs_basis分别比之前下降了0.03和0.02，这说明不确定性分析可有效提升估算结果的验证精度。

2.2 优化前后AOD空间覆盖度对比

AOD空间覆盖度定义为区域内有效像元数量（指原始卫星遥感观测和不确定性约束STRE建模优化AOD数据值大于0的像元）占区域内总像元数量的百分比。

图 3表示2013-01—2017-12东部地区逐日MODIS AOD和不确定性约束下STRE建模优化AOD空间覆盖度统计结果图。对比图 3中两种AOD数据的空间覆盖度分布可知，研究期间不确定性约束下STRE建模优化AOD数据空间覆盖度相较于原始MODIS AOD产品有明显提升。就年均而言，原始MODIS AOD产品空间覆盖度约为15%，而不确定性约束下STRE建模优化AOD数据年均空间覆盖度保持在50%左右。

就季节而言，由于秋季湿度降低、降水减少、大气中云量减少等因素的影响^[18]，研究区内秋季原始MODIS AOD平均空间覆盖度最高，为19.42%，与之对应的不确定性约束下STRE建模优化AOD空间覆盖度平均提升幅度最大，达36.47%；而冬季受地表雪覆盖且植被落叶的影响，地表反射率提高^[19]。因此，原始MODIS AOD数据冬季平均覆盖度最低，仅为11.43%，不确定性约束下STRE建模优化AOD空间覆盖度平均提升幅度最小，为30.02 %。

对于日均粒度而言，与原始MODIS AOD日均产品空间覆盖度相比，2013—2017年不确定性约束下STRE建模优化AOD空间覆盖度提升平均值达到33.57%，单日覆盖度提升最大值为75.58%，日均空间覆盖度最大值达97.77%。

图 4展示了2013年8月27日MODIS AOD、STRE建模优化AOD（STRE AOD）、STRE建模优化估算结果标准差（STRE标准差）以及不确定性约束下STRE建模优化AOD（C_C_STRE AOD）空间分布情况。从图 4(a)中可知，2013年8月27日东部地区MODIS AOD空间覆盖度非常低, 仅为11.52%，而STRE AOD空间覆盖度可提升至96.01%。

图 4(c)表明, 云南西部、广西西部以及陕西、山西等区域像元AOD估算结果标准差较大、存在较大的不确定性。这一现象主要原因是原始MODIS AOD产品在较大时间和空间变程内严重缺失，导致STRE模型所获取的AOD时空变化信息非常有限。为减少STRE建模优化中不确定性较大估算结果对整体建模优化结果的影响，本文采用不确定性约束对STRE AOD估算值质量进行了控制，经过参数动态调整，最终将像元位置标准差大于0.045的AOD估算结果视为异常估算结果。经过不确定性约束后的STRE AOD空间分布如图 4(d)所示，其空间覆盖度为78.65%。

2.3 优化前后AOD时间完整度对比

AOD时间完整度为指定时段区域内某空间位置像元AOD值有效天数占总天数的百分比。图 5展示了不确定性约束下STRE建模优化前后AOD时间完整度的空间分布情况。其中，研究区内2013—2017年原始MODIS AOD产品像元时间完整度平均值为16.04%；不确定性约束下STRE建模优化后像元时间完整度平均值为50.59%。对比图 5(a)和5(b)可知，研究区域内42 450个像元时间完整度在不确定性约束下STRE建模优化后都有一定提升，最小值为13.47%，最大值为79.14%，平均值为34.56%。

各省、市原始MODIS AOD与不确定性约束下STRE建模后AOD时间完整度统计结果见表 1。

表 1表明，分省市而言，东部地区除北京市（30.20%）和天津市（30.28%）以外，其他省份和直辖市原始MODIS AOD产品时间完整度平均值均在30%以下。其中，四川省的平均时间完整度最低，仅为6.65%，主要原因是该地区全年云量覆盖率较高，常年总云量可达75%以上，MODIS AOD反演缺失率严重。不确定性约束STRE建模优化后AOD时间完整度得到了明显提升。其中，北京市、天津市、河北省以及辽宁省等13个省和直辖市的时间完整度平均值达到50%以上，北京市平均时间完整度最高为80.36%；贵州省平均时间完整度为33.10%，在所有省和直辖市中最小，气候湿润和云雨天气较多等气象条件是该地区MODIS AOD严重缺失的主要原因。

图 6为AERONET Beijing观测站2013-03-01—2014-02-28原始MODIS AOD、不确定性约束下STRE建模优化后AOD以及地基AERONET观测AOD日均数据时间序列对比图。

从图 6可知，该时段共364 d，其中AERONET AOD日均观测数据为69个，MODIS AOD日均数据总数为215个，而C_C_STRE AOD日均数据总数为360个。同时，MODIS AOD与AERONET AOD匹配数量为53个，其绝对偏差Abs_basis为0.29；而C_C_STRE AOD与AERONET AOD匹配数量为69个，其绝对偏差Abs_basis为0.23，从上述指标的对比可以发现，不确定性约束下STRE建模优化不仅可以显著提升AOD时间序列完整度，还能够提升优化AOD的整体验证精度。

3 结语

本文针对现有AOD空间覆盖优化研究中存在融合过程时空联合关联特征考虑不全面（仅考虑AOD产品间的时间关联，而忽略了AOD产品间的空间关联关系）、估算结果不确定性未加以约束等问题，提出了一种不确定性约束下的AOD空间覆盖STRE建模优化方法，并采用全球AERONET地基观测AOD数据对STRE建模优化前后以及不确定性约束下的STRE AOD数据精度进行了评价和结果分析。研究结果表明，通过综合顾及AOD数据时空关联关系，不确定性约束下STRE建模可有效地提升原始AOD产品空间覆盖度和时间完整度。研究时段实验区卫星遥感AOD空间覆盖度由15.96%提升至49.53%，单日覆盖度最大提升75.58%；时间完整度由16.04%提升至50.59%，提升幅度最大为79.14%。同时，不确定性约束下的STRE建模优化在提升原始MODIS AOD产品空间覆盖度和时间完整度的基础上，还可进一步优化提升STRE建模优化结果的精度，降低原始AOD长时间和大范围连续缺失所造成的不确定性影响。