EMD与分形相结合的遥感影像水体信息提取方法

周林滔， 杨国范， 赵福强， 杜娟

(沈阳农业大学水利学院,沈阳 110866)

0 引言

近年来遥感数据量正呈几何级数增长，如何快速有效地从海量数据中提取有效信息已成为遥感应用的关键。目前国内外很多学者在研究遥感影像水域目标提取方面时，主要使用的方法多基于地物光谱特征。如Barton等[1]使用AVHRR的第2波段与第1波段的比值识别水体; Steven等[2]采用非监督分类方法将地物分为10类，然后再聚类为水体和陆地2类; 陈蕾等[3]利用不同水质的光谱反射率差别较大的特征提取水体信息,有效地将各种水质类型的水体与山体阴影区分开来; 宋启帆等[4]通过多种融合试验认为采用归一化差异水体指数法提取水体的精度最高; 郭振亚等[5]提出了将各波段灰度值加权并设定阈值提取水体信息的方法，有效避免了阴影干扰。像元的亮度值代表其中地物的平均辐射值，随地物的成分、纹理及形态等的变化而变化，因此实际应用中通常会遇到地物类别光谱反射率混合的问题。空间特征相结合的研究逐步得到重视。如Argialas等[6]提出了基于面向对象的分类方法，综合利用空间位置和光谱特征提取了水域信息； 都金康等[7]提出采用决策树方法提取水体信息，结合空间特征信息对其进行分类。纹理特征是地物空间特征的重要组成，目前研究中分形方法的应用尚不成熟，大多还停留在可行性分析试验阶段，尚未达到对整幅影像应用分形方法提取的程度。

本文研究的主要问题为： 一是如何有效避免“同物异谱”现象； 二是如何有效结合光谱特征和纹理特征提取水体信息。本文尝试使用经验模态分解(empirical mode decomposition,EMD)的方法分解光谱信息，同时使用分形理论计算图像每个像元的分形维数，生成分维图，从而综合利用光谱信息和纹理信息提高地物识别精度。EMD分解可将图像信号分解为一组“由细到粗”的细节信息和一个大尺度趋势信息,每个尺度都具有明显的物理意义，即信号在某一频率和振幅范围内的特征信息。通常大尺度分量代表影像的整体轮廓特征,小尺度分量反映影像的纹理特征等，有助于增强地物特征差异性。像元的分形维数可反映影像的平滑程度，一般分形维数小的区域影像较平滑，分形维数大的区域影像变化越大，通常为物体的边缘[8]。

1 研究方法

1.1 EMD影像分解

EMD是由美国NASA提出的一种信号分析方法[9]。该方法不需要事先选择基函数，而是根据信号本身特征将原始信号分解为有限个本征模态函数(instrinsic mode function,IMF)和一个趋势项rn，每个IMF是信号在不同频率和不同振幅下的分量信号，能够很好地反映信号在不同尺度下局部频率特征，这与傅里叶分解和小波分解方法有着本质性的差别[10-11]，可表示为

(1)

式中：I为原始信号；IMFj为第j个模态分量；rn为单调残差函数；m为模态分量个数；t为时间。

每个本征模态函数需满足如下条件之一[12-13]： ①在整个信号长度上，一个IMF的极值点和过零点数目必须相等或至多只相差一点； ②由所有极大值点构成的包络线和由极小值点构成的包络线的均值和接近于零，即IMF的上下包络线对称于时间轴。

EMD分解图像实现过程有以下3个步骤。

1)找到信号的极大值和极小值，即图像曲面所有局部极值点。

2)对所有的极大值和极小值分别进行曲面拟合，经插值后分别得到极大值点和极小值点曲面的包络面Emax和Emin，包络面的均值为

(2)

3)设原始信号为I(t)，则h1=I(t)-Emean。理论上h1为第一个模态分量IMF1，但不一定满足上述条件，再将h1作为原始信号，重复上述步骤，直到满足IMF条件或得到终止条件，即

(3)

式中：SD为标准差；M，N分别为图像x，y方向的像元个数；k为分解出的二维固态模函数(bidimensional intrinsic mode function,BIMF)个数。

重复步骤1—3最终得到IMF1。第一层残差Ires=I-IMF1，将残差作为原始信号进行分解可得到IMF2，以此类推，得到所有模态分量。

1.2 分形维数计算

本文借鉴朱骥等[14]提出的单个像元分维数的计算方法，得到整幅影像的分形图。以任一像元为中心，分别将横向、垂向、东北向及西北向共4个方向的分维数取平均值后作为该像元的分维数(图1)。

图1 分形维数方向维

选定窗口大小L，任选p[0,L],令r=b-a,0≤a

(4)

式中：x，y为像元的平面坐标；f为相应的灰度值。遍历所有可能的a值和对应的b值计算n(p)的平均值，即覆盖L像元×L像元窗口的盒子数为

(5)

2 模型建立

使用EMD方法和分形理论建立结合光谱信息和纹理特征的水体信息提取模型，主要步骤包括： ①对环境小卫星的4波段数据做主成分(principal component,PC)分析，使用降噪后的第一主分量； ②对第一主成分作EMD分解，获取前3个尺度的模态函数IMF1，IMF2，IMF3； ③使用上述分形维数计算方法对第一主分量每个像元的分维数进行计算，生成分形图； ④在ENVI中使用“LayerStack”工具波段组合原始4波段、IMF1，IMF2，IMF3及分形图； ⑤使用极大似然法分类器将研究区内地物分为水体、居民地及其他未分类地物3类。

2.1 EMD分解

现将2012年9月18日第一主分量影像进行二维EMD分解，结果如图2所示。细节信息包含影像的纹理特征等，有助于增强地物特征差异性。

(a) 原始影像(b) IMF1(c) IMF2(d) IMF3

从图2得知，原始影像为第一主分量影像，左下方黑色区域为柳河彰武地区局部河段；IMF1表征原始影像的最小尺度细节信息；IMF2表征次小尺度细节信息，以此类推，大尺度分量代表影像的整体轮廓特征。

2.2 窗口大小对分形维数计算的影响

已有研究认为，分形维数算法是n2型的[15]，即运算时间与窗口大小的平方成正比。因此窗口大小的选择非常重要，在保证回归方程有较好的直线性的同时也要兼顾算法的时间复杂度。

本文分别使用4像元×4像元、5像元×5像元、8像元×8像元及10像元×10像元大小的窗口进行分形维数的计算试验(图3)，发现分形维数值随窗口增大而减小，这与朱骥等[14]研究的结论相符。小窗口计算得到的分维值图能够反映地物的细节特征。大窗口计算得到的分形图能够反映地物大范围的整体特征。本文最终选择窗口大小为5像元×5像元。

图3 窗口大小为5像元×5像元(左)和10像元×10像元(右)的分形图

3 模型应用

3.1 研究区概况与数据源

研究区位于辽宁省西北部，科尔沁沙地的南部，E121°31′48″～121°34′48″，N42°4′12″～42°30′36″之间，年降水量510 mm，平均相对湿度61%，年平均气温7.2℃，属于半干旱地区。流经研究区的柳河为辽河的一级支流，发源于内蒙古自治区，于新民市汇入辽河，年平均地表径流量为3.28亿m3/a。研究区内的柳河河段位于闹得海水库下游，丰水期平均水深约0.5 m，枯水期平均水深约0.1 m。

遥感影像数据来源于环境减灾小卫星HJ-1A的CCD1相机，空间分辨率为30 m，探测谱段范围为蓝、绿、红、近红外波段，时相为2012年9月18日。本研究用的影像数据经过几何纠正及大气校正等预处理，制成了4(R)3(G)2(B)合成图像(图4)，共含1 673像元×1 234像元，包括彰武县县城、柳河河道及坑塘水库等地物。图4中北缘白色斑块为厚云，会对图像分类结果造成影响，但由于该处不存在水域，因此不影响水体信息提取结果。

图4 研究区示意图

3.2 精度评价

利用上述模型提取水体信息，包括柳河河道和坑塘水库等水体。其中柳河河道平均水深很浅，年均水深最大值约0.5 m，最小值约0.1 m。由试验可知，其光谱特征与居民地极为相似，如图5所示。

图6分别给出利用极大似然法、EMD极大似然法、EMD分形极大似然法及其聚类、过滤处理后得到的分类结果。

图5 水体与居民地的波谱曲线

(a) 极大似然法(b) EMD极大似然法(c) EMD分形极大似然法(d) 对(c)聚类、过滤处理

比较图6可以看出，单纯使用极大似然法分类器提取水体信息时，其绝大部分被误分类为居民地，同时未分类地物被误分为居民地的情况也很严重(图6(a))； 使用EMD方法后上述情况有所改善，但仍存在水体信息被误分为居民地的情况(图6(b))； 加入分形纹理特征后该情况得到有效改善，总体的分类精度得到提高，但河流水体的连续性却不如之前，不连续处为漏提取部分，这可能与部分河道周边像元分维数小、图像变化不突出有关； 此外，居民地被误分类为水体的情况较少(图6(c))。对EMD分形极大似然法的分类结果(图6(c))进行聚类处理和过滤处理，最终得到图6(d)中连续性较好的河道水体信息，水体误分减少，地物分类更加清晰。这表明，本文提出的利用结合光谱信息和纹理信息的水体提取模型可以较好地避免将河流水体或未分类地物误分为居民地，从而提高水体信息的自动提取精度。

为定量评价EMD分形-极大似然法与极大似然法提取精度，使用验证样本建立混淆矩阵，分别计算2种方法提取水体信息的生产者精度、用户精度以及Kappa，并将其进行比较，如表1所示。

表1 水体的提取精度比较

由表1看出，EMD分形-极大似然法精度比传统监督分类极大似然法精度有了明显的提高，其中生产者精度分别提高了24.39%，用户精度提高了11.47%，Kappa提高了0.129 9。

4 结论

1)本文建立的水体信息提取模型充分结合了地物光谱特征和纹理特征。结合EMD算法重组波谱的模态分量后增强了光谱特征的差异性，从而有效避免了由于“同物异谱”引起的误分现象。依据分形理论，通过大量试验选定了合适的窗口大小，计算像元的分维数，较好地提取了地物的纹理特征。利用本文建立的模型可以在保证提取精度的前提下，有效地提高水资源调查、监测及保护的效率。

2)EMD分解方法理论上具有通用性，适用于解决各种遥感数据源中的“同物异谱”现象。本文研究发现，EMD极大似然法加入分形维数后虽然能够有效地避免将水体误分为居民地等情况，但同时会在不同程度上影响水体提取的完整性，其具体原因和改进方法还有待进一步研究。

参考文献(References)：

[1] Barton I J,Bathols J M.Monitoring floods with AVHRR[J].Remote Sensing of Environment,1989,30(1):89-94.

[2] Steven M K,Patrick L,Brezonik L G,et al.A procedure for regional lake water clarity assessment using Landsat multispectral data[J].Remote Sensing of Environment,2002,82(1):38-47.

[3] 陈蕾,邓孺孺,陈启东,等.基于水质类型的TM图像水体信息提取[J].国土资源遥感,2012,24(1):90-94.

Chen L,Deng R R,Chen Q D,et al.The extraction of water body information from TM imagery based on water quality types[J].Remote Sensing for Land and Resources,2012,24(1):90-94.

[4] 宋启帆,王少军,张志,等.基于WorldView II图像的钨矿区水体信息提取方法研究——以江西大余县为例[J].国土资源遥感,2011,23(2):33-36.

Song Q F,Wang S J,Zhang Z,et al.A water information extraction method based on WorldView II Remote sensing image in tungsten ore districts:A case study of Dayu County in Jiangxi Province[J].Remote Sensing for Land and Resources,2011,23(2):33-36.

[5] 郭振亚,王心源,王传辉,等.巢湖流域水体信息提取方法研究[J].遥感技术与应用,2012,27(3):443-448.

Guo Z Y,Wang X Y,Wang C H,et al.The research on extraction method of water body information in Chaohu lake basin[J].Remote Sensing Technology and Application,2012,27(3):443-448.

[6] Argialas D,Tzotsos A.Automatic extraction of physiographic features and alluvial fans in Nevada,USA from digital elevation models and satellite imagery through multiresolution segmentation and object-oriented classification[C]//Proceedings of Asprs 2006 Annual Conference,2006.

[7] 都金康,黄永胜,冯学智,等.SPOT卫星影像的水体提取方法及分类研究[J].遥感学报,2001,5(3):214-218.

Du J K,Huang Y S,Feng X Z,et al.Study on water bodies extraction and classification from SPOT image[J].Journal of Remote Sensing,2001,5(3):214-218.

[8] Zheng Q N,Laironald J,Huang N E,et al.Observation of ocean current response to 1998 hurricane georges in the Gulf of Mexico[J].Acta Oceanologica Sinica,2006,25(1):1-9.

[9] 杨玉静,冯建辉.纹理特征提取及辅助遥感影像分类技术研究[J].海洋测绘,2008,28(4):37-40.

Yang Y J,Feng J H.Research on extraction and assistant classification of remote sensing for texture feature[J].Hydrographic Survey and Charting,2008,28(4):37-40.

[10]Liu Z S,Wu S H,Li H,et al.Operational observations of three dimensional wind with incoherent Doppler wind LiDAR[C]//Proc of 25th International Laser Radar Conference,2010.

[11]张毅坤,麻晓畅,华灯鑫,等.基于EMD-DISPO的Mie散射激光雷达回波信号去噪方法研究[J].光谱学与光谱分析,2011,31(11):2996-3000.

Zhang Y K,Ma X C,Hua D X,et al.The Mie scattering LiDAR return signal denoising research based on EMD-DISPO[J].Spectroscopy and Spectral Analysis,2011,31(11):2996-3000.

[12]Yin S R,Wang W R.LiDAR signal denoising based on wavelet domain spatial filtering[C]//2006 CIE International Conference on Radar.Shanghai,China:Institute of Electrical and Electronics Engineers Inc,2007:1-3.

[13]李卿,张国平,刘洋.基于EMD的拉曼光谱去噪方法研究[J].光谱学与光谱分析,2009,29(1):142-144.

Li Q,Zhang G P,Liu Y.A study of Raman spectra denoising based on empirical mode decomposition[J].Spectroscopy and Spectral Analysis,2009,29(1):142-144.

[14]朱骥,林子瑜,王昂生.数字图像单个像元分形维数的特征与计算方法[J].光电工程,2005,32(2):24-27.

Zhu J,Lin Z Y,Wang A S.Features and algorithm of fractal dimension of single pixel in digital image[J].Opto-Electronic Engineering,2005,32(2):24-27.

[15]王娟,张军,吕兆峰.基于分形纹理的遥感影像土地覆盖的分类方法研究[J].测绘科学,2008,33(2):15-17,32.

Wang J,Zhang J,Lü Z F.Study on classification of land cover with remote sensing image based on fractal texture[J].Science of Surveying and Mapping,2008,33(2):15-17,32.