引言

近年来, 由于我国沙尘暴频繁发生, 以及人们生活水平的提高和对环境问题的日益关注, 沙尘暴研究已经渗透到自然科学和社会科学研究的许多领域, 包括农业气象、环境生态、交通运输、生活健康、区域和全球气候等等^[1]。对沙尘暴发生区域、强度、输送和沉降的准确监测是所有研究必须首先解决的问题。通过地面设立监测站来监测沙尘暴, 受到人力、物力、财力和技术的限制, 很难进行大规模布点, 目前比较完整的地面观测主要依赖于气象观测站, 但现有地基观测站缺乏沙尘气溶胶观测的专用仪器。大范围系统性监测沙尘暴的发生、移动和沉降, 以及定量遥感空气中的载沙量只有依靠卫星遥感来解决。

卫星遥感对于大范围监测沙尘区域分布和沙尘强度, 研究沙尘气溶胶与地区气候的相互作用是比较理想的手段之一。利用气象卫星的可见光通道进行海上沙尘信息的反演有许多成功的研究^[2-3]。海洋表面均一性好, 且可见光的反射率低, 沙尘云作为亮目标叠加在暗的海洋背景上, 利用它对太阳光的后向散射可以反演出沙尘层光学厚度, 并转换成沙尘浓度。相似的方法用在陆地上空就会出现问题, 因为下垫面背景季节和地域的变化大, 地表光学特性及其变化经常难以提供。由于裸露地表和沙尘层都是高反射目标, 在干旱和半干旱地区上空的沙霾与背景难以分辨, 尽管有些研究者曾试图进行沙漠亮目标上空的气溶胶卫星反演试验^[4], 但没有形成业务产品。Herman等^[5]曾经利用Nimbus 7卫星上TOMS紫外探测通道计算陆地和海上的气溶胶指数, 包括沙漠上空的沙尘气溶胶, 他们的研究是基于地表在紫外波段反射率较低的特征, 而且陆地和海洋几乎是常数, 所以其研究结果得到了广泛的认同。

欧洲开发了一种陆地上空沙尘遥感监测的新方法, 该方法利用气象卫星的红外探测通道, 如Meteosat卫星的10.5~12.5 μm通道观测数据。这一方法的原理是:空中沙尘加大了大气对太阳短波辐射的后向散射, 同时阻挡了地球向外太空的热辐射, 从而减少了卫星对地观测的红外辐射。欧洲科学家利用这一红外亮温衰减作为反映沙尘强度的间接定量参数, 并对Meteosat长时间序列资料进行处理, 处理结果用于沙漠气候的相关分析^[6-7]。

沙尘目标在红外窗区还有一个显著特征, 就是红外两个分裂窗通道亮温差 (BT₁₁与BT₁₂差) 在沙尘区是负值, 它是识别沙尘目标的关键因子, 很多研究者对这一特征进行了阐述和应用^[8-9], 不过某些地表和云也有类似的负值特征, 因此需要联合其他特征参数进行综合判别。除热红外通道外, Ackerman^[10]的早期研究发现, 中红外3.7 μm通道和热红外窗区通道的亮温差值在沙尘区比其他目标亮温差值要大, 特别是在白天更加明显 (主要由于沙尘粒子对太阳光中红外的后向散射较强), 这也是用来跟踪遥感沙尘目标的重要特征参量。

国内的沙尘暴卫星遥感监测应用研究从20世纪90年代开始, 主要利用极轨气象卫星NOAA/AVHRR的多通道数据进行沙尘目标判识^[11-13], 其基本原理是采用多光谱分类和图像解译的方法判识沙尘目标, 采用的通道主要以可见近红外通道为主, 并通过人机交互的方式生成沙尘暴监测图像, 同时还开展了卫星遥感沙尘光学厚度和载沙量初步研究^[14]。这些研究大部分是基于极轨卫星数据, 监测方法多处于准业务非自动化阶段, 需要人工干预才能生成沙尘暴监测产品, 而且极轨气象卫星在白天只有一次过境观测, 对于沙尘这种演变较快且移动迅速的目标存在明显的局限性。胡秀清等^[15-17]自2001年着手开展用静止气象卫星进行自动沙尘暴监测研究, 本文是在以上工作基础上的补充和完善, 重点引入了反映沙尘定量信息的参数IDDI, 主要介绍利用沙尘目标在红外分裂窗通道的两个特征, 结合多通道光谱聚类方法, 并综合采用其他一系列遥感资料处理技术, 开发了一套静止气象卫星沙尘暴自动监测业务算法, 该算法不仅应用于GMS-5和GOES-9等卫星资料, 而且已经成为我国新发射的静止气象卫星FY-2C的业务产品。

1.   原理

1.1   分裂窗亮温差

根据物理实际测量结果^[8], 如图 1所示颗粒状的石英沙在两个粒径范围 (0~75 μm, 75~250 μm) 的比辐射率随波长变化, 图中显示空气中的沙尘 (粒径一般小于5 μm) 在热红外波段11 μm附近比辐射率相对较小, 也就是吸收性弱, 这一特征也叫透明特征, 随着粒子尺度减小而透明性增强。另外悬浮在空气中的沙尘粒子之间间隔一般数倍于粒子半径 (强沙尘暴除外), 而地表沙土粒子一般紧密相邻。研究表明, 粒子至少被3倍于粒子半径的空间分离才会产生独立散射。对于不透明紧密填充的介质在热红外波段不产生独立散射, 此时Mie散射计算的反照率和相函数就不准确, 这种悬浮颗粒的密集效应直接影响到它的比辐射率光谱对比度, 沙尘粒子越密集比辐射率光谱差异越小, 越疏松比辐射率光谱差异越大。因此悬浮在空气中的沙尘在11 μm和12 μm通道的比辐射率差异与沙尘粒子半径和沙尘强度有关, 但总趋势是11 μm通道比辐射率小于12 μm通道, 因此这两个通道的亮温差通常是负值, 但差值绝对值大小影响因素较多。

图  1  纯硅酸盐石英沙热红外光谱比辐射率^[9]

Figure  1.  Thermal infrared emissivity spectrum of pure quartz 75~250 μm and 0~75 μm of particle size^[9]

下载: 全尺寸图片 幻灯片

卫星观测到的沙尘目标信号可以表示为:

式 (1) 中, R_λ为卫星观测到的目标红外辐射, B_λ为Plank函数, T_s与T_a分别代表地表温度和有效大气温度, τ_λ为沙尘气溶胶消光光学厚度, exp (-τ_λ) 相当于沙尘层的透过率, 1-exp (-τ_λ) 就是沙尘层的吸收率, 也是沙尘的发射率。式 (1) 表明:卫星观测到的红外辐射由两部分组成, 一部分是由地表发射, 经大气消光后到达卫星传感器, 另一部分是大气本身发射的红外辐射。当沙尘气溶胶光学厚度不大时, 其发射作用相对较弱, 卫星观测到的主要是透过沙尘气溶胶层到达卫星的地表发射辐射。由于沙尘气溶胶在两个分裂窗通道上的透过率存在差异, 卫星探测到沙尘目标的这两个通道亮温便产生了差异。通过对这种分裂窗亮温差的分析就可以推算沙尘的存在。这就是分裂窗监测沙尘天气的基本原理。

图 2是利用MODTRAN4对沙尘性气溶胶的模拟结果。图 2说明不同能见度的大气条件可以产生-2~0.4 K的分裂窗亮温差。图中的模拟结果说明沙尘强度变化对分裂窗亮温差的影响, 当沙尘在开始由弱变强 (能见度降低, 但大于0.5 km) 时, 分裂窗亮温差是负值, 绝对值由小变大, 这主要是由于穿过半透明沙尘层的地面信号越来越小, 但在11 μm通道减小得更快。云图图像分析表明, 某些卷云同样会产生很小的分裂窗负亮温差。单纯通过分裂窗图像亮温差还不能完全准确地区分沙尘和某些云目标。

图  2  中纬度冬季标准大气条件下不同强度的沙尘分裂窗亮温差随地面能见度的变化

Figure  2.  Brightness temperature difference (BTD) in split window dependence of dust loading in several winter surface temperatures in mid-latitudes

下载: 全尺寸图片 幻灯片

1.2   红外差值沙尘指数

大气中悬浮的沙尘因为对太阳辐射的后向散射, 会减少入射地面的太阳短波辐射, 使地表接受的太阳辐射减弱升温减缓, 同时地面向外空发射的红外辐射被沙尘层部分吸收, 降低地表向外空的红外辐射, 这两个效应都造成地表向空间的红外辐射减少。悬浮在空中的沙尘温度通常比地面低, 沙尘层的高度越高, 温度就越低, 沙尘层的自身发射不能补偿它对地表向空间红外辐射的衰减。这些过程的综合影响是, 沙尘粒子将导致大气外界的红外辐射较强衰减, 红外差值沙尘指数 (infrared difference dust index, IDDI) 就是表征大气沙尘气溶胶造成地气系统亮温的衰减量, 它可以通过GMS, METEOSAT, FY-2等卫星的红外通道进行计算。IDDI定义为卫星观测到实时目标亮温减去地表背景亮温, 地表背景亮温就是同时刻晴空大气地表亮温。

大气中沙尘导致的亮温衰减量 (相对于晴空红外辐射), 可以作为沙尘的强度指数。模式研究表明:地面和500 hPa之间的沙尘层, 光学厚度在0.5~3将会带来10~30 K的红外亮温衰减^[6]。对于0.25反照率且光学厚度为1的沙尘, 无论粗粒子 (>1 μm) 还是细粒子 ( < 1 μm) 都会导致红外辐射和亮温的减少, 因此IDDI可以解译为沙漠和半干旱地区大气中沙尘强度的表征参量, IDDI值越大代表空气的沙尘越多, IDDI越小代表空气中的沙尘越少。IDDI产品现在已经成为METEOSAT业务产品, 并积累了20年的历史数据集, 在日产品基础上合成月平均产品。IDDI数据集为欧洲研究沙尘的起源和沙尘的气候效应提供基础数据。

IDDI是地气系统由于沙尘气溶胶存在所带来的亮温衰减, 根据这一原理利用亚洲区域的静止气象卫星数据可以生成中国北方地区的IDDI日产品。IDDI的生成过程包括如下步骤:首先建立地表背景亮温图像, 它代表一段连续时间相对晴空条件下的热红外图像, 这一段时间间隔要足够短, 以避免地表的季节影响, 但也要持续时间足够长, 以确保每个像元至少有1次晴空观测数据, 为了满足这两个条件, 选用10 d作为生成背景图像的参考时段。对于某一个像元, 将参考时段内单天的最大亮温值作为目标的背景亮温。第二步, 背景图像减去当天的观测图像得到每日的亮温差值图像, 它表示大气成分 (气溶胶、水汽或云) 造成的亮温衰减。第三步, 滤除云的覆盖区域, 剩下区域的IDDI值就是大气气溶胶的定量参数, 排除通常的气溶胶和水汽的影响, 就是沙尘的IDDI值。这里云检测是一个关键步骤, 下文将描述云检测算法思路。晴空大气水汽带来的地表亮温的衰减量, 在北方干燥区域相对于沙尘的亮温衰减是一个小量 (0~3 K), 对于沙尘的判视来说可以忽略不计。因此剔除云区后得到的就是某个时刻IDDI图像, 多天的IDDI合成平均值就得到月平均IDDI。

1.3   IDDI检验和误差分析

为了使IDDI能定量表征大气中的沙尘含量, 有必要建立IDDI同现有其他表征大气气溶胶参数的函数关系, 如气溶胶光学厚度AOD或地面能见度。Legrand^[6]利用西北非洲39个站的能见度资料对IDDI进行比对, 比对时以观测站为中心3×3像元阵的IDDI。他们设定了7个能见度等级, IDDI数据与能见度对应, 对于每个间隔的能见度计算出IDDI均值和中值, 结果见图 3。能见度小于10 km被认为与大气沙尘有关, 能见度小于5 km与严重沙尘天气条件相关, IDDI与能见度的关系表明IDDI大于9.8对应于有沙尘情况, IDDI值大于14.6对应于严重沙尘天气。这些阈值的时空变化随着站点、季节和地域不同而改变, IDDI阈值10和15可以分别确定为对应的沙尘和严重沙尘条件。通常能见度10 km和5 km大约分别对应气溶胶的光学厚度为0.3~0.5和0.6~0.9。

图  3  地面能见度和红外差值沙尘指数对应关系^[6]

(其中能见度对应7个等级[0, 2.5], [2.5, 5], [5, 7.5], [7.5, 10], [10, 15], [15, 20], [20, 30] km; 括号内数字为样本数)

Figure  3.  Climatological relation between IDDI and visibility^[6]

(visibility measurements being classified into seven categories defined by the intervals [0, 2.5], [2.5, 5], [5, 7.5], [7.5, 10], [10, 15], [15, 20], [20, 30] km; the number inside bracket denotes the sample number)

下载: 全尺寸图片 幻灯片

IDDI计算的假设前提是地表温度特征在一定时段内没有发生改变, 而且也认为最大亮温的那一天大气是干洁的, 实际情况一般没有那么理想, 所以计算IDDI会存在系统的误差源, 主要表现为两点: ①在背景图像合成过程中, 如果最晴空的天气里还残留着有沙尘的信息, 会导致最后的IDDI值低估; ②由于冷空气的到来引起地面的突然大幅度降温, 这时用前10 d的合成背景亮温会高于当天的实际地表亮温, 这样会造成IDDI的高估。另外长时间持续云覆盖区最大亮温值依旧不是地表亮温, 或者在10 d合成亮温时出现逆温雾也会将逆温雾温度当成了地表亮温, 不过这两种情况往往是在温暖潮湿的中低纬度地区出现比较多, 对于干旱和半干旱地区出现这种情况的几率较少, 可以通过先验知识对其进行过滤。

2.   沙尘暴遥感监测数据处理

2.1   云检测

在沙尘判识之前, 云检测质量的好坏直接影响到沙尘的监测结果, 沙尘在相当程度上很容易被判识为云, 所以云检测在沙尘监测应用上要十分谨慎小心, 如果把沙尘判识成低云就会造成沙尘的漏判。同时如果云区的漏判会造成沙尘的错判, 它会把云造成的亮温衰减误认为是沙尘的贡献。云检测既可以直接使用业务的云检测产品, 也可以在本算法中单独进行云检测处理, 不过在使用业务云检测产品前必须对其进行严格检验, 尤其注意云检测产品中可能将沙尘判识为云。这里主要采用了单独的云检测处理, 首先利用可见光通道的反射率门槛值, 结合IDDI门槛, 主要检测白亮的整片云区。接下来进行红外窗区通道亮温的3×3像元矩阵偏差和IDDI的3×3像元矩阵偏差检测, 这对于检测云边缘比较有效。分裂窗通道亮温差T₁₁-T₁₂也做为云检测的重要依据, 因为通常卷云的分裂窗亮温差为正, 且亮温较低, 这对于检测高卷云比较有效。综合使用这些云检测算法, 给出云检测掩码, 在沙尘判识和IDDI计算时不云区进行处理。

2.2   背景参考图像生成

卫星观测的原始图像 (OI) 包含着地表和大气的辐射信息, 多天合成的背景参考图像 (RI) 应是干洁的, 且清晰地显示出地表的热辐射背景, 实时图像与背景图像相减得到差值图像 (DI), 显示的是变化大气的辐射影响。要分离OI图像中地表的稳定贡献和变化的大气负贡献。稳定的贡献是指在一段时间内 (10 d) 每个OI像元出现的最大亮温, 构建出晴空干洁 (热的) 参考图像 (RI)。用下式表示,

式 (3) 中, n为背景合成天数, i, j代表像元行列号, r代表P_ij^k中的最大值。图 4显示了某个时段04:25 (世界时, 下同) 的红外背景图像, 表示前10 d地表亮温值, 为后面的IDDI计算做准备。

图  4  GMS-5 10 d合成红外背景亮温图像

Figure  4.  GMS-5 integrated brightness temperature image of background during ten days

下载: 全尺寸图片 幻灯片

2.3   沙尘目标判识

沙尘目标判识是沙尘暴遥感监测自动算法的核心部分, 基于上述的物理原理, 结合静止卫星全部通道的多光谱观测信息进行多光谱聚类分析, 进行沙尘目标的判识。

最早的沙尘暴卫星监测方法就是采用多光谱聚类的方法, 主要采用AVH RR的可见光和红外两个通道进行聚类分析, 随着研究和应用的深入, 以及新型的多通道遥感器上天 (如MODIS), 逐渐采用更多光谱通道进行沙尘目标的卫星图像识别^[13-14]。研究者在使用原始卫星遥感通道的基础上, 开始利用通道组合来提取某些遥感特征量作为沙尘识别参量, 比如AVHRR的3, 4通道的亮温差、比值和分裂窗通道亮温差等。沙尘目标在不同的观测通道会显示出一个光谱聚类区, 这是单一目标识别过程, 而且沙尘在热红外分裂窗有非常明显的负亮温差特征, 结合沙尘的其他通道阈值就能比较准确的识别沙尘目标。采用已知的沙尘样本进行卫星图像的监督分类, 主要采用如下步骤来进行沙尘光谱参数设定。

在算法开发阶段, 通过提取静止气象卫星 (如GMS) 观测的沙尘样本进行统计分析, 并反复试验, 确立沙尘判识的判别因子和判别规则:分裂窗亮温差、可见光通道反射率、红外通道亮温、水汽通道亮温、IDDI值、可见光反射率与背景差值以及非云像元。通过多个样本统计分析建立起这些判别因子的阈值范围和权重, 再经过多次沙尘天气和非沙尘天气的监测进行检验和修正。多通道沙尘样本的直方图分析时, 还结合了其他非沙尘目标的联合分析, 最终确立一系列沙尘目标的判识门槛, 作为沙尘目标的判别规则。这一系列阈值形成一个参数文件, 在业务算法运行时直接读取。它们是判识沙尘目标的关键参数, 它随季节和观测时间的变化而有所不同, 另外类似遥感器通道的光谱响应函数差异也影响到这些阈值的选取, 只有选取精确可靠的阈值才能无误地利用卫星观测数据判识出沙尘目标。

3.   监测个例分析

采用上述的沙尘暴卫星遥感监测原理和方法, 对2001, 2002年3—5月GMS-5数据进行了处理, 并与地面气象观测站观测结果进行比较检验。彩图 5和彩图 6是2001年4月7日11:32 GMS-5沙尘暴监测结果图像, 它们分别表示沙尘暴监测结果的两种不同产品形式。彩图 5是沙尘暴监测结果格点图; 彩图 6是晴空区的IDDI结果产品, 它是沙尘暴监测产品的一个副产品, 该图是保存每个时次的无云区域的IDDI值, IDDI不仅是监测沙尘暴区域和面积, 而且还监测没有发生沙尘暴天气里, 大气中的微弱浮尘和气溶胶含量特征参量, 不过正如前面所描述, 该参数还不是一个完全的定量产品。

图  5  GSM-5沙尘暴监测结果格点图

（图中红点代表地面有沙尘报告的气象观测站）

Figure  5.  GSM-5 dust detection grid image

(red points denote the weather station with the report of dust occurrence)

下载: 全尺寸图片 幻灯片

图  6  GSM-5计算单天单时次的IDDI图

Figure  6.  GSM-5 IDDI image of single observation time

下载: 全尺寸图片 幻灯片

采用地面气象观测资料来验证卫星监测的准确性和可靠性。2001年4月7日沙尘暴天气过程检验的结果如彩图 5所示, 红点站表示有地面沙尘天气记录, 除了被云遮掩的区域, 其他区域都被精确判识出来。不仅如此, 自动监测算法还能较全面的监测到出现扬沙的区域。目前的算法可以实现从09:00—15:00沙尘连续的跟踪监测, 多个时次监测结果生成沙尘天气过程的演变动画。彩图 7显示了2001年4月中国北方IDDI月平均值, 它能清晰地表示这一时段沙尘发生的强度区域分布情况, 从中分析沙尘发生的源区和影响区域, 最明显的3个强沙尘区是塔克拉玛干沙漠、巴丹吉林—腾格里沙漠和毛乌素沙地。

图  7  GSM-5计算中国北方IDDI月平均值

Figure  7.  GSM-5 monthly mean IDDI in North China

下载: 全尺寸图片 幻灯片

在每天各个时次监测结果的基础上, 进行长时间的统计合成得到某一时段的沙尘发生频次统计图, 彩图 8显示2001年春季沙尘发生频次统计图, 它显示了这个春季沙尘发生的主要区域和影响路径。从多年的沙尘发生统计频次可以清楚分析沙尘发生的年际变化, 卫星遥感的优势就是能够得到比地面气象观测资料^[18-19]更精细的地域分布。

图  8  2001年春季GSM-5监测沙尘天气发生频次统计图

（图中色标表示相对频次大小）

Figure  8.  Frequency of dust occurence in spring of 2001 by GSM-5 detection

(color bar means relative frequency)

下载: 全尺寸图片 幻灯片

4.   小结与讨论

本文基于空气中沙尘在红外通道的两种响应特征, 结合传统的多通道光谱聚类方法, 开发了一套静止气象卫星监测我国北方沙尘暴业务算法, 该算法在国家卫星气象中心业务试验运行了近3年, 并不断补充完善, 现已在FY-2C遥感产品中投入业务运行, 成为一套客观自动的沙尘暴卫星遥感监测业务。从监测的结果分析, 它对于大的沙尘过程能够进行比较准确的监测, 并能充分发挥静止气象卫星高时间频次观测的优势, 监测时空变化较快的沙尘灾害天气过程。同时该系统输出的沙尘强度的半定量产品IDDI, 为沙尘气溶胶对局地和全球气候的影响提供了遥感数据。多时次长序列卫星监测结果有利于全面了解我国北方沙尘源区及其传播路径, 为沙尘暴治理决策提供客观科学依据。

红外通道的IDDI值和分裂窗通道亮温差不仅为沙尘遥感监测提供判据, 同时为空气中沙尘的定量参数反演提供了科学借鉴, 但要更好地解决这一问题, 需要精确知道沙尘的微物理参数 (复折射指数和粒子谱分布), 结合辐射传输计算模拟, 建立光学厚度查算表。沙尘定量遥感结果的真值检验是一个比较困难的工作, 中国气象局正在建设的沙尘暴监测预警服务系统中的太阳光度计观测数据是一个比较理想的地面真值资料, 不过严重沙尘暴发生时太阳光度计不能有效工作, 得不到可用的观测数据。

在实际的沙尘暴业务遥感监测过程中, 上述的原理方法尽管对大部分的沙尘目标有很好的跟踪监测效果, 但是笔者发现有少量的云和沙尘的区分存在一定困难, 通常的云检测算法也经常将沙尘归类为云目标, 有些薄卷云和强降水云会出现类似沙尘的分裂窗亮温差负值信号响应, 要准确区分这两类目标, 需要严格从云与沙的辐射机理及相关的先验知识去解决, 避免云与沙的误判。IDDI虽然可以作为沙尘强度的半定量参数, 但它受地表背景亮温和气象因子的影响不容忽视, 地表背景亮温在一天的日循环中存在变化 (正午过后最高, 子夜凌晨最低), 同样的沙尘强度在不同背景亮温情况下, 会形成不同IDDI值, 这就造成了一天中不同时间IDDI值的不可比性, 尤其正午与子夜之间的差别最大, Legrand^[6]对这种影响进行了比较全面的辐射传输模拟计算。另外分裂窗通道亮温差和IDDI两者都受大气中水汽量的影响, 大气中水汽量的差异会造成它们的不同, 而这种影响对于较大空间范围沙尘监测的阈值确定带来困难, 但在我国北方干旱和半干旱地区大气中水汽量含量不高, 水汽量变化也不剧烈。另外, 大气温度廓线和非沙尘类气溶胶的影响也是一个值得注意的影响因子。因此要更加准确监测和定量遥感沙尘, 需要进一步模拟分析以上影响因子, 并在实际遥感监测业务中引入相应的辅助参数, 这正是下一步工作的研究方向。