算法流程如图1所示, 首先针对暗目标法不适用于GF-4数据的情况, 对GF-4/PMI传感器的数据进行了分析, 并根据已有数据和资料对暗目标法进行调整。 利用长时序的AERONET地基数据和PSO聚类算法对京津冀地区的气溶胶类型进行了聚类分析, 并结合实验验证选取了符合实际情况的气溶胶类型。 将聚类得到的气溶胶类型代入6S辐射传输模型建立查找表, 结合调整之后的暗目标法建立类三元一次方程组。 最后选取京津冀区域GF-4卫星影像进行AOD反演, 并与MOD04_L2产品和AERONET地基观测数据进行对比来验证算法的可行性。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 算法流程图Fig.1 Algorithm flow diagram

GF-4卫星是我国第一颗静止轨道遥感卫星, 卫星采用地球同步轨道, 其位与赤道106° E上空36 000 km处, 作为我国高分辨率对地观测重大专项的卫星之一, GF-4卫星上的PMI传感器兼具中空间分辨率和高时间分辨率的成像特点, 相对于其他同步轨道卫星, 不仅具有更加细腻的成像效果, 还能对某一大面积区域进行长期而连续的观测。 GF-4卫星上搭载的PMI传感器参数如表1所示。

表1

Table 1

表1(Table 1)

表1 GF-4卫星PMI传感器参数 Table 1 Parameters of PMI sensor for GF-4 satellite 谱段 谱段范围 /μ m 空间分 辨率/m 幅宽 /km 重访 时间/s 可见光 近红外 (VNIR) 1 0.45~0.90 400 20 2 0.45~0.52 3 0.52~0.60 50 4 0.63~0.69 5 0.76~0.90 中波红外 (MWIR) 6 3.50~4.10 400

表1 GF-4卫星PMI传感器参数 Table 1 Parameters of PMI sensor for GF-4 satellite

假设地球为均匀郎伯表面, 则星载传感器接收到的大气顶层反射率(即表观反射率)ρ _TOA可以表示为^[8]

ρTOA(θs, θv, φ)=ρ0(θs, θv, φ)+T(θs)T(θv)ρs(θs, θv, φ)1-ρs(θs, θv, φ)S(1)

式(1)中, θ _s为太阳天顶角, θ _v为卫星天顶角, φ 为相对方位角; ρ ₀是路径辐射项等效反射率; ρ _s为地表二项反射率; S为大气下界的半球反照率; T(μ _s)和T(μ _v)分别是向下和向上整层的大气透过率, 两者的乘积代表大气总的透射率。

表观反射率受到大气环境下的等效反射率和地表反射率影响。 因此, 地表反射率ρ _s可以根据表观发射率ρ _TOA计算得到, 简化如式(2)

ρs=ρTOA-ρ0(ρTOA-ρ0)S+T(2)

暗目标法也叫浓密植被法, 传统暗目标法利用短波红外波段与红蓝波段之间地表反射率之间的相关关系。 由表1可知GF-4卫星PMI传感器缺少暗目标法所需的短波红外波段, 因此需要对暗目标法进行调整和改进。 假设GF-4卫星PMI多光谱数据中暗像元的红蓝波段地表反射率间具有良好的相关性, 即

ρreds=kρblues+b(3)

式(3)中, k为相关系数; b为相应偏差。

利用2018年9月内蒙古试验场地测得的植被光谱数据, 对其中212条不同的植被光谱数据分析, 得到PMI传感器红蓝波段等效地表反射率的散点分布和拟合结果如图2所示, 拟合之后两者的关系为f(x)=1.643x-0.007, 相关系数达0.926。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 红蓝波段等效地表反射率的拟合曲线Fig.2 Fit curve of red-blue band equivalent surface reflectivity

将红、 蓝波段的表观反射率 ρredTOA, ρblueTOA代入式(2), 再联立式(3)可得类三元一次方程组, 分析可知方程组中未知量为AOD和红、 蓝波段的地表反射率 ρreds和 ρblues。 基于6S模型, 针对研究区域设置相应的大气模式、 气溶胶类型和观测几何条件等参数, 建立不同观测条件和气溶胶类型下AOD和ρ ₀, T, S之间的查找表, 通过插值得到AOD后再代入模型获得地表反射率, 对类三元一次方程组进行求解即可完成AOD的反演。 但是在查表过程中, 每个像元的红蓝波段的地表反射率并不总是相等的关系, 因此可将等式化为求偏差绝对值的问题, 寻找最小的那一组AOD数据即可, 则方程组转化为式(4)

ρreds=ρredTOA-ρ0(ρredTOA-ρ0)S+Tρblues=ρblueTOA-ρ0(ρblueTOA-ρ0)S+Tδ=ρreds-(1.643ρblues-0.007)(4)

粒子群优化算法是由Kennedy和Eberhar对鸟群觅食行为的研究中分析总结出来的寻优算法。 算法的原理如下: 在预设搜索空间中初始化一群具有潜在最优解的粒子, 这些粒子以一定的速度飞行, 粒子搜索的方向和距离由这些粒子的速度决定, 粒子的好坏由适应度值筛选, 粒子在搜索空间中飞行时, 会结合自身和其他粒子的适应度值不断更新自身的适应度值, 完成目标寻优^[9]。

聚类可以理解为样本的排列与组合, 这种组合使得同一类中的样本特征相似度最大, 而不同类之间的样本差异也最大。 PSO聚类算法的过程中, 每个粒子都看作是对样本间聚类的划分, 都可以看成需求最佳聚类中心的可行解, 因此可将粒子结构描述为预定分类数K个聚类中心C_k=(c_k1, c_k2, …, c_kN)(k=1, 2, …, K)的连接, 并编码构造如式(5)

Xi=(c11, c12, …, c1N, c21, c22, …, c2N, …, ck1, ck2, …, ckN)(5)

粒子i可看成K× N维向量, 当所有粒子的位置确定时, 其代表的聚类中心随之确定。 假设共有J个样本, 对粒子进行拆分, 第j个样本的所属类的划分采用最近邻域法

‖Sj-Ck'‖=mink=1, 2, …, K‖Sj-Ck‖(6)

Ck=1mumk∑s∈Sks(k=1, 2, …, K)(7)

则样本S_j属于k'类, S_k指的是第k个聚类中的所有样本, num_k是该类中的样本数。 对该粒子处的各样本进行分类之后, 计算适应度函数值fitness(t), fitness(t)采用类内离散度之和Jc来表示

fitness=Jc=∑k=1K∑h=1numk‖Skn-Ck‖2(8)

式(8)中, S_kn指第k类中的第n个样本, C_k指第k类的聚类中心, t是迭代次数。 当粒子的适应度值小于其个体和全局历史最优适应度值时, 则更新粒子的适应度值。

6S模型模拟了辐射在太阳-地表目标-接收器这一传输过程。 假设研究区域为单一地表, 使用6S模型中的城市型和海洋型、 以及OPAC数据集^[10]中定义的大陆平均型和大陆污染型气溶胶类型, 研究不同气溶胶类型对表观反射率的影响, 如图3所示。 可知当AOD较小时, 四种气溶胶类型对表观反射率影响的差异不大, 当AOD较大时, 差异逐渐增大。 大陆平均型、 大陆污染型和海洋型三种气溶胶类型得到的表观反射率均随着AOD值的增大而逐渐增大, 且增加的幅度不相同, 但是城市型气溶胶类型得到的表观反射率随AOD值的增大而减小。 实验结果说明气溶胶类型对AOD的敏感程度有很大差异, 这种差异来自于气溶胶类型中各种类型粒子不同的吸收、 散射特性。 进一步证明, 气溶胶类型会对AOD反演造成很大误差。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 四种气溶胶类型对表观反射率的影响Fig.3 Effect of four aerosol types on apparent reflectivity

根据长时序持续观测的AERONET地基数据发现, 使用气溶胶微物理特性和光化学特性参数共同表征气溶胶类型, 并利用聚类算法得出气溶胶类型是可行的^{[11, 12]}。 以京津冀地区为研究区域, 选择AERONET地基网络中2015年1月1日到2017年1月1日的北京站点(39.977° N, 116.381° E)和香河站点(39.754° N, 116.962° E)气溶胶的特性数据进行统计分析。

最佳聚类数量是聚类分析的关键参数之一。 聚类数量m不断增加时, 使用PSO聚类算法充分迭代求得各样本到自身所对应聚类中心的欧氏距离之和的最小值, 用能量E表示(即最优适应度fitness), E与m的关系为E=f(m), 根据m对于E的敏感性和拟合关系进行分析, 可得到最佳聚类数量m。 图4反映了两者在实数平面坐标上的散点分布情况, 将其放到对数坐标轴上(图4右上方), 可知能量E和m直线拟合度较好, 证明了在实数坐标空间内E随m变化而指数衰减。 由图4可以看出, 聚类数量m足够大时, 新增聚类所占比例越来越小, 能量E的衰减程度也越来越弱, 表明新增聚类对于气溶胶类型差异性的影响越来越小。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 能量E与聚类数目m的关系Fig.4 Relationship between energy E and cluster number m

引用半衰期序列来进一步观测^[13], 以E₀=f(1)为基准, 即所有的样本只划分为一种类型时对应的能量, 当能量依次减半时取各聚类数量对应的E最为接近的情况, 进而可以得到处于能量半衰期位置的序列p, p对应聚类数目的下限。 由于聚类分析之前要把所有的气溶胶样本归一化, 使得E大幅减少, 因此选择p=2, 即第二个半衰期对应的m=8作为最佳聚类数量。 PSO聚类算法的参数设置如表2所示, 根据PSO聚类算法得到8种气溶胶类型聚类结果。 试验发现, 这8种气溶胶类型中的5种占比高达89.4%, 满足研究区域大多数的大气情况, 故只对聚类结果中的这5种进行讨论。

表2

Table 2

表2(Table 2)

表2 PSO聚类算法的参数 Table 2 Parameters of PSO clustering algorithm Number of categories/m Maximum/ Minimum Inertiaweight (ω max/ω min) Learning factors (c1/c2) Number of trainings (T) 8 0.9/0.4 1.2/1.2 300

表2 PSO聚类算法的参数 Table 2 Parameters of PSO clustering algorithm

该地区气溶胶粒子的体积尺度谱适合用双模态对数分布描述^{[14, 15]}, 所以描述气溶胶物理特性的粒子谱分布参数主要有粒子半径、 体积浓度等, 表征气溶胶光学特性的参数有复折射指数、 单次散射反照率以及不对称因子。 聚类所得5种气溶胶类型体积尺度谱分布特征如图5所示, 复折射指数随波长变化情况如图6所示。 C1, C2和C4这三种类型占据了该地区这两年气溶胶类型的70%。 C1和C2的体积谱相似, 二者中两种模态粒子的体积浓度都比较小, 而且二者细模态粒子峰值半径均在0.15 μ m左右, 5种类型中最小。 C2聚类结果占比最大, 其细粒子和粗粒子的体积占比最低, 复折射指数的虚部最大, 表示很强的吸收能力, 是清洁天气的气溶胶类型; C1的占比较C2小, 其粗细粒子体积浓度比背景值均略有增加, 而且粗粒子浓度增幅大于细粒子, 为晴朗天气的气溶胶类型。 C4和C5体积谱均近似于单峰分布, 但表现相反, C4中细模态粒子占优, 所以细模态粒子的消光贡献远大于粗模态粒子, 且复折射指数实部和虚部也明显降低, 即C4代表着雾霾、 阴霾天气; 而C5体积尺度谱以粗模态为主导峰, 这主要来源于冬季和春季的沙尘型气溶胶, 这种气溶胶的实部较高, 具有较强的散射能力。 对于C3来说, 其体积尺度谱呈现出驼型分布, 即细粒子体积和粗粒子体积相当, 浓度峰值均在0.15 μ m³· μ m^-2以上, 代表严重的灰霾天气。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 5种气溶胶类型体积尺度谱分布Fig.5 Volumetric spectral distribution of five aerosol types

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 5种气溶胶模式复折射指数 (a): 实部; (b): 虚部Fig.6 (a) The birefringence index real part (a) and birefringence index imaginary part (b) of five types of aerosol

基于对5种气溶胶类型特点的分析, C1, C2和C4型更符合京津冀地区夏、 秋季节的气溶胶类型特征。 经过实验验证, C1和C2型气溶胶类型的反演结果相似, 且二者的物理和光学特性非常相似, 为此不分析C2, 以C1, C4以及大陆型气溶胶类型进行AOD的反演。