利用独立分量分析进行火山灰云遥感检测

李成范 戴羊羊 赵俊娟 尹京苑 周时强

(上海大学计算机工程与科学学院，上海 200444)

0 引言

一次大型的火山喷发往往能够形成巨大的火山灰云。火山灰云成分主要包括火山灰气体和火成岩碎屑两部分，其中火成岩碎屑通常是由直径＜2mm的火山碎屑物构成，火山灰气体主要是由 H2O、SO2、H2S、CO2等气体构成(Miffre et al.，2012;Donovan et al.，2013;Zaksek et al.，2013)。火成岩碎屑和火山灰气体混合能够在平流层中形成酸性气溶胶，消散周期较长，往往长达几周。此外，火山灰云与气象云特征比较接近，在高空中都呈现亮色调，不易区分，然而火山灰云的形成高度、光谱特征、组成成分等则与气象云具有明显的差异。大量的火山灰云覆盖在地球上空，不但削弱了到达地面的太阳辐射，引起臭氧层破坏、大气污染、温室效应、酸雨以及气温和降水异常等全球气候和环境系统的重大变化(Luke et al.，2001)，而且火山灰云的飘浮高度(一般处于平流层)恰好也是航空器飞行的高度，很容易腐蚀飞机发动机，导致飞机发动机熄火，引发航空安全事故，造成重大经济损失和人员伤亡(Mecikalski et al.，2007;Mccarthy et al.，2008;Prata，2009;尹京苑等，2013)。近年来，随着中国国际航线的大量开辟，中国飞行器遭遇火山灰云的可能性也随之增大。目前，如何对火山灰云进行检测已成为火山和安全领域科研工作者面前的重要任务之一。

遥感技术具有快速、空间覆盖范围广和时间分辨率高的特点(屈春燕等，2006;季灵运等，2009;许建东等，2009)，能够实时、准确地获取火山喷发形成火山灰云的海量信息，为火山灰云检测研究提供了有效的技术支撑。从20世纪80年代起，各国研究人员曾利用紫外吸收法、分裂窗亮温差法和多波段法等对遥感数据中的一些典型的火山灰云进行了检测研究(Andronico et al.，2009;Gangale et al.，2010;朱琳等，2011;Marzano et al.，2012)。然而，由于遥感数据存在波段相关性和数据冗余，在一定程度上限制了火山灰云的检测精度，这就需要在火山灰云遥感检测中引入新的数据处理方法。主成分分析(Principal Component Analysis，PCA)能够将大量复杂信息有效地压缩到少数几个主成分中，因此被广泛应用于数据压缩和高光谱领域。独立分量分析(Independent Component Analysis，ICA)是近年来发展起来的一种新的数据处理方法，能够将观测数据通过某种线性分解成为相互独立的特征成分，在遥感检测火山灰云中具有较好的应用潜力。在利用PCA预处理MODIS热红外波段数据的基础上，本文引入ICA方法用于遥感检测火山灰云信息，使得获取的各个成分信息之间相互独立，以便能够更加有效地表达不同成分信息和处理后续数据，从而提高了火山灰云遥感检测的精度。

1 研究方法和研究概况

利用ICA方法对研究区的火山灰云信息进行遥感检测可分以下4步进行(图1):1)利用PCA方法对MODIS传感器的热红外波段数据进行处理，将绝大多数有效信息压缩到少数几个主成分图像中;2)根据获取的火山灰云PC图像，分别对不同的MODIS热红外波段进行火山灰云敏感性分析，并从中选取适合于进行检测的波段;3)利用ICA方法进行火山灰云遥感检测，将火山灰云信息分离成相互独立的独立成分图像;4)利用图像分割检测出火山灰云信息。

图1 数据处理与方法流程Fig.1 Flow chart of data processing and methods.

1.1 独立分量分析方法

ICA作为一种新兴的信号分离技术，自1988年首次被提出以来，经过不断改进和完善，其应用已扩展到图像处理和模式识别等多个领域(Hyvärinen et al.，2000;Li et al.，2013)。ICA的计算公式为

式(1)中:X为可观测到的随机变量，A为未知混合矩阵，S为未知隐藏变量，且需要在尽可能少的假设条件下同时估计出未知混合矩阵A和向量s。

ICA具有如下假设:s是统计独立的，且为非高斯分布;s的数目与观测变量的数目相同，即混合矩阵是方阵。于是，在估计出未知混合矩阵A后，可以进一步计算A的逆矩阵W，则独立成分s为

对于标准ICA而言，当满足上述假设条件时，就能实现混合矩阵A和独立成分s的求解。ICA方法是在PCA基础上发展而来的。PCA通过将大量信息压缩到少数主成分中来实现数据降维和消除冗余，得到的主成分是不相关的，而ICA则得到的是尽可能相互独立的成分，且更加符合实际情况。与PCA相比，ICA具有更好的适应性，在火山灰云遥感检测中的优势更加明显。本文在PCA预处理MODIS热红外波段数据的基础上，利用ICA对火山灰云进行遥感检测，并将检测结果与美国地质调查局(USGS)标准光谱数据库和火山灰云SO2浓度分布进行对比分析。

1.2 数据源和研究区概况

搭载在Terra/aqua卫星上的MODIS传感器于1999年12月18日成功发射，共有36个离散波段，其中1～19波段和26波段为可见光、近红外波段，20～36波段为热红外波段，空间分辨率分别为250m、500m和1000m，最大扫描宽度为2330km;光谱覆盖了从可见光到热红外(0.4～1.44μm)的范围，且具有较高的时间分辨率，每天可过境4次，能够有效地对各种突发性灾害进行快速、实时监测。目前，MODIS数据已被广泛应用于大气污染、土地利用和生态环境变化以及地震与火山热红外异常监测等领域。

冰岛位于北大西洋中部，紧贴北极圈，面积约为10.3万km2，岛内冰川和火山广泛分布，有着“极圈火岛”的美誉。其中，位于首都雷克雅未克东南部的艾雅法拉火山在2010年3月20日晚喷发，接着在4月14日再次剧烈喷发，第2次喷发的强度和规模远大于前一次，并形成了巨大的火山灰云团(图2)。在风力的作用下，火山灰云扩散到欧洲大陆和北大西洋上空，造成欧洲大多数机场关闭，1000多个航班被取消，空中交通几近瘫痪，直接经济损失达上千亿欧元。

图2 艾雅法拉火山灰云的MODIS假彩色合成图像(a)与放大的火山灰云图像(b)Fig.2 The MODIS false color image of Eyjafjallajokull volcanic ash cloud(a)and the enlarged image of volcanic ash cloud(b).

2 遥感检测2010年艾雅法拉火山灰云

2.1 主成分分析处理

由PCA基本原理可知，PC1信息量最大，集中了大约80%的火山灰云信息，PC2次之，随后PC图像信息量依次递减。在本实验中，分别利用遥感图像处理软件ENVI4.6对MODIS传感器的热红外波段数据进行PCA处理，依次输入MODIS传感器的20～36热红外波段数据(第26波段除外)，分别得到相应的PC图像(图3)。在获得的PC图像中，仅有PC1，2，3，4中噪声信息较少，如图3所示，基本上能够满足遥感检测火山灰云信息的需要，在其余PC图像中，由于噪声信息较大或者基本上都是噪声信息，在本实验中忽略不计。

图3 PC图像，a，b，c，d分别为PC1～PC4图像Fig.3 PC images，a，b，c and d are the images of PC1，2，3 and 4，respectively.

从图3中看出，在PC1图像(图3a)中，火山灰云、气象云和陆地的边界对比最明显，其次是PC2图像(图3b)。PC1和PC2图像中火山口NE方向的小块火山灰云都被检测出来。PC3图像(图3c)和PC4图像(图3d)中火山灰云与气象云和陆地的对比相对较弱，也没有检测出火山口NE方向的小块火山灰云。据分析，PC图像中主要背景包括气象云和地表覆盖。对于气象云，云层越高，颜色越白，这是因为云层越高，其在热红外波段范围内的亮度温度越低，颜色就越白;相反，当云层较低时，由于吸收大量的地面辐射热量，其在热红外波段范围内亮度温度就较高，颜色也就较深。对于地表覆盖，地表亮度温度逐渐从浅灰色到黑色，浅灰色表明亮度温度较低，黑色表明亮度温度较高，并呈渐变关系。总的来说，PC1～4图像基本上增强了火山灰云与气象云和地表覆盖等背景的对比度，较好地突出了火山灰云信息，有利于后续的遥感检测研究。

2.2 火山灰云敏感性分析和波段选取

为了能够定量分析火山灰云PC1-4的图像效果，分别计算其可释方差、信噪比(SNR)，并根据火山灰云信息量递减的原则进行排列，结果如表1所示。

从表1中看出，PC1和PC2图像中火山灰云信息的可释方差分别为0.13和0.20，信噪比分别为1.66和1.33，火山灰云呈现出明显的亮色调，与地物背景对比最为明显，检测效果也最好。PC3和PC4图像的火山灰云信息量相对较少，可释方差分别为0.008和0.011，信噪比分别为4.55和1.58，火山灰云呈现出暗色调，与地物背景对比相对较弱。由此可见，并不是可释方差和信噪比的值越大，火山灰云对比效果就越好。

此外，为了验证MODIS热红外波段对火山灰云的敏感性，分别计算了不同热红外波段对PC1～4图像的贡献，结果如表2所示。

表1 PC图像对比Table1 Comparison of PC images

表2 MODIS热红外波段对PC图像的贡献(%)Table2 Contribution of MODIS thermal infrared bands to PC images(%)

从表2中看出，对于PC1图像，30波段的贡献率最大，其次是36波段;对于PC2图像，贡献率最大的36波段，其次是25波段;对于PC3图像，贡献率最大的是31和32波段，其他波段的贡献率则相对较小;对于PC4图像，贡献率最大的是29和32波段。综合来看，对于PC1～4图像，贡献最大的为36波段，总贡献率达到了72%，其次为31和30波段，分别达到了67%和65%。于是，经过对MODIS热红外波段进行火山灰云敏感性分析，最终选取36、31和30波段进行独立分量分析处理。

2.3 独立分量分析处理

分别将MODIS热红外波段36、31和30调入遥感图像处理软件ENVI4.6中进行ICA处理，即可得到相应的IC图像，结果如图4所示。

根据ICA基本原理，ICA处理后得到的IC图像是随机产生的，并不是根据火山灰云信息量的大小依次出现。因此，图4中的IC1(图4a)、IC2(图4b)和IC3(图4c)图像也是随机出现的。从图4(a～c)中看出，IC2图像的火山灰云信息量最大，影像最为清晰，且火山灰云与气象云、地表覆盖等背景的区分度最高，能够较好地反映出研究区内的火山灰云信息。此外，在IC2图像中，火山灰云呈现白色，纹理比较均匀，与周围呈现白色且纹理混乱的散碎云层对比明显，较好地突出了火山灰云的纹理和光谱特性。而在IC1和IC3图像中，火山灰云信息量较少，且影像相对模糊，与周围地物的区分度较小，不利于火山灰云的遥感检测。因此，主要选取IC2图像进行后续的火山灰云遥感检测。

图4 IC图像a，b，c分别为IC1，2，3图像，d为ICA检测出的火山灰云信息Fig.4 IC images a，b and c are the images of IC1，2 and 3，respectively，d is the detected volcanic ash cloud information using ICA.

2.4 图像分割和去噪处理

在得到火山灰云灰度IC图像之后，需要进行图像分割处理，以便将火山灰云信息检测出来。图像分割通常利用(X+kδ)来确定进行图像分割的下限和分割强度等级。其中，X为火山灰云信息的统计均值，δ为火山灰云信息的标准差，k为分割系数，通常取值范围在［1，3］之间。通过在遥感图像中选取100个火山灰云统计样本点，分别统计出IC图像中最大值、最小值、均值和标准差，结果如表3所示。经过多次验证，最终选取分割系数k=1、分割阈值为1.06时能够取得最佳的火山灰云信息分割效果。此外，为了消除分割后火山灰云中的少量气象云噪声，还需要对IC图像进行3×3滤波处理。经过图像分割和去噪处理后的火山灰云信息如图4d所示。

表3 IC灰度图像统计信息Table3 Statistical information of IC gray image

从图4d中看出，火山灰云的主体部分和火山口NE方向的小块火山灰云信息都被检测出来。这表明利用ICA方法既能有效地从MODIS热红外波段数据中检测出大片的火山灰云主体部分，又能检测出小块的火山灰云信息，在实际应用中优势较为明显。

3 火山灰云检测结果分析

冰岛艾雅法拉火山在2010年4月喷发形成规模巨大的火山灰云，并逐渐扩散，造成多个机场关闭，对全球环境和航空业产生的影响是近年来少见的。鉴于此，有关研究人员分别利用不同方法对火山灰云团进行了识别和动态追踪，并取得了一定的成果。图5a为德国研究人员Zaksek等(2013)从MODIS遥感数据中获得的艾雅法拉火山灰云信息，图5b为中国研究人员朱琳等(2011)从FY-3A遥感数据中获得的火山灰云信息。

图5 MODIS图像中检测出的火山灰云信息(a)(Zaksek et al.，2013)和FY-3A图像中检测出的火山灰云信息(b)(朱琳等，2011)Fig.5 The detected volcanic ash cloud information by MODIS image(a)(Zaksek et al.，2013)and the detected volcanic ash cloud information by FY-3A image(b)(Zhu et al.，2011).

由于不同研究中采用的地理坐标和投影系统往往是不一样的，因此从遥感数据中识别出的火山灰云形状也并不完全一致，甚至会存在一定的偏差和变形。从图5中看出，从MODIS和FY-3A数据都能够遥感检测出火山灰云信息的主体部分，且FY-3A数据还能更加准确地遥感检测出小块的火山灰云信息。综合图4d和图5，本实验利用ICA方法从MODIS热红外波段数据中检测出的火山灰云信息与已有的研究结果具有较好的一致性。但是，本实验中遥感检测出的火山灰云分布区域相对较窄，据分析这主要是由于分布在外围的火山灰云信息浓度较为稀薄而没有被ICA有效地检测出来造成的。

此外，为了定量地验证ICA方法检测火山灰云的准确性，下面将分别从光谱匹配度和火山灰云SO2浓度两方面进行分析。

3.1 光谱匹配度分析

根据MODIS数据中各波段的波长区间，将USGS标准光谱数据库中的反射率数据重采样至MODIS数据相对应的波段区间，并建立其在MODIS数据中的光谱曲线。通过将大气校正后的光谱曲线与重采样后的火山灰矿物标准光谱库进行相似性匹配(图6)。结果表明，火山灰的光谱曲线与光谱库中的火山灰的匹配程度达到83.55%，取得了较好的效果。

ICA方法能够降低遥感数据不同波段之间的相关性和冗余，通过分离出不同的地物特征信息，有利于火山灰云的遥感检测，并且检测出的火山灰云也更加贴近实际分布情况。

3.2 火山灰云SO2浓度分布分析

SO2是目前较为常用的火山灰云遥感检测指标之一。虽然火山类型和组成成分多种多样，但是无论何种类型火山喷发时，总会含有一定数量的SO2气体成分，这就使得火山灰云的SO2浓度要明显高于普通大气中的含量。

采用2010年4月19日艾雅法拉火山爆发期间的OMI卫星遥感数据，利用波段残差差分(Band Residual Difference，BRD)法计算 SO2浓度。其基本原理是利用 SO2在波长310.8～314.4nm范围内的吸收波峰与波谷位置(310.8nm、311.9nm、313.2nm、314.4nm)来计算SO2有效信息。然后，对获取的SO2柱浓度数据进行预处理，亦即利用数据属性文件中的Quality Flages参数和Radiative Cloud Fraction参数信息去掉高云量像元(Radiative Cloud Fraction＞0.2)，进而得到有效的SO2浓度像元信息(图7)。

图6 光谱曲线匹配结果Fig.6 The result of spectral curve matching of volcanic ash cloud.

图7 火山灰云SO2浓度分布Fig.7 SO2concentration distribution of volcanic ash cloud.

从图7中看出，火山灰云的SO2浓度分布是从北向南、从中心向外逐渐降低的。据分析，这可能是由于火山灰刚喷发进入高空时，携带了大量的热量和一定量的SO2，使得大气对流层中SO2的浓度相对较为积聚，要明显地高于周围的SO2浓度造成的。而且火山喷发初期的SO2浓度的扩散较慢，随着时间的延长，在风力作用下，扩散速度加快，SO2浓度迅速降低，引起SO2浓度分布与火山灰云的变化趋势较为相似。通过对比图5和图6可知，火山灰云SO2浓度分布与检测出的火山灰云信息存在较高的一致性。

4 结论

利用遥感技术能够实时检测和追踪火山灰云变化信息，为防灾减灾、生态环境和航空安全领域提供预警和技术支持。本研究提出在对MODIS数据进行PCA处理的基础上，利用ICA方法对2010年4月19日形成的冰岛艾雅法拉火山灰云进行遥感检测，并引入光谱匹配度和SO2浓度分布进行验证。结果表明，ICA既能从MODIS图像中检测出火山灰云信息的主体部分，又能有效地检测出小块区域的火山灰云信息，获取的火山灰光谱特征与USGS标准光谱数据库的光谱匹配度达到83.55%，且与火山灰云SO2浓度分布和相关研究成果具有较好的一致性。

尽管本实验取得了较好的火山灰云遥感检测效果，但是仍然存在一些不足之处:1)在利用ICA方法进行火山灰云遥感检测时获取的独立成分出现的次序是随机的，并不是根据火山灰云信息量的大小依次出现的;2)ICA方法从遥感数据中仅检测出了浓度较大的火山灰云信息，而对于浓度较为稀薄的外围火山灰云信息则并未检测出来。目前，这些问题还都很少涉及，这需要在今后的研究中进行更加深入的探讨。

季灵运，许建东，林旭东，等.2009.利用卫星热红外遥感技术监测长白山天池火山活动性［J］.地震地质，31(4):617—627.doi:10.3969/j.issn.0253-4967.2009.04.005.

JI Ling-yun，XU Jian-dong，LIN Xu-dong，et al.2009.Application of satellite thermal infrared remote sensing in monitoring Changbaishan Tianchi volcano activity［J］.Seismology and Geology，31(4):617—627(in Chinese).

屈春燕，单新建，马瑾.2006.卫星热红外遥感在火山活动性监测中的应用［J］.地震地质，28(1):99—110.

QU Chun-yan，SHAN Xin-jian，MA Jin.2006.Application of satellite thermal infrared remote sensing in detection of volcano activity［J］.Seismology and Geology，28(1):99—110(in Chinese).

许建东，栾鹏，樊笑英，等.2009.基于遥感影像光谱与纹理分析的地物分类:以长白山天池火山地区为例［J］.地震地质，31(4):607—616.doi:10.3969/j.issn.0253-4967.2009.04.004.

XU Jian-dong，LUAN Peng，FAN Xiao-ying，et al.2009.Analysis of spectrum and texture information on Changbaishan Tianchi volcano caldera and its application［J］.Seismology and Geology，31(4):607—616(in Chinese).

尹京苑，沈迪，李成范.2013.卫星遥感技术在火山灰云监测中的应用［J］.地震地质，35(2):347—362.doi:10.3969/j.issn.0253-4967.2013.02.013.

YIN Jing-yuan，SHEN Di，LI Cheng-fan.2013.Application of satellite remote sensing in volcanic ash cloud monitoring［J］.Seismology and Geology，35(2):347—362(in Chinese).

朱琳，刘健，刘诚，等.2011.复杂气象条件下火山灰云遥感监测新方法［J］.中国科学(D辑)，41(7):1029—1036.

ZHU Lin，LIU Jian，LIU Cheng，et al.2011.Satellite remote sensing of volcanic ash cloud in complicated meteorological conditions［J］.Science in China(Ser D)，41(7):1029—1036(in Chinese).

Andronico D，Spinetti C，Cristaldi A，et al.2009.Observations of Mt.Enta volcanic ash plumes in 2006:An integrated approach from ground-based and polar satellite Noaa-AVHRR monitoring system ［J］.Journal of Volcanology and Geothermal Research，180(2-4):135—147.

Donovan D P，Apitley A.2013.Practical depolarization-ratio-based inversion procedure:Lidar measurements of theEyjafjllajokull ash cloud over the Netherlands［J］.Applied Optics，52(11):2394—2415.

Gangale G，Prata A J，Clarisse L.2010.The infrared spectral signature of volcanic ash determined from high-spectral resolution satellite measurements［J］.Remote Sensing of Environment，114(2):414—425.

Hyvärinen A，Oja E.2000.Independent component analysis:Algorithm and application［J］.Neural Networks，13(4-5):411—430.

Li C F，Yin J Y.2013.Variational Bayesian independent component analysis-support vector machine for remote sensing classification［J］.Computers and Electrical Engineering，39(3):717—726.

Luke P F，Andrew J L H，Robert W.2001.Improved identification of volcanic features using Landsat 7 ETM+［J］.Remote Sensing of Environment，78(1-2):180—193.

Marzano F S，Picciotti E，Vulpiani G，et al.2012.Synthetic signatures of volcanic ash cloud particles from X-band dual-polarization radar［J］.IEEE Transactions on Geoscience and Remote Sensing，50(1):193—211.

Mccarthy E B，Bluth G J S，Watson I M，et al.2008.Detection and analysis of the volcanic clouds associated with the 18 and 28 August 2000 eruption of Miyakejima volcano，Japan［J］.International Journal of Remote Sensing，29(22):6597—6620.

Mecikalski J R，Feltz W F，Murray J J，et al.2007.Aviation applications for satellite-based observations of cloud properties，convection initiation，in-flight icing，turbulence，and volcanic ash ［J］.Bulletin of the American Meteorological Society，88(10):1589—1607.

Miffre A，Davic G，Thomas B，et al.2012.Interpretation of accurate UV polarization lidar measurements:Application to volcanic ash number concentration retrieval［J］.Journal of Atmospheric and Oceanic Technology，29(4):558—568.

Prata A J.2009.Satellite detection of hazardous volcanic clouds and the risk to global air traffic［J］.Natural Hazards，51(2):303—324.

Zaksek K，Hort M，Zaleteij J，et al.2013.Monitoring volcanic ash cloud top height through simultaneous retrieval of optical data from orbiting and geostationary satellite ［J］.Atmospheric Chemistry and Physics，13(5):2589—2606.