对中红外波段(2.5~5.0 μ m)的大气组分的透过率进行了仿真。 中红外大气总透过率受水汽和CO₂影响最大[图1(a)], 其他组分仅在特定的吸收位置有一定程度影响(O₃, 4.75 μ m; CH₄, 3.3 μ m; CO, 4.6 μ m; N₂O, 4.5 μ m; N₂, 4.3 μ m), 气溶胶的光谱透过率变化平稳, 保持在较高水平[图1(b)]。 水汽和CO₂的强吸收是导致2.7和4.3 μ m附近大气低透过率(甚至为0)的原因。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 中红外波段大气组分的透过率和光学厚度Fig.1 Transmittance and optical depth of atmospheric components in mid-infrared band

光学厚度(optical depth, OD)表征气体分子的吸收强度, 为了探究两个强吸收谱段的具体范围, 仿真计算了MODTRAN中6种标准大气模式的大气组分OD。 图1(c)和(d)展示了中纬度夏季大气模式下的中红外各气体组分的OD, 落在灰色区域范围内的散点代表OD> 1。 将OD> 1作为判断2.7和4.3 μ m附近处强吸收波段的准则, 统计6种标准大气模式下水汽和CO₂的计算结果(表1)。 对不同标准大气模式下统计的强吸收波段区间取交集, 最终确定中红外大气强吸收波段为2.52~2.83和4.18~4.47 μ m。

表1

Table 1

表1(Table 1)

表1 中红外强吸收波段的统计结果 Table 1 Statistical results of mid-infrared strong absorption bands 标准大气模式 水汽OD> 1的波段区间/μ m CO2 OD> 1的波段区间/μ m 强吸收波段区间/μ m 热带大气 2.50~2.99 2.66~2.83、 4.18~4.48 2.50~2.99、 4.18~4.48 中纬度夏季大气 2.50~2.96 2.66~2.83、 4.18~4.48 2.50~2.96、 4.18~4.48 中纬度冬季大气 2.52~2.87 2.67~2.83、 4.18~4.48 2.52~2.87、 4.18~4.48 亚北极夏季大气 2.51~2.93 2.66~2.83、 4.18~4.48 2.51~2.93、 4.18~4.48 亚北极冬季大气 2.52~2.83 2.67~2.83、 4.18~4.47 2.52~2.83、 4.18~4.47 1976美国标准大气 2.51~2.89 2.67~2.83、 4.18~4.48 2.51~2.89、 4.18~4.47

表1 中红外强吸收波段的统计结果 Table 1 Statistical results of mid-infrared strong absorption bands

1.2.1 地表产品数据

(1) 地表反/发射率数据

采用2022年11月份的MODIS地表反/发射率产品MYD09A1(波段7)和MOD11C3(波段20、 22和23), 重建出两个强吸收波段的地表反射率。 其中, MYD09A1产品提供了大气校正之后的波段1至7的地表反射率估计值; MOD11C3产品通过合成和平均MOD11C1每日文件相应月份的值, 提供了包含中红外波段的月度地表发射率估计值。 图2展示了基于非负矩阵分解方法重建的2.52~2.83和4.18~4.47 μ m两个强吸收波段的地表平均反射率。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 中红外强吸收波段全球平均地表反射率重建结果Fig.2 Reconstruction result of global mean surface reflectance in mid-infrared strong absorption band

(2) 地表温度数据

考虑到地表温度产品在全球尺度空间分布和时间上的连续性, 现有卫星产品很难同时满足时空要求。 采用ERA5逐小时大气廓线再分析数据的最底层地表温度作为温度数据源(https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels?tab=form)。 ERA5数据的空间分辨率为0.25° × 0.25° , 下载2022年11月1日至30日全球区域范围内的大气廓线数据, 根据气压大小将最底层温度提取出来, 然后对30 d的全球地表温度数据取均值, 生产出11月份的温度产品[图3(a)]。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 本研究使用的数据融合产品Fig.3 The data fusion products used in this study

1.2.2 大气产品数据

在中红外大气强吸收波段, 水汽和CO₂是影响卫星入瞳辐亮度的两个关键因素。 除此之外, 云对辐射传输的影响也至关重要。 采用云光学厚度(COT)作为云的物理参数指标, 输入到MODTRAN中进行辐射传输解算, 用以表征云对入瞳辐亮度的影响。 通过对2022年11月1日至30日之间的5km分辨率的MOD05_L2和MOD06_L2产品进行拼接、 融合处理, 获得水汽和COT月均数据产品[图3(b)和(c)]; 值得一提的是, 在处理水汽和COT数据时, 不可避免地会出现经纬度缺失、 数据横跨东西半球等数据处理问题, 这需要有灵活的代码调试经验来解决这些问题。 CO₂数据产品[图3(d)]采用OCO-2卫星数据衍生的月均模式产品, 获取时间为2022年11月份, 空间分辨率为0.5° × 0.625° 。 至此, 本文所使用的月均产品均已生产完毕, 详细的生产流程见图4。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 月均产品生产流程图Fig.4 Flowchart of monthly average product production

1.3.1 基于非负矩阵分解的地表光谱重建方法

非负矩阵分解(non-negative matrix factorization, NMF)算法是一种新的矩阵分解方法^[11], 它通过对基向量与对应的系数增加非负约束, 以确保分解结果均为正值。 通过NMF分解得到的端元向量所组成的非负矩阵与计算得到的权重系数矩阵相乘, 便能重构原始光谱矩阵。 基于NMF进行地物光谱重建^[12], 相比基于主成分分析(PCA)^{[13, 14, 15]}的光谱重建的方法更具有优势, 所提取和采用的端元向量具有明确的物理意义。 对于一个a× b的非负矩阵N_{a× b}(a表示光谱库中数据的个数, b表示每条光谱的波段个数), 存在非负矩阵V_{a× k}和C_{k× b}, 使得

Na×b≈Va×kCk×b(1)

式(1)中, V_{a× k}为基矩阵(端元向量为列向量所组成的矩阵), C_{k× b}为系数矩阵(也称丰度矩阵), k为分解得到矩阵的秩, 且满足k=min(a, b); N_{a× b}中的列向量可解释为对基矩阵V_{a× k}中所有行向量(即端元向量)的加权之和, 且权重系数为系数矩阵C_{k× b}中对应列向量中的元素(丰度)。

当采用欧式距离作为光谱重建的目标函数时, 非负矩阵分解问题被表示成如式(2)最优化问题:

minf(V, C)=N-VC‖F2s.t. V≥0, C≥0(2)

式(2)中, ‖ · ‖ _F表示矩阵Frobenius范数; s.t.为subject to的缩写, 代表约束条件。

1.3.2 中红外辐射传输仿真计算方法

中红外谱段介于可见光-近红外和热红外谱段之间, 兼具短波辐射的反射特性和长波辐射的发射特性, 辐射传输仿真计算具有很强的复杂性^[16]。 太阳和地球的普朗克发射光谱在中红外存在“ 重叠” 区域, 这会导致白天卫星观测的入瞳辐亮度不仅包含地表反射太阳的能量, 还包括地表自身的发射能量^[17]。 图5展示了星载中红外卫星传感器在白天无云条件下的辐射能量传输过程。 传感器在大气层顶接收到的总入瞳辐亮度I_i可以表示为^[18]

Ii(θv, φ)={εi(θv, φ)Bi(Ts)+[1-εi(θv, φ)]Ratm_i↓(θv, φ)+[1-εi(θv, φ)]Ratm_i↓sca(θv, φ)+ρb_i(θv, φ, θs, φs)Eicos(θs)τi(θs, φs)}τi(θv, φ)+Ratm_i↑(θv, φ)+Ratm_i↑sca(θv, φ)(3)

式(3)中, θ _v和φ 分别表示卫星观测路径的天顶角和方位角; θ _s和φ _s分别表示太阳的天顶角和方位角; τ _i和ε _i为通道i的有效透过率和发射率; ρ _{b_i}和T_s为通道i的表面双向反射率和地表温度; R_{atm_i↑} (θ _v, φ )为通道i中大气向上的路径热辐射(图5中的路径②); Ratm_i↑sca(θ _v, φ )为通道i中向上的太阳漫射辐射(图5中的路径③); ε _i(θ _v, φ )B_i(T_s)表示地表直接发射的辐射(图5中的路径④); [1-ε _i(θ _v, φ )]R_{atm_i↓} (θ _v, φ )为地表反射的大气向下的路径热辐射(图5中的路径⑤); [1-ε _i(θ _v, φ )] Ratm_i↓sca(θ _v, φ )为地表反射的向下的太阳漫射辐射(图5中的路径⑥); ρ _{b_i}(θ _v, φ , θ _s, φ _s)E_icos(θ _s)τ _i(θ _s, φ _s)为地表反射的太阳直射辐射(图5中的路径⑦)。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 中红外谱段大气辐射传输过程示意图Fig.5 Schematic diagram of atmospheric radiative transfer process in the mid-infrared spectrum

因此, 中红外辐射传输仿真计算不仅要考虑地表自身的辐射, 太阳辐射同样也不可以忽略。 最终中红外入瞳辐亮度应该由总的地表反射辐射、 地表发射辐射、 大气路径热辐射和向上的太阳漫射辐射四部分辐射分量构成。 辐射分量的空域变化, 会直接引起背景杂波的空域变化, 并影响变化区域的杂波起伏水平。

1.3.3 背景杂波起伏水平计算方法

利用11× 11像元的方窗邻域统计法计算单块小区域杂波起伏DN_area的数值, 表达式为

DNarea=∑i=111×11(DNi-DN¯)211×11(4)

DN¯=∑i=111×11DNi11×11(5)

式(4)和式(5)中, DN_i为单像元i的背景辐亮度, DN¯为11× 11方窗范围内背景辐亮度均值水平。

将11× 11的方窗逐行列在背景杂波图像上进行平移(从左至右, 从上至下), 迭代计算每个方窗内的杂波起伏水平。 之后, 统计DN_area的概率直方图, 并换算到95%概率下的杂波起伏。

2.1.1 有云情形

(1) 2.52~2.83 μ m波段杂波强度仿真计算结果

2.52~2.83 μ m波段的背景杂波强度主要受到云和水汽的耦合影响, CO₂的影响非常微弱, 来自地表的能量几乎完全被抑制。 云对杂波强度的影响主要是“ 散射增强效应” , 即杂波强度与COT呈现非线性的正比关系。 非洲中部沙漠区、 阿拉伯半岛沙漠区(尤其南部沙漠区)以及印度北部污染区这些几乎无云区域的杂波强度非常低, 在图6(a)中显示为“ 深蓝色” 。 而在欧洲东部、 北美洲部分COT的高值区, 杂波强度较大。 水汽的强吸收使得杂波强度分布整体上与水汽含量呈现非线性反比关系, 热跌区域水汽高值区的杂波强度较低, 中纬度区域水汽低值区展现较高的杂波强度。 值得一提的是, 青藏高原和南美洲安第斯山脉的杂波强度异常高, 这是高海拔、 高大气透过率和极低水汽含量等因素导致的。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 全球背景杂波强度分布Fig.6 Global background clutter intensity distribution

(2) 4.18~4.47 μ m波段杂波强度仿真计算结果

随着波长向长波方向移动, 云对杂波的影响由“ 散射增强” 转变成了“ 消光减弱” 。 4.18~4.47 μ m波段的背景杂波强度主要受到云消光作用影响。 在非洲中部、 阿拉伯半岛南部、 印度中部和安第斯山脉中部等COT的低值区, 杂波强度非常高, 其他的COT高值区呈现较低的杂波强度[图6(b)]。

2.1.2 无云情形

(1) 2.52~2.83 μ m波段杂波强度仿真计算结果

无云情形下的辐射传输计算结果能够清晰地表征背景数据自身对杂波强度的影响。 无云时, 2.52~2.83 μ m波段的背景杂波强度与水汽含量呈现非线性的反比关系。 从全球尺度来看, 青藏高原、 安第斯山脉以及高纬度区域, 因降水较少, 会出现明显的杂波强度高值区[图6(c)]。

(2) 4.18~4.47 μ m波段杂波强度仿真计算结果

经过详尽的敏感性测试, 水汽在4.18~4.47 μ m波段的敏感性非常弱, 几乎可以忽略不计; 在4.20~4.38 μ m波段, CO₂含量与背景杂波强度正相关, 而在4.18~4.20和4.38~4.47 μ m波段处, CO₂含量与背景杂波强度负相关, 整体上抵消了CO₂的影响。 因此无云时, 4.18~4.47 μ m波段的背景杂波强度主要与地表要素相关, 与地表温度呈现非线性的正比关系[图6(d)]。

2.2.1 有云情形

(1) 2.52~2.83 μ m波段杂波起伏水平计算结果

2.52~2.83 μ m波段全球背景杂波起伏水平具有“ 点状峰值、 区域增强、 整体低值” 的明显特征[图7(a)]。 青藏高原附近的杂波起伏最大, 并且呈现从中心向周围递减的趋势; 非洲北部、 阿拉伯半岛、 印度半岛、 澳大利亚西南部、 安第斯山脉中部等地区呈现“ 区域增强” 的特征。 而全球其他区域呈现“ 整体低值” 的特征, 大小基本保持在10^-4量级附近。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 有云情形下全球背景杂波起伏水平分布Fig.7 Global background clutter fluctuation level distribution under cloud condition

(2) 4.18~4.47 μ m波段杂波起伏水平计算结果

4.18~4.47 μ m波段全球背景杂波起伏水平呈现“ 区域高值、 片状增强、 高纬低值” 的特征[图7(b)]。 非洲北部、 阿拉伯半岛、 印度半岛、 青藏高原、 澳大利亚西南部和安迪斯山脉中部的杂波起伏有“ 区域高值” 的特征。 太平洋中部、 大西洋中部、 科迪勒拉山系及沿岸海域、 古巴岛附近及沿岸海域等呈现“ 片状增强” 的趋势。 而在纬度较高的上下两侧区域呈现“ 高纬低值” 特征。

(3) 概率直方图统计结果

图8是有云情形时, 换算95%概率下的全球背景杂波起伏直方图。 在0至2× 10^-3 W/m²/sr/μ m的低水平杂波起伏范围内, 2.52~2.83 μ m波段的杂波起伏比例(约85%)大于4.18~4.47 μ m波段的杂波起伏比例(约60%); 在2× 10^-3~6.5× 10^-3 W/m²/sr/μ m的较高水平杂波起伏范围内, 2.52~2.83 μ m波段的杂波起伏比例小于4.18~4.47 μ m波段的杂波起伏比例。 因此, 从全球尺度杂波抑制的角度, 在有云情形下, 2.52~2.83 μ m波段的红外目标探测性能整体优于4.18~4.47 μ m波段。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 有云情形下全球背景杂波起伏水平直方图Fig.8 Histogram of global background clutter fluctuation level under cloud condition

2.2.2 无云情形

(1) 2.52~2.83 μ m波段杂波起伏水平计算结果

无云时, 2.52~2.83 μ m波段全球背景杂波起伏水平在水汽低值区有较大的增幅。 青藏高原地区的杂波起伏水平出现峰值, 并从中心向周围减弱; 安第斯山脉中部、 高纬地区等水汽低值区则出现了不同程度的起伏水平增强[图9(a)]。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 4.18~5.47 μ m波段全球背景杂波起伏水平分布Fig.9 Global background clutter fluctuation level distribution in 4.18~4.47 μ m band

(2) 4.18~4.47μ m波段杂波起伏水平计算结果

无云时, 4.18~4.47 μ m波段的整体杂波起伏水平较低, 仅在安第斯山脉及沿岸海域、 非洲的西南部及沿岸海域、 澳洲西南部及沿岸海域出现高值区, 其他区域的杂波起伏基本控制在2× 10^-4 W/m²/sr/μ m以内[图9(b)]。

(3) 概率直方图统计结果

图10是无云情形时, 换算95%概率下的全球背景杂波起伏直方图。 2.52~2.83 μ m波段的杂波起伏有近80%落入0~10^-3范围内, 而4.18~4.47 μ m波段的杂波起伏大小量级均在10^-4。 因此, 在无云情形下, 4.18~4.47 μ m波段的红外目标探测性能整体优于2.52~2.83 μ m波段。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 无云情形下全球背景杂波起伏水平直方图Fig.10 Histogram of global background clutter fluctuation level under cloud-free condition