高光谱遥感图像频域自适应同态滤波薄云去除

维普资讯 http://www.cqvip.com 第３０卷第２期　物探化探计算技术　２００８年３月　文章编号：１ｏｏ１—１７４９（２００８）Ｏ２＿－０１６３—０６　高光谱遥感图像频域自适应　同态滤波薄云去除　刘　星，胡光道　（中国地质大学资源学院数地所，武－２７，．４３００７４）　摘要：针对高光谱遥感图像波段多，各波段受云层影响程度不一致，云层本身不均匀等特点，　提出了一种频域自适应同态滤波方法。该方法首先检测某个图像云层的范围和对应厚度，根据　云层处于低频特点，利用变差函数确定滤波窗口，把云层按照空间邻域大小转换为傅立叶空间，　根据云层厚薄调整高通滤波截止频率，然后反变换到空间域。该方法可解决频域空间同态滤波　方法破坏无云区的信息，克服空间域滤波方法效果不佳，难以控制和解释的缺陷。结果表明，该　方法在保持无云区信息不破坏前提下，较好地改善了不同厚薄云层的质量，计算机内存开销很　小，对大数据处理速度较快。　关键词：高光谱图像；频率域；自适应；云层去除　中图分类号：ＴＰ　７５　文献标识码：Ａ　同态滤波　Ｊ。基于图像本身，实用性和操作性较　０前言　强，通常在频率域进行，有些也拓展到空间域－８　Ｊ。　但很多文献　．８　指出，傅立叶变换消耗大量内存，　遥感图像获取通常受天气的干扰，尤其在我国　难以完成大数据量的图像处理。同时，频率域或小　西南地区，由于卫星通常在上午十点左右采集数　波域滤波去云以损害无云区信息为代价；而空间域　据，而数据受云层的影响较大，严重影响了遥感技　滤波通常需经验来调整参数　】，其结果很难解释，　术在该地区的解译和地质应用¨】。云层去除研究　难以控制去云结果。作者在本文中，从地质统计学　是遥感图像各种应用的前提，具有重要的现实意　空间相关原理，来决定滤波窗口的大小，在频域内　义　。］。目前，高光谱遥感以其波段多，波谱分辨　进行薄云去除，减少了内存开支，又避免了无云区　率高而广泛应用于地学领域。以往针对某个或少　信息的损失。　数几个波段的云层去除方法，难以适应大气窗口对　云层透射率不同特点　，而且，云层去除方法要么　１　同态滤波去云原理　用其它时相图像来替代局部图像－５】，或者用变换　的方法，例如Ｋ—Ｔ变换，去除喊噪声的分量，操作　１．１薄云区成像原理［　．１０川　性和效果不太现实，信息通常产生损失　Ｊ。对于　在有云区域，遥感卫星所接收到的图像信号，　小块薄云，利用地质统计学原理进行克立格插值，　是由太阳辐射经云层反射部分，太阳辐射经地物或　也能起到良好的效果，但对大块云区效果却不理　地面反射后，再穿透云层这二个部分组成，用公式　相［７】　ｔ：ｎ　ｏ　表达为　】　目前，对于大面积薄云的去除方法，通常采用　，Ｙ）＝　（　，Ｙ）・ｆｒ（　，Ｙ）　（１）　基金项目：中国地质大学优秀青年基金资助（ＣＵＧＱＮＬ０５１８）；地质过程与矿产资源国家重点实验室开放基金资助（ＭＧＭ　Ｒ一２００２一Ｏ８）　收稿日期：２００７—０４—１０　改回日期：２００７—１２—１８　维普资讯 http://www.cqvip.com

１６４　物探化探计算技术　３０卷　式中　（　，Ｙ）为照射分量，其频谱特性集中在低　频波段，是云层分布的函数；　（　，Ｙ）为反射分量，　是地物的函数，其频率处于高频区域。去云过程就　是去除ｆｉ（ｘ，Ｙ），求　（　，Ｙ）的过程，一般的滤波　器无法把ｆｉ（ｘ，Ｙ）去掉¨。。。　１．２同态滤波　云层的判别主要借鉴谢华美　的方法，首先　将高光谱反射率图像线性拉伸到０灰阶～２５５灰　阶（不拉伸需要调整参数），统计１７，×１７，邻域窗口内　的平均值和方差，并求出邻域窗口的最小平均值和　最大方差。对每个像素为中心的邻域，求出该邻域　平均值与最大平均值的差，以及方差与最小方差的　差值。因为有云分布的区域，平均值偏大，方差偏　小，把这二个值取绝对值相加，即对薄云分布打分。　求出分值最大和最小的值，利用下式确定云层厚　同态滤波把频率过滤作为图像处理的核心。　其原理是，对式（１）二边求对数，则有　ｌｎ／－（　，Ｙ）＝Ｉｎ　（　，Ｙ）＋ｌｎ　（　，Ｙ）　（２）　这样就把照射分量和反射分量分离，并转换到　度。　傅立叶空间。　（　）＝　×　Ｆ（Ｈ，　）＝ｆｌ　ｌｎｆ（ｘ，Ｙ）ｌ＝ｆｌｌｎｆｉ（ｘ，Ｙ）ｌ＋　／／，／，／，／／／，ｒＯＷ　（７）　厂Ｉｌ　（　，Ｙ）ｌ＝，（Ｈ，　）＋　式中ｓｃ口２ｅ（　，Ｙ）是各个像素对应的云层厚度；　Ｒ（Ｈ，　）　（３）　ｓｃ口　是按照邻域求出的各个像素的打分；／ｔ＂１．ｔ／ｎ一　然后，利用高通滤波，提取高频成份，过滤低　删是最大的截止频率；ｍｉｎ—ｓｃ口　和ｍａｘ—ｓｃ口ｒｅ分　频，即可以把反映地物细节的信息加强，从而降低　别是有云区像素的最大和最小得分，分值越小，云　薄云层照射噪声。　层越厚。　Ｇ（Ｈ，　）＝，（Ｈ，　）日（Ｈ，　）＋　２．２利用变差函数确定滤波窗口　Ｒ（Ｈ，　）日（Ｈ，　）　（４）　变差函数是地质统计学中研究分布于空间，并　再把傅立叶滤波结果反变换到空间域　呈现一定结构性和随机性的自然现象，具有揭示周　ｌｎｆ（　，Ｙ）＝　｛Ｇ（ｕ，　）｝＝　｛，（ｕ，　）　期性和无周期性生态参数本质的能力－ｌ引。几乎所　Ｈ（Ｈ，　）｝＋厂　｛Ｒ（Ｈ，　）日（Ｈ，　）｝　（５）　有的空间数据，都具有某种空间依赖，即空间自相　对滤波结果图像求自然指数　关或空间关联　一个区域单元上的某种地理现　Ｇ（　，Ｙ）＝ｅｘｐ［１ｎｆ（ｘ，Ｙ）］＝　象，或某一属性值与邻近区域单元上的同一现象或　属性值相关－ｌ引。遥感图像的反射率同样服从这个　￣ｘｐＥ２　｛，（Ｈ，　）日（Ｈ，　）｝＋　规律，即在一定邻域内的值相互依赖，范围之外则　７　｛Ｒ（ｕ，　）日（ｕ，　）｝］　（６）　没有空间相关性，所以对有云区的滤波可以限定在　这是同态滤波的基本原理，按照传统路线，对　变程之内，而不用考虑整个图像－１　，这样就可以节　整体图像滤波，存在以下问题¨　Ｊ引：　约计算机内存，有利于大数据量运算。　（１）傅立叶变换消耗大量的时间和空间，难以　邻域大、小的确定要逐步进行，通常截取１００　完成大数据量运算。　×１００的数据，就能达到确定变程的目的。如果变　（２）对有云区域进行处理的时候，同时处理了　程大于截取的数据，则要增加截取数据范围计算变　无云的区域，破坏了无云区域的图像。　程。本次研究采用的数据变程在各个方向和不同　（３）用同态滤波的方法去云，会产生严重的边　区域中，都处于１０个像素左右（如下页图１所示，　界效应，边界亮度增高，会影响图像判读。　横坐标为变程，纵坐标为变差函数，虚线为基台值，　点线为实验变差曲线、实线为拟合变差曲线），因　２云区判别和滤波窗口确定　此认为云层的高频和低频部份相关性也在１０×１０　的窗口之内，这个窗口也是上节云层判别滑动窗口　２．１云层判别　大小的依据。　由于传统同态滤波方法是针对整个图像进行　处理，这将对无云区产生破坏，也增大了计算机内　３薄云的去除　存开销。如果能够判别云层，又能确定影响云区的　空间结构范围，则能够有效保护无云区信息，也能　３．１滤波器设计　节约计算机的开销。　由于地表景物占据了邻域的高频成份，所以对　维普资讯 http://www.cqvip.com

２期　刘星等：高光谱遥感图像频域自适应同态滤波薄云去除　１６５　４０　０　３０　２０　ｌ０　０　０　３　６　９　１２　ｌｈｌ　ｌ５　ｌ８　２ｌ　０　３　６　９　１２　　ｌｈｌ　１５　１８　２１　（ａ）Ｄｉｒｅｃｔｉｏｎ　０　（ｂ）Ｄｉｒｅｃｔｉｏｎ４５　６４　５６　４８　三４２　０　２８　ｌ４　０　０　３　６　９　１２　ｌｈｌ　１５　１８　２ｌ　３２　ｌ６　０　０　３　６　９　１２　　ｌｈｌ　１５　１８　２Ｉ　（ｃ）Ｄｉｒｅｃｔｉｏｎ　９０　（ｄ）Ｄｉｒ￣ｆｉｏｎ１３５　图１　图像各个方向变程　Ｆｉｇ．１　Ｖａｒｉａｂｌｅ　ｄｉｓｔａｎｃｅ　ｏｆ　ｆｏｕｒ　ｄｉｒｅｃｔｉｏｎ　ｉｎ　ｔｈｅ　ｉｍａｇｅ　薄云覆盖的去除，主要是需过滤掉低频成份，保留　高频成份，通常采用巴特沃思高通滤波器　”　（取　１　ｌｎ　ｓｃ）＝一ｅｘｐ（１ｎ（ｓｃａｌｅ））　南１　＋　Ｉｋ　Ｘ　　＝Ｑ　Ｉ　（８）　式中Ｄ０为截止频率，以往Ｄ０为固定值，考虑到　通常，式（１２）中ｋ取０．１４１或者１，滤波阶次　为ｌ。该滤波器无振铃，图像模糊程度低。　根据上述原理，薄云去除的基本步骤为：　３．２　自适应同态滤波薄云去除　云的厚薄程度不同，产生的照射和反射分量不同，　（１）对经大气校正的反射率图像做灰度拉伸，　Ｄ０应随云层厚度自适应变化。由公式（２）可知　拉伸范围０—２５５。　ｆ　（　，Ｙ）一　（　（　，Ｙ））　（９）　（２）截取１／２或１／４块图像，转成ＡＳＳＣＩＩ。计　而　（　，Ｙ）则与云的厚度系数ｓｃａｌｅ相关，即，　算图像四个方向上的变差函数，即０。、４５。、９０。、　（　，），）＝￣（ｓｃａｌｅ（ｘ，），）），由于傅立叶变换是线性　１３５。四个方向，取最大变程作为滑动窗口邻域大　的，由公式（４）进一步可知　小。如果变程较大，则考虑截取更大区域图像。　，（　，　）一　（Ｉｎ（ｓｃａｌｅ））　（１Ｏ）　（３）识别薄云并判定云层的厚度。通常云区　因此，公式（８）中的截止频率Ｄ０应该根据　均值取大于６Ｏ，得分小于０．７的邻域。　ｓｃａｌｅ的变化进行调整，才能有效抑止，（　，　）分　量。增强　（　，　）分量。既然，（　，ｔ，）一　（Ｉｎ（ｓｃａｌｅ）），那么　日（　，　）　（　１　ｌｎ　（４）根据邻域大小，截取云区图像，取对数值，　进行傅立叶变换。利用自适应滤波器对该邻域的　傅立叶变换进行高通滤波，再进行反变换到空间　域，并做自然指数变换。　）　）：　‘，　１　（１１）　（５）由于傅立叶变换会产生边界效应，不同邻　日（　，ｔ，）×　１　、　域间数值差异较大，各个邻域的均值比在变换前偏　大，因此要对各邻域滤波结果进行调整。通常采用　以下公式　ｆｔｗ：ｆｗ　ｉ—ｔ—ｔｅｄ—＿＿ｗｉ—ｎ＿－＿『＿ｍｉｎｔ＋ｍｉｍｌｎｏ＝—————　————■　一十　一　ｓｃａｒ，ｅ—Ｄ（　，　）　（１２）　。　ｍａｘｔ—ｍｌｎｔ　（１３）１ｊ　　式中　是变换后的像素亮度；ｓｃａｌｅ是云厚度；　ｉｔｌｔｅｄ—ｗｉｎ是经傅立叶滤波后的像素值；ｍａｘｏ、　维普资讯 http://www.cqvip.com

ｌ６６　物探化探计算技术　３０卷　ｍｉｎｏ、ｍａｘｌ、ｍｉｎ１分别表示原始图像邻域的最大最　信息提高许多，并且未破坏无云区。　３．４高光谱图像光谱保真　高光谱图像处理的一个基本要求，就是谱线高　保真度。由于高光谱各个波段对地物的反射率是　不一致的，各个波段受云层影响程度也不一致，处　理后会产生部分光谱变形。为了保证光谱的保真　小值、滤波变换后邻域的最大、最小值。经过处理，　边界效应减小，云区的对比度和未处理区保持一致　（见图２）。　（６）由于高光谱图像波段多，各个波段受云层　的影响差别不一致，需逐个处理波段，但无须调整　参数。根据算法编制的程序，能够逐个处理多波段　图像。　３．３　自适应同态滤波薄云去除与传统同态滤波　性，需要把滤波后图像重新拉伸为原来的数值范围　（见式（１３）），同时降低平均亮度。为了检验本方　法处理结果是否保持有光谱真实性，本研究从１８３　个波段中，按等距抽取了１０个波段，在有云区随机　选择了三个像素点，形成波普曲线，与未处理之前　光谱作相关分析。从对比结果图３可以看出，在检　验水平为０．Ｏ１时，本方法在去云后光谱特征与处　理前相关程度较高。　方法相比较　按照传统的方法，含有Ⅳ个像素的图像进行　傅里叶变换，其时问复杂度为０（　），另外还需要　１０　Ｎ的辅助空间　ｊ。该方法消耗大量的时问和空　间，对于大数据，一幅１００　Ｍ字节的存储图像经变　换后的结果将达到８００　Ｍ　，照此，目前的计算机　难以支持大于１Ｇ的数据处理。而本文的方法内　存开销复杂度为０（ｙ　）（ｙ为变程），时间开销为　０（ｙ　）×ｎｕｎ—ｗｉｎ（ｎｕｎ—ｗｉｎ为滑动窗口个数），基　４结论　作者在本文中，从薄云区成像原理出发，推导　本上不存在大数据无法计算问题。　本方法在处理结果上与传统的方法对比，细节　出一种随云层变化的自适应频率域滤波方法，对薄　云的影响进行去除。利用变差函数，确定了图像的　（ａ）原始２８波段图像　（ｂ）２８波段应用　（ｃ）１００波段图像　（ｄ】１００波段应用　本文方法去云的结果　本文方法去云的结果　（ｃ）２８波段应用传统同态滤　波方法（截止频率１　５０）得　到的结果　（Ｄ　２８波段利用传统方法　（截止频率１Ｏ）的结果　（ｇ）１００波段应用传统方法　取截止频率为１５Ｏ的结　果［ｍ１】　‘ｈ）１Ｏ０波段应用传统方法　取截止频率为１Ｏ的结　果［ｍ１】　图２　自适应同态滤波结果图像与原始图像对比　Ｆｉｇ．２　Ｃｏｎｔｒａｓｔｓ　ｂｅｔｗｅｅｎ　ｔｗｏ　ｏｒｉ　ｎａｌ＿ｂａｎｄｓ　ａｎｄ　ｔｈｅ　ｒｅｓｕｌｔｓ　ｐｒｏｃｅｓｓｅｄ　ｂｙ　ｔｈｅ　ｓｅｌｆ－ａｄａｐｔｉｖｅ　ｈｏｍｏｍｏ￣ｅ　ｆｉｌｔｅｒ　维普资讯 http://www.cqvip.com ２期　刘星等：高光谱遥感图像频域自适应同态滤波薄云去除　１６７　譬　兰　享　享　４　６　８　１０　２　６　８　１０　ＢａｎｄＮｕｍｂｅｒ　ＢａｎｄＮｕｍｂｅｒ　（ａ）处理前像素５０—２８０处光谱曲线，均值０．５９７　（ｂ）经处理后均值０．６１，相关系数０．９０６　８　墨　４　６　８　１０　２　４　６　８　１０　ＢａｎｄＮｕｍｂｅｒ　ＢａｎｄＮｕｍｂｅｒ　（ｃ）处理前像素９０—１ｌＯ处光谱曲线，均值０．４１０　２　（ｄ）处理后均值０．４０９　２，相关系数０．９２０　譬　ｉ　呈　２　４　６　８　１０　２　４　６　８　１０　ＢａｎｄＮｕｍｂｅｒ　ＢａｎｄＮｕｍｂｅｒ　（ｅ）像素１００．２００Ｊ［１２光谱曲线，均值０．５４５　８　（０处理后均值０．５４４　５，相关系数０．９６０　５　图３　薄云去除前、后光谱对比　Ｆｉｇ．３　Ａ　ｃｏｎｔｒａｓｔ　ｓｋｅｔｃｈ　ｂｅｔｗｅｅｎ　ｓｐｅｃｔｒｕｍ　ｏｆ　ａｌｌ　ｏｒｉｇｉｎａｌ　ａｎｄ　ａ　ｃｌｏｕｄｓ　ｌ＇￣ｍｏｖａｌ　ｉｍａｇｅ　空间自相关区域，并在该邻域内进行薄云去除，减　［２］宋小宁，赵英时．ＭＯＤＩＳ图像的云检测及分析［Ｊ］．　少了计算机内存开销，使得以往不能完成的大数据　中国图象图形学报，２００３，８Ａ（９）：１０７９．　量处理变成现实。同时，本算法还能够自动检测云　［４］杨凯，卢健．遥感图象处理原理和方法［Ｍ］．北京：测　绘出版社，１９８８．　层分布及其厚度，根据云层影响程度，调整滤波截　［５］　王惠，谭兵，沈志云．多源遥感影像的去云层处理　止频率，能够有效地处理云层变化较大的薄云区。　［Ｊ］．测绘学院学报，２００１，１８（３）：１９５．　运算速度很快，薄云去除效果有所改进。无须调整　［６］王润，刘洪斌，宫瑞．多光谱遥感图像去云方法［Ｊ］．　参数，能够适合同一地区高光谱不同波段数据的处　计算机与现代化，２００５，（６）：１３．　理，并且减小了傅立叶变换后的边界效应。　［７］　陈姚，王金亮，李石华．遥感图像中云层遮挡影响消　本文算法基于ＩＤＬ语言编写，实验数据来自云　除方法研究述评［Ｊ］．国土资源遥感，２００３，１（６７）：　南普郎地区的美国ＥＯ一１　Ｈｙｐｅｒｉｏｎ　Ｌｅｖｅｌｌ高光谱　６１．　影像，数据经过大气校正，并转换为反射率。　［８］　谢华美，李荣艳，田艳琴，等．基于大数据量遥感图像　的薄云去除［Ｊ］，北京师范大学学报（自然科学版），　参考文献：　２００６，４２（１）：４２．　［９］武平．用ＥＮＶＩ软件去除薄云的方法探讨［Ｊ］．有　［１］冯春，马建文，戴芹，等．一种改进的遥感图像薄云快　色金属，２００５，２９．　速去除方法［Ｊ３．国土资源遥感，２００４，４（６２）：１．　［１０］　赵忠明，朱重光．遥感图象中薄云的去除方法［Ｊ］．　维普资讯 http://www.cqvip.com

１６８　环境遥感，１９９６，１１（３）：１９５．　物探化探计算技术　３０卷　［１５］　吴刚，陈慧琼，杨敬安．基于图像变差函数的多尺度　边缘检测［Ｊ］．合肥工业大学学报（自然科学版），　２０００，２３（２）：２２７．　［１１］吴为禄．遥感图像中的云层消除处理［Ｊ］．铁路航　测，２００３，（１）：６．　［１３］谢华美，何启翱，郑宁，等．基于ＥＲＤＡＳ二次开发的　遥感图像同态滤波薄云去除算法的改进［Ｊ］．北京　师范大学学报（自然科学版），２００５，４１（２）：１５０．　［１４］侯景儒，尹镇南，李维明，等．实用地质统计学［Ｍ］．　北京：地质出版社，１９９８．　［１６］丰茂森．遥感图像数字处理［Ｍ］．北京：地质出版　社，１９９２．　作者简介：刘星（１９７４一），男，讲师，博士，从事成　矿预测研究。　征稿启事　《物探化探计算技术》是由四川省教育厅主管、成都理工大学与中国地质科学院物化探研究所共同主　办的科学技术刊物（双月刊），国内、外公开发行，其办刊宗旨是坚持理论与实践相结合的方针，分析研究　有关地球探测信息处理和计算机技术及其应用与发展理论。　１．征稿范围：　（１）地球物理、地球化学、遥感地质、数学地质、矿床地质、能源地质、水文地质、工程地质、环境地质、　灾害地质等探测与信息处理技术。　（２）国内、外地球探测信息处理和计算机技术及其应用的综合评价、发展动态和软件的开发、应用介　绍等。　２．来稿要求：　（１）文稿论点明确，论据充分、文字精炼、数据准确、文稿字数６　０００字左右为宜（含图表），并通过网　络将文章的电子邮件发到编辑部邮箱内。　（２）论文内容应包括：中、英文论文名称、摘要（含“研究的问题、过程和方法，以及应用的结果”等内　容，约５００字左右）；中、英文关键词；中、英文的图题、表题；对文中的同一项目术语，前后应统一；计算单　位一律使用中华人民共和国法定计量单位，并用外文符号表示，如ｋｍ、ｒｆｌ、Ｐａ。　（３）图、表应明确、规范、清晰、简洁，并请从网上提供图件。　（４）凡引用他人公开发表的成果，必须列入参考文献中。在正文引用的部位顺序标注，在正文顺序列　出参考文献目录，编排格式如下：　期刊：序号著者．文章题目．刊物名称，出版年份，卷（期），页次。　专著：序号著者．书名．出版地：出版者，出版年，页次。　凡引用未公开发表的内部资料，不得列入参考文献。文献中作者在三人以内，应全部列出；若超过三　人，请在写出前三人姓名之后再用“等”字表示。中、外作者姓名均需按“姓在前，名在后”的方式排列。　（５）请提供第一作者简介，其中包括第一作者的姓名、性别、出生年份、职称、学历、从事工作、重要著　作。通讯地址，邮编，联系电话，电子邮件。　（６）若论文有项目资助，应注明科研项目名称，以及项目编号。　（７）作者请自留底稿，稿件一经审查通过，本刊将函告作者。来稿采用与否，将在收稿三个月内回告　作者。稿件刊用后即付稿酬，并赠寄当期刊物二册。　来稿地址：四川省、成都市、成都理工大学（内）《物探化探计算技术》编辑部　邮政编码：６１００５９　电话号码：（０２８）８４０７８９９５　电子邮箱：ｗｔｈｔ＠ｃｄｕｔ．ｅｄｕ．ｃｎ　网　址：ｈｔｔｐ：／／ｗｔｈｔｊｓｊｓ．ｐｅｒｉｏｄｉｃａｌｓ．ｎｅｔ．ｃｎ　ｈｔｔｐ：／／ｗｔｈｔ．ｃｈｉｎａｊｏｕｒｎａ１．ｎｅｔ．ｃｎ