李成范 戴羊羊 趙俊娟 尹京苑 周時強

(上海大學計算機工程與科學學院，上海 200444)

0 引言

一次大型的火山噴發(fā)往往能夠形成巨大的火山灰云。火山灰云成分主要包括火山灰氣體和火成巖碎屑兩部分，其中火成巖碎屑通常是由直徑＜2mm的火山碎屑物構成，火山灰氣體主要是由 H2O、SO2、H2S、CO2等氣體構成(Miffre et al.，2012;Donovan et al.，2013;Zaksek et al.，2013)?；鸪蓭r碎屑和火山灰氣體混合能夠在平流層中形成酸性氣溶膠，消散周期較長，往往長達幾周。此外，火山灰云與氣象云特征比較接近，在高空中都呈現亮色調，不易區(qū)分，然而火山灰云的形成高度、光譜特征、組成成分等則與氣象云具有明顯的差異。大量的火山灰云覆蓋在地球上空，不但削弱了到達地面的太陽輻射，引起臭氧層破壞、大氣污染、溫室效應、酸雨以及氣溫和降水異常等全球氣候和環(huán)境系統(tǒng)的重大變化(Luke et al.，2001)，而且火山灰云的飄浮高度(一般處于平流層)恰好也是航空器飛行的高度，很容易腐蝕飛機發(fā)動機，導致飛機發(fā)動機熄火，引發(fā)航空安全事故，造成重大經濟損失和人員傷亡(Mecikalski et al.，2007;Mccarthy et al.，2008;Prata，2009;尹京苑等，2013)。近年來，隨著中國國際航線的大量開辟，中國飛行器遭遇火山灰云的可能性也隨之增大。目前，如何對火山灰云進行檢測已成為火山和安全領域科研工作者面前的重要任務之一。

遙感技術具有快速、空間覆蓋范圍廣和時間分辨率高的特點(屈春燕等，2006;季靈運等，2009;許建東等，2009)，能夠實時、準確地獲取火山噴發(fā)形成火山灰云的海量信息，為火山灰云檢測研究提供了有效的技術支撐。從20世紀80年代起，各國研究人員曾利用紫外吸收法、分裂窗亮溫差法和多波段法等對遙感數據中的一些典型的火山灰云進行了檢測研究(Andronico et al.，2009;Gangale et al.，2010;朱琳等，2011;Marzano et al.，2012)。然而，由于遙感數據存在波段相關性和數據冗余，在一定程度上限制了火山灰云的檢測精度，這就需要在火山灰云遙感檢測中引入新的數據處理方法。主成分分析(Principal Component Analysis，PCA)能夠將大量復雜信息有效地壓縮到少數幾個主成分中，因此被廣泛應用于數據壓縮和高光譜領域。獨立分量分析(Independent Component Analysis，ICA)是近年來發(fā)展起來的一種新的數據處理方法，能夠將觀測數據通過某種線性分解成為相互獨立的特征成分，在遙感檢測火山灰云中具有較好的應用潛力。在利用PCA預處理MODIS熱紅外波段數據的基礎上，本文引入ICA方法用于遙感檢測火山灰云信息，使得獲取的各個成分信息之間相互獨立，以便能夠更加有效地表達不同成分信息和處理后續(xù)數據，從而提高了火山灰云遙感檢測的精度。

1 研究方法和研究概況

利用ICA方法對研究區(qū)的火山灰云信息進行遙感檢測可分以下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是統(tǒng)計獨立的，且為非高斯分布;s的數目與觀測變量的數目相同，即混合矩陣是方陣。于是，在估計出未知混合矩陣A后，可以進一步計算A的逆矩陣W，則獨立成分s為

對于標準ICA而言，當滿足上述假設條件時，就能實現混合矩陣A和獨立成分s的求解。ICA方法是在PCA基礎上發(fā)展而來的。PCA通過將大量信息壓縮到少數主成分中來實現數據降維和消除冗余，得到的主成分是不相關的，而ICA則得到的是盡可能相互獨立的成分，且更加符合實際情況。與PCA相比，ICA具有更好的適應性，在火山灰云遙感檢測中的優(yōu)勢更加明顯。本文在PCA預處理MODIS熱紅外波段數據的基礎上，利用ICA對火山灰云進行遙感檢測，并將檢測結果與美國地質調查局(USGS)標準光譜數據庫和火山灰云SO2濃度分布進行對比分析。

1.2 數據源和研究區(qū)概況

搭載在Terra/aqua衛(wèi)星上的MODIS傳感器于1999年12月18日成功發(fā)射，共有36個離散波段，其中1～19波段和26波段為可見光、近紅外波段，20～36波段為熱紅外波段，空間分辨率分別為250m、500m和1000m，最大掃描寬度為2330km;光譜覆蓋了從可見光到熱紅外(0.4～1.44μm)的范圍，且具有較高的時間分辨率，每天可過境4次，能夠有效地對各種突發(fā)性災害進行快速、實時監(jiān)測。目前，MODIS數據已被廣泛應用于大氣污染、土地利用和生態(tài)環(huán)境變化以及地震與火山熱紅外異常監(jiān)測等領域。

冰島位于北大西洋中部，緊貼北極圈，面積約為10.3萬km2，島內冰川和火山廣泛分布，有著“極圈火島”的美譽。其中，位于首都雷克雅未克東南部的艾雅法拉火山在2010年3月20日晚噴發(fā)，接著在4月14日再次劇烈噴發(fā)，第2次噴發(fā)的強度和規(guī)模遠大于前一次，并形成了巨大的火山灰云團(圖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圖像中主要背景包括氣象云和地表覆蓋。對于氣象云，云層越高，顏色越白，這是因為云層越高，其在熱紅外波段范圍內的亮度溫度越低，顏色就越白;相反，當云層較低時，由于吸收大量的地面輻射熱量，其在熱紅外波段范圍內亮度溫度就較高，顏色也就較深。對于地表覆蓋，地表亮度溫度逐漸從淺灰色到黑色，淺灰色表明亮度溫度較低，黑色表明亮度溫度較高，并呈漸變關系?？偟膩碚f，PC1～4圖像基本上增強了火山灰云與氣象云和地表覆蓋等背景的對比度，較好地突出了火山灰云信息，有利于后續(xù)的遙感檢測研究。

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圖像的火山灰云信息量最大，影像最為清晰，且火山灰云與氣象云、地表覆蓋等背景的區(qū)分度最高，能夠較好地反映出研究區(qū)內的火山灰云信息。此外，在IC2圖像中，火山灰云呈現白色，紋理比較均勻，與周圍呈現白色且紋理混亂的散碎云層對比明顯，較好地突出了火山灰云的紋理和光譜特性。而在IC1和IC3圖像中，火山灰云信息量較少，且影像相對模糊，與周圍地物的區(qū)分度較小，不利于火山灰云的遙感檢測。因此，主要選取IC2圖像進行后續(xù)的火山灰云遙感檢測。

圖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為火山灰云信息的統(tǒng)計均值，δ為火山灰云信息的標準差，k為分割系數，通常取值范圍在［1，3］之間。通過在遙感圖像中選取100個火山灰云統(tǒng)計樣本點，分別統(tǒng)計出IC圖像中最大值、最小值、均值和標準差，結果如表3所示。經過多次驗證，最終選取分割系數k=1、分割閾值為1.06時能夠取得最佳的火山灰云信息分割效果。此外，為了消除分割后火山灰云中的少量氣象云噪聲，還需要對IC圖像進行3×3濾波處理。經過圖像分割和去噪處理后的火山灰云信息如圖4d所示。

表3 IC灰度圖像統(tǒng)計信息Table3 Statistical information of IC gray image

從圖4d中看出，火山灰云的主體部分和火山口NE方向的小塊火山灰云信息都被檢測出來。這表明利用ICA方法既能有效地從MODIS熱紅外波段數據中檢測出大片的火山灰云主體部分，又能檢測出小塊的火山灰云信息，在實際應用中優(yōu)勢較為明顯。

3 火山灰云檢測結果分析

冰島艾雅法拉火山在2010年4月噴發(fā)形成規(guī)模巨大的火山灰云，并逐漸擴散，造成多個機場關閉，對全球環(huán)境和航空業(yè)產生的影響是近年來少見的。鑒于此，有關研究人員分別利用不同方法對火山灰云團進行了識別和動態(tài)追蹤，并取得了一定的成果。圖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).

由于不同研究中采用的地理坐標和投影系統(tǒng)往往是不一樣的，因此從遙感數據中識別出的火山灰云形狀也并不完全一致，甚至會存在一定的偏差和變形。從圖5中看出，從MODIS和FY-3A數據都能夠遙感檢測出火山灰云信息的主體部分，且FY-3A數據還能更加準確地遙感檢測出小塊的火山灰云信息。綜合圖4d和圖5，本實驗利用ICA方法從MODIS熱紅外波段數據中檢測出的火山灰云信息與已有的研究結果具有較好的一致性。但是，本實驗中遙感檢測出的火山灰云分布區(qū)域相對較窄，據分析這主要是由于分布在外圍的火山灰云信息濃度較為稀薄而沒有被ICA有效地檢測出來造成的。

此外，為了定量地驗證ICA方法檢測火山灰云的準確性，下面將分別從光譜匹配度和火山灰云SO2濃度兩方面進行分析。

3.1 光譜匹配度分析

根據MODIS數據中各波段的波長區(qū)間，將USGS標準光譜數據庫中的反射率數據重采樣至MODIS數據相對應的波段區(qū)間，并建立其在MODIS數據中的光譜曲線。通過將大氣校正后的光譜曲線與重采樣后的火山灰礦物標準光譜庫進行相似性匹配(圖6)。結果表明，火山灰的光譜曲線與光譜庫中的火山灰的匹配程度達到83.55%，取得了較好的效果。

ICA方法能夠降低遙感數據不同波段之間的相關性和冗余，通過分離出不同的地物特征信息，有利于火山灰云的遙感檢測，并且檢測出的火山灰云也更加貼近實際分布情況。

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

SO2是目前較為常用的火山灰云遙感檢測指標之一。雖然火山類型和組成成分多種多樣，但是無論何種類型火山噴發(fā)時，總會含有一定數量的SO2氣體成分，這就使得火山灰云的SO2濃度要明顯高于普通大氣中的含量。

采用2010年4月19日艾雅法拉火山爆發(fā)期間的OMI衛(wèi)星遙感數據，利用波段殘差差分(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濃度分布是從北向南、從中心向外逐漸降低的。據分析，這可能是由于火山灰剛噴發(fā)進入高空時，攜帶了大量的熱量和一定量的SO2，使得大氣對流層中SO2的濃度相對較為積聚，要明顯地高于周圍的SO2濃度造成的。而且火山噴發(fā)初期的SO2濃度的擴散較慢，隨著時間的延長，在風力作用下，擴散速度加快，SO2濃度迅速降低，引起SO2濃度分布與火山灰云的變化趨勢較為相似。通過對比圖5和圖6可知，火山灰云SO2濃度分布與檢測出的火山灰云信息存在較高的一致性。

4 結論

利用遙感技術能夠實時檢測和追蹤火山灰云變化信息，為防災減災、生態(tài)環(huán)境和航空安全領域提供預警和技術支持。本研究提出在對MODIS數據進行PCA處理的基礎上，利用ICA方法對2010年4月19日形成的冰島艾雅法拉火山灰云進行遙感檢測，并引入光譜匹配度和SO2濃度分布進行驗證。結果表明，ICA既能從MODIS圖像中檢測出火山灰云信息的主體部分，又能有效地檢測出小塊區(qū)域的火山灰云信息，獲取的火山灰光譜特征與USGS標準光譜數據庫的光譜匹配度達到83.55%，且與火山灰云SO2濃度分布和相關研究成果具有較好的一致性。

盡管本實驗取得了較好的火山灰云遙感檢測效果，但是仍然存在一些不足之處:1)在利用ICA方法進行火山灰云遙感檢測時獲取的獨立成分出現的次序是隨機的，并不是根據火山灰云信息量的大小依次出現的;2)ICA方法從遙感數據中僅檢測出了濃度較大的火山灰云信息，而對于濃度較為稀薄的外圍火山灰云信息則并未檢測出來。目前，這些問題還都很少涉及，這需要在今后的研究中進行更加深入的探討。

季靈運，許建東，林旭東，等.2009.利用衛(wèi)星熱紅外遙感技術監(jiān)測長白山天池火山活動性［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.衛(wèi)星熱紅外遙感在火山活動性監(jiān)測中的應用［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.基于遙感影像光譜與紋理分析的地物分類:以長白山天池火山地區(qū)為例［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.衛(wèi)星遙感技術在火山灰云監(jiān)測中的應用［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.復雜氣象條件下火山灰云遙感監(jiān)測新方法［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.