CN107862667A - 一种基于高分辨率遥感影像的城市阴影检测与去除方法 - Google Patents

一种基于高分辨率遥感影像的城市阴影检测与去除方法 Download PDF

Info

Publication number
CN107862667A
CN107862667A CN201711185046.5A CN201711185046A CN107862667A CN 107862667 A CN107862667 A CN 107862667A CN 201711185046 A CN201711185046 A CN 201711185046A CN 107862667 A CN107862667 A CN 107862667A
Authority
CN
China
Prior art keywords
shadow
mrow
region
shade
block
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201711185046.5A
Other languages
English (en)
Other versions
CN107862667B (zh
Inventor
张永军
岳照溪
毛炜青
郭功举
潘琛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
SHANGHAI INSTITUTE OF SURVEYING AND MAPPING
Wuhan University WHU
Original Assignee
SHANGHAI INSTITUTE OF SURVEYING AND MAPPING
Wuhan University WHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by SHANGHAI INSTITUTE OF SURVEYING AND MAPPING, Wuhan University WHU filed Critical SHANGHAI INSTITUTE OF SURVEYING AND MAPPING
Priority to CN201711185046.5A priority Critical patent/CN107862667B/zh
Publication of CN107862667A publication Critical patent/CN107862667A/zh
Application granted granted Critical
Publication of CN107862667B publication Critical patent/CN107862667B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/70Denoising; Smoothing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/13Edge detection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/136Segmentation; Edge detection involving thresholding
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/194Segmentation; Edge detection involving foreground-background segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/40Analysis of texture
    • G06T7/41Analysis of texture based on statistical description of texture
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10032Satellite or aerial image; Remote sensing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20024Filtering details
    • G06T2207/20028Bilateral filtering
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30181Earth observation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Probability & Statistics with Applications (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于高分辨率遥感影像的城市阴影检测去除方法,首先通过阈值对双边滤波后的影像计算得到特征值进行阴影检测,并利用图割法对该结果进行区域增长得到阴影掩码;同时用双边滤波前后的亮度图相减得到细节图,进而得到纹理较弱的区域,结合NDWI完成城市水域的检测,在阴影掩膜中去掉水域的部分,并分别对阴影与非阴影进行分割,计算影像光谱与纹理信息,匹配得到每块阴影区对应的同质非影区;最后用矩匹配和直方图匹配的方法分别去除非水域与水域的阴影。本发明针对城市高分影像地表的复杂性,大大改善了阴影检测与去除的效果,可用于改进城市变化检测与地物识别等工程应用问题。

Description

一种基于高分辨率遥感影像的城市阴影检测与去除方法
技术领域
本本发明属于遥感和摄影测量技术领域,涉及一种基于高分辨率遥感影像水域检测与城市阴影的检测去除方法,尤其是涉及一种基于高分辨率遥感影像优化矩匹配的阴影检测与去除方法。
背景技术
阴影的检测与去除在高分辨率遥感影像分析与应用中是一个关键问题。随着遥感技术的高速发展,影像空间分辨率已经达到较高的水平,城市遥感影像分析中,阴影提取可以辅助建筑物等高程信息的提取,阴影去除可以丰富影像信息,改善地物识别、变化检测以及影像匹配等问题。
高分辨率城市遥感影像包含地物信息丰富,城市水域复杂、阴影区域零散且分布无规则,增加了阴影检测与去除的难度。城市阴影与水域的光谱特征非常接近,需要将水域与地面的阴影进行区分处理才能得到更好的效果。水域提取一般用归一化水指数NDWI进行判别,但是城市水质复杂,单一的光谱特征不能很好的提取水域。阴影检测主要有色彩空间特征提取并结合阈值的方法,不需要过多的场景参数描述,简单可行,但是效果并不一定理想,也有结合区域增长与形态学处理的阴影检测方法的研究,但是因为在增长和形态学处理的过程中往往不能保证阴影边缘的准确性。阴影去除主要有线性拉伸、直方图匹配以及同态滤波等,但是存在纹理信息恢复不完整以及光谱特征保真度不高等问题,总之,传统的阴影检测与去除方法往往存在信息丢失,在高分辨率遥感影像中难以达到理想的效果,不能很好的用于后续的变化检测与地物识别。
近年来基于机器学习与图割法的检测法与矩匹配的阴影补偿算法已有一些研究和成果,如结合支持向量机与Grabcut的阴影检测方法,还有在影像分割基础上进行区域的矩匹配,从而去除阴影方法等,在自然影像的阴影去除方面都有了一定的成果,但是在高分辨率遥感影像中,由于地物类型复杂性的提升,自然影像的阴影去除方法在应用中,会出现分割结果不理想与分割块间色彩过渡的不平滑等问题,不能达到理想的效果。
发明内容
本发明主要是解决现有技术所存在高分辨率城市遥感影像地物复杂度高,以及分割结果不理想、分割块间色彩过渡的不平滑以及光谱特征保真度不高等问题;提供了一种结合水域检测与优化矩匹配的阴影检测与去除方法,可以有效改善城市高分辨率遥感影像阴影检测与去除的效果。
本发明所采用的技术方案是,一种基于高分辨率遥感影像城市阴影检测与去除方法,包括以下步骤:
步骤1,对原始高分辨率影像进行初步阴影检测,并用图割法进行区域增长,得到初始的阴影区和非阴影区;
步骤2,对原始高分辨率影像的亮度分量I进行直方图均衡化得到图IHE,并对IHE进行双边滤波;
步骤3,用步骤2中的IHE亮度图与步骤2中的双边滤波结果进行差值处理,得到细节纹理图,利用Otsu阈值法在细节图上得到纹理较弱的区域,并结合NDWI得到水域;
步骤4,步骤1中的初始阴影区域去除水域的部分,得到城市非水域阴影掩膜,并确定阴影边缘;
步骤5,计算原始影像的光谱特征值与纹理特征值,并将阴影与非阴影区分别进行SLIC影像分割,得到阴影与非阴影区分割块;
步骤6,根据步骤5中计算得到的光谱特征值与纹理特征值,对每个阴影区内的分割块,在周围临近的非阴影区中计算出特征最相近的分割块,然后分别两个对应分割块的各波段的均值和方差;
步骤7,对步骤6的计算结果,用距离作为权值整合每个阴影区与周围非阴影区所有对应分割块的波段特征值,计算出该阴影区的整体校正参数,并对阴影区进行分波段逐像素的校正;
步骤8,校正好阴影区域之后,对于步骤4中的阴影边缘,计算该阴影边缘两侧阴影分割块与对应非阴影分割块平均的均值与方差,并对阴影边缘进行分波段逐像素校正;
步骤9,对于步骤3中检测得的水域部分,计算亮度直方图,获得水域的阴影部分,并对水域中阴影与非阴影部分的直方图进行匹配,得到无阴影的水域。
进一步的,步骤7的具体实现方式如下,
步骤7.1,首先计算整个阴影区域内的均值与方差,然后用距离作为权值,整合非阴影对应区域的均值与方差,如下式所示,
其中,C(i,j)表示该阴影区所对应的非阴影区整体的加权均值和方差,为非阴影区每个分割块的均值和方差组成的向量,CB(m)是每个分割块的校正参数,m是阴影区内对应的分割块数,为非阴影区不同分割块与对应阴影分割块的距离组成的权向量;
步骤7.2,按照上一步的结果,用整体阴影和非阴影区的均值和方差进行计算,按照下式进行分波段逐像素的校正,
其中,Is′(i,j)与Is(i,j)分别为阴影区校正前后各像素的值,μn与μs分别为非阴影区与阴影区整体的均值,σn-s为阴影区与非阴影区的协方差,与σs为阴影区的方差。
进一步的,步骤8的具体实现方式如下,
首先在与该阴影边缘临近的同质分割块中寻找阴影和非阴影两个分割块作为参考,然后计算两参考块的均值与方差的平均值,并用按照下式进行校正,
其中Ie′(i,j)与Ie(i,j)分别为阴影边缘校正前后各像素的值,μavg和σavg分别为阴影边缘两侧的阴影分割块与对应非阴影分割块平均的均值和方差,μe和σe分别为阴影边缘的均值和方差。
进一步的,步骤1的具体实现方式如下,
步骤1.1,对影像进行分波段的双边滤波,并利用四波段数据计算得到的B’、I、Q、A四种特征值,分别用Otsu阈值法进行处理,得到的四个结果的交集,即为初步阴影检测结果,其中B′=B/(R+G+B),R、G、B分别表示红、绿、蓝三波段数据,B′为蓝波段的归一化值,是阴影的主特征值,I为HSI色彩空间下的亮度值,Q为复合特征Q=B′-I,A为复合特征,
其中NDVI为归一化植被指数,NDVI=(NR-R)/(NR+R),t为判定植被的阈值;
步骤1.2,用图割法对阴影进行区域增长,认为待增长区是前景,非阴影区是可能背景,图割后得到的前景与可能前景都视为阴影区,背景与可能背景都视为非阴影区。
进一步的,步骤2中对IHE进行双边滤波的实现方式如下,
用双边滤波对IHE图进行处理得到滤波影像,双边滤波器中:
式中IBF是双边滤波后得到的影像,i,j分别为每个像素的横纵坐标,权重系数ω(i,j)=Ws(i,j)×Wr(i,j),其中Ws(i,j)和Wr(i,j)分别为IHE图的相似性因子和空间距离判定因子,σs和σr表示高斯函数标准偏差,(i,j)表示待滤波中心像素坐标,(i’,j’)表示窗口其他像素坐标,IHE(i’,j’)表示双边滤波后中心像素的亮度值。
进一步的,步骤4中通过对城市非水域阴影掩膜影像各进行一次腐蚀和一次膨胀,然后将两次的结果相减,得到的结果视为阴影边缘。
进一步的,步骤5的具体实现方式如下,
步骤5.1,用Hog算子提取影像的弱纹理强度值,并用Gabor变化提取不同尺度和方向的纹理向量式中σklkl分别是每个分割区块的小波变换系数的方差与均值,k、l分别表示小波变换的尺度和方向;
步骤5.2,将原始影像从RGB颜色空间转换到CIE-Lab颜色空间,对应每个像素的(L,a,b)颜色值和(x,y)坐标组成一个5维向量V[L,a,b,x,y];
步骤5.3,按照城市非水域阴影掩膜将影像分为阴影和非阴影两张影像,并分成K个种子点,根据步骤5.2中计算所得像素间V向量的距离大小,对每个种子点的周围空间进行K-means聚类,得到最终SLIC分割结果。
进一步的,步骤6中对每个阴影区内的分割块,在周围临近的非阴影区中计算出特征最相近的分割块,实现方式如下,
针对每个阴影区内的分割块,计算其周围非阴影区各分割块的Gabor向量间的距离dG,结合L光谱分量、Hog算子,计算相似性,并用距离作为权值,找到与该阴影分割块地物类型最接近的非阴影区分割块,
式中Dist(s,n)表示非阴影区与阴影区内分割块特征值的欧氏距离,距离最小的表示非阴影区地物类型最接近阴影区的分割块,s表示阴影区分割块,n表示非阴影区分割块,为表示两个分割块距离的权值,(Xs,Ys)、(Xn,Yn)分别为阴影分割块与非阴影分割块的中心像素坐标。
本发明具有如下优点:城市水域检测中综合考虑纹理和光谱特征,改善了城市水质复杂难以检测的问题,同时将地面阴影与水面阴影分别去除,改善了阴影去除的结果;阴影与非阴影同质区域匹配的时候,同时考虑了大范围特征Hog算子和细节特征如Gabor纹理特征和色彩特征,增强了同质区的匹配效果;在阴影去除的过程中,改善了以分割块为单位的亮度拉伸,而是以阴影块为整体计算校正参数,增强了结果图色彩的平滑性。
附图说明
图1为本发明实施例的总体流程图;
图2为本发明实施例的水域检测结果示意图。
图3为本发明实施例的阴影检测结果示意图。
图4为本发明实施例的阴影区与非阴影区分别分割示意图。
图5为本发明实施例的非影区对应的拉伸参数的整合示意图。
图6为本发明实施例的阴影去除结果示意图。
具体实施方式
为了便于本领域普通技术人员理解和实施本发明,下面结合附图及实施例对本发明作进一步的详细描述,应当理解,此处所描述的实施示例仅用于说明和解释本发明,并不用于限定本发明。
请见图1,本发明提供的一种基于高分辨率遥感影像城市阴影检测与去除方法,包括以下步骤:
步骤1:对高分辨率影像四波段数据(R、G、B、NR)进行双边滤波,并用四波段数据计算B’、I、Q、A四种特征值,用Otsu阈值法[1]得进行初步阴影检测,并用图割法[2]进行区域增长,
[1]高贤君,万幼川,杨元维,等.高分辨率遥感影像阴影的自动检测与自动补偿[J].2014.
[2]Carsten Rother,Vladimir Kolmogorov,Andrew Blake,“GrabCut”-Interactive Foreground Extraction using Iterated Graph Cuts.[J]2004.
具体实现包括以下子步骤:
步骤1.1:首先对影像进行分波段的双边滤波(对每一波段进行双边滤波时参数一样),去除不必要的细节信息,然后利用光谱信息进行阴影检测(计算得到的B’、I、Q、A四种特征值,需要分别用Otsu阈值法进行处理,得到的四个结果的交集,就是阴影检测结果),其中B′=B/(R+G+B)(等式右边RGB分别表示红绿蓝三波段数据),即为蓝波段的归一化值,是阴影的主特征值,I为HSI色彩空间下的亮度值,Q为复合特征Q=B′-I,A为复合特征:
其中NDVI为归一化植被指数,这个指数是已经广泛认可的植被检测指数,计算方法是(NR-R)/(NR+R),t为判定植被的阈值,一般设为0.4。
步骤1.2:用图割法对阴影进行区域增长,认为待增长区(即阴影区)是前景,非阴影区是可能背景,图割后得到的前景与可能前景都视为阴影区,背景与可能背景都视为非阴影区。
步骤2:对影像的亮度分量I进行直方图均衡化得到图IHE,并对IHE进行双边滤波;
用双边滤波对IHE图进行处理得到滤波影像,双边滤波器中:
式中IBF是双边滤波后得到的影像,i,j分别为每个像素的横纵坐标,权重系数ω(i,j)为两部分的乘积,在这里定义ω(i,j)=Ws(i,j)×Wr(i,j),其中Ws(i,j)和Wr(i,j)分别为IHE图的相似性因子和空间距离判定因子,σs和σr表示高斯函数标准偏差,本发明实施例中σs取值2,σr取值为IHE(MAX)/10,IHE(MAX)表示均衡化亮度最大值,e表示数学常量,(i,j)表示待滤波中心像素坐标,(i’,j’)表示窗口其他像素坐标,IHE(i’,j’)表示双边滤波后中心像素的亮度值。
步骤3:用步骤2中的IHE亮度图与步骤2中的双边滤波结果进行差值处理,得到细节纹理图,用Otsu阈值法得到纹理较弱的区域,并结合NDWI得到水域;
即先将IHE与双边滤波结果相减得到细节图,并用Otsu提取弱纹理区(是在细节图上用阈值提取弱纹理区,计算方法就是直接对细节图做Otsu阈值处理,将细节图分为“强纹理区”和“弱纹理区”,因为细节图的结果是灰度数据,所以通过选择合适的阈值就能得到弱纹理区,得到的结果为二值图,即强纹理区取值为0,弱纹理区取值为1),再利用四波段信息计算NDWI(这个指数是已经广泛认可的植被检测指数,计算方法是(NR-R)/(NR+R)),用阈值法去除纹理同样较弱的道路、平台及操场等人工地物,弱纹理区是包括水域、道路、平台等,但是这些地物的NDWI,即归一化水指数很小,在弱纹理区检测结果的基础上,将NDWI值小的像素值改为0,最终得到的结果图中,取值为0的就是非水域,取值为1的就是水域,此结果作为水域检测结果。
请见图2水域检测结果示意图,a图为有水域的城市遥感影像,b图为人工绘制的水域掩膜,c图为只用NDWI检测的水域结果,可以看出水域和阴影混在一起无法区分,d图为IHE亮度图与步骤2中的滤波结果进行差值处理,得到细节纹理图,可以看出水域和道路的部分纹理很弱,e图是细节信息和NDWI结合得到的水域检测结果,f图为去除小面积面元后的阴影检测结果。
步骤4;步骤1中的阴影区域去除水域的部分(水域就是步骤3的最后结果),得到城市非水域阴影掩膜,确定阴影边缘,即对影像(这个影像就是城市非水域阴影掩膜影像,膨胀腐蚀都是对这个结果进行的)各进行一次腐蚀和一次膨胀,然后将两次的结果相减,得到的结果视为阴影边缘,然后针对城市非水域阴影掩膜进行小面积面元(本发明实施例中,小面积面元为面积小于50像素的空洞或者碎块)的去除,防止像素级阴影检测导致的细小空洞和过于细碎的部分,减小误检和错检的现象,得到更为完整且整洁的阴影掩膜。
请见图3阴影检测结果示意图,a图是城市遥感影像,b图是利用传统特征法,提取影像特征并结合阈值得到的阴影检测结果,可以看出非常破碎,误检和漏检的情况都有,c图是本发明阴影检测的初步结果,可以看出已经有了较好的改善,d图是对在c图结果的基础上,去除小面积,得到最终的阴影检测结果。
步骤5:计算影像的光谱特征值与纹理特征值,并将阴影与非阴影区分别进行SLIC影像分割,请见图4阴影区与非阴影区分别分割示意图:
步骤5.1:用Hog算子提取影像的弱纹理强度值,并用Gabor[3]变化提取不同尺度和方向的纹理向量其中k、l取值都为0~4,式中σklkl分别是每个分割区块的小波变换系数的方差与均值,k、l分别表示小波变换的尺度和方向。此步得到的纹理和特征信息会用于步骤6中的同质分割块计算与匹配。
步骤5.2:将原始影像从RGB颜色空间转换到CIE-Lab颜色空间,对应每个像素的(L,a,b)颜色值和(x,y)坐标组成一个5维向量V[L,a,b,x,y];
步骤5.3:按照阴影掩膜将影像分为阴影和非阴影两张影像,并分成K个种子点,对每个种子点的周围空间进行K-means聚类,聚类的依据就是步骤5.2中计算所得像素间V向量的距离大小,然后多次迭代得到最终SLIC分割结果;因为高分辨率影像里面地物类型比较多,如果对整块阴影区进行亮度的拉伸,效果不佳,为了避免一整块都变亮显得很奇怪,本发明按照地物类型进行分块,进行进一步分析。
步骤6:对每个阴影区内的分割块(即SLIC分割结果),在周围临近的非阴影区中计算出特征最相近的分割块(即同质分割块),然后分别计算对应分割块的各波段的均值和方差;
针对每个阴影分割块,计算其周围非阴影区各分割块的Gabor向量间的距离dG(城市建筑物阴影一般是块状的,每块阴影区周围有非阴影区环绕,对阴影与非阴影区都进行分割得到小的分割块,针对每一个阴影区内的小分块,遍历周围非阴影区的小分块,就可以在非阴影区找到与阴影区内小分块地物类型最相近的分割块,地物类型最接近的判别依据就是式叁,进行计算得到的Dist(s,n)最小的非阴影与阴影分块,就是地物类型最接近的),并且结合L光谱分量(L光谱分量就是步骤5.2中五维向量V[L,a,b,x,y]里面的L)、Hog算子[4],计算相似性,并用距离作为权值,找到与该阴影块地物类型最接近的非阴影区分割块:
式中Dist(s,n)表示非阴影区与阴影区内分割块特征值的欧氏距离,距离最小的表示非阴影区地物类型最接近阴影区的分割块,s表示阴影区分割块,n表示非阴影区分割块,为表示两个分割块距离的权值,(Xs,Ys)、(Xn,Yn)分别为阴影区分割块与非阴影区分割块的中心像素坐标。
[3]朱健翔,苏光大,李迎春.结合Gabor特征与Adaboost的人脸表情识别[J].光电子·激光,2016.
[4]郭金鑫,陈玮.基于HOG多特征融合与随机森林的人脸识别[J].计算机科学,2014.
步骤7:对步骤6的计算结果,用距离作为权值整合每个阴影区与周围非阴影区所有对应分割块的波段特征值,计算出该阴影区的整体校正参数(即步骤7.1的C(i,j)),并进行分波段(R/G/B/NR四波段)逐像素的校正,具体实现包括以下子步骤:
图5为非影区对应的拉伸参数的整合示意图,针对阴影区内的每一个分割块,在周围非阴影区寻找纹理信息最相近的对应非阴影分割块,图中的箭头表示计算所得的阴影区与非阴影区对应的分割块,然后分别计算两分割块的均值与方差,并计算两个分割块的中心像素的距离以用作步骤7中的权值;整个阴影区的所有对应均值与方差记录到步骤7中的中,用以进行整块阴影区阴影去除。
步骤7.1:首先计算整个阴影区域内的均值与方差,然后用距离作为权值,整合非阴影对应区域的均值与方差,如下式所示:
其中C(i,j)表示该阴影区所对应的非阴影区整体的加权均值和方差,为非阴影区每个分割块的均值和方差组成的向量(是包含均值和方差的一个二维向量,因为在阴影去除这一步,需要对阴影区域进行亮度的拉伸,落实到计算上,就是需要对每一波段的均值与方差进行调整,使得阴影区的均值与方差,与同质的非阴影区(就是在非阴影区内,与阴影区地物类型相同或相似的区域)均值与方差相同或者接近,从而达到去除阴影的效果),CB(m)是每个分割块的校正参数,包括该块的均值与方差,m是阴影区内对应的分割块数,为非阴影区不同分割块与对应阴影区分割块的距离(这个距离就是对应的阴影分割块与非阴影分割块中心像素的欧氏距离,分割块的中心像素计算方式就是分割块内的所有像素横纵坐标的均值)组成的权向量;
步骤7.2:按照上一步的结果,用整体阴影和非阴影区的均值和方差进行计算,按照下式进行分波段逐像素的校正:
其中Is′(i,j)与Is(i,j)分别为阴影区校正前后各像素的值,μn与μs分别为非阴影区与阴影区整体的均值(是整体的均值),σn-s为阴影区与非阴影区的协方差,与σs为阴影区的方差,根据式肆中的C(i,j)获得。
步骤8:校正好阴影区域之后,对于步骤4中的阴影边缘,计算该的阴影边缘两侧阴影分割块与对应非阴影分割块平均的均值与方差,并将阴影边缘逐像素拉伸为阴影与非阴影分割块的平均值;即首先在与该区域临近的同质分割块中寻找阴影和非阴影两个分割块作为参考,然后计算两参考块的均值与方差的平均值,并用按照下式分波段进行校正:
其中Ie′(i,j)与Ie(i,j)分别为阴影边缘校正前后各像素的值,Ie(i,j)是经过步骤7.2处理之后的值,μavg和σavg分别为阴影边缘两侧的阴影分割块与对应非阴影分割块平均的均值和方差,μe和σe分别为阴影边缘的均值和方差。
步骤9:对于步骤3中检测所得的水域部分,计算亮度值,并对水域中的阴影部分与非阴影部分
进行直方图匹配。计算原影像水域部分的亮度直方图,然后用Otsu阈值法得到水域的阴影部分,分别计算水域中阴影与非阴影的直方图,将阴影区直方图与非阴影区直方图进行匹配[5],得到无阴影的水域。其中,直方图匹配或叫做直方图规定化,是把原图像的直方图按照给定的直方图加以映射,使得新的图像直方图分布类似于给定的函数,这里把水域非阴影区的作为给定直方图,水域非阴影区作为原影像直方图然后进行处理。处理的具体步骤是:1.求给定的函数的累积直方图s;2.求原图像的累积直方图G;3.求s中每一个值G中距离最小的位置index;4.求原图像每个像素通过index映射到新像素的值。
[5]吴铁洲,熊才权.直方图匹配图像增强技术的算法研究与实现[J].湖北工业大学学报,2005.
请见图6阴影去除结果示意图,a为城市遥感影像,b为城市非水域阴影掩膜影像,c为传统矩匹配法阴影校正结果,d是对阴影区域进行gamma校正的结果,d是利用本发明阴影校正后的结果。可以看出c图传统矩匹配校正中,由于匹配方法的不完善,部分区域会出现过校正和欠校正的情况;d图gamma校正也是影像亮度调整的一个常用方法,但是因为没有顾及与周围非阴影区的协调性,所以效果也不佳;e图本发明的阴影校正结果与周围非阴影区过渡自然,整体协调,并且没有出现过校正和欠校正的问题。
应当理解的是,本说明书未详细阐述的部分均属于现有技术。
应当理解的是,上述针对较佳实施例的描述较为详细,并不能因此而认为是对本发明专利保护范围的限制,本领域的普通技术人员在本发明的启示下,在不脱离本发明权利要求所保护的范围情况下,还可以做出替换或变形,均落入本发明的保护范围之内,本发明的请求保护范围应以所附权利要求为准。

Claims (8)

1.一种基于高分辨率遥感影像的城市阴影检测与去除方法,其特征在于,包括如下步骤:
步骤1,对原始高分辨率影像进行初步阴影检测,并用图割法进行区域增长,得到初始的阴影区和非阴影区;
步骤2,对原始高分辨率影像的亮度分量I进行直方图均衡化得到图IHE,并对IHE进行双边滤波;
步骤3,用步骤2中的IHE亮度图与步骤2中的双边滤波结果进行差值处理,得到细节纹理图,利用Otsu阈值法在细节图上得到纹理较弱的区域,并结合NDWI得到水域;
步骤4,步骤1中的初始阴影区域去除水域的部分,得到城市非水域阴影掩膜,并确定阴影边缘;
步骤5,计算原始影像的光谱特征值与纹理特征值,并将阴影与非阴影区分别进行SLIC影像分割,得到阴影与非阴影区分割块;
步骤6,根据步骤5中计算得到的光谱特征值与纹理特征值,对每个阴影区内的分割块,在周围临近的非阴影区中计算出特征最相近的分割块,然后分别计算对应分割块的各波段的均值和方差,
步骤7,对步骤6的计算结果,用距离作为权值整合每个阴影区与周围非阴影区所有对应分割块的波段特征值,计算出该阴影区的整体校正参数,并对阴影区进行分波段逐像素的校正;
步骤8,校正好阴影区域之后,对于步骤4中的阴影边缘,计算该阴影边缘两侧阴影分割块与对应非阴影分割块平均的均值与方差,并对阴影边缘进行分波段逐像素校正;
步骤9,对于步骤3中检测的水域部分,计算亮度直方图,获得水域的阴影部分,并对水域中阴影与非阴影部分的直方图进行匹配,得到无阴影的水域。
2.如权利要求1所述的一种基于高分辨率遥感影像的城市阴影检测与去除方法,其特征在于:步骤7的具体实现方式如下,
步骤7.1,首先计算整个阴影区域内的均值与方差,然后用距离作为权值,整合非阴影对应区域的均值与方差,如下式所示,
<mrow> <mi>C</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>=</mo> <mover> <msub> <mi>C</mi> <mi>B</mi> </msub> <mo>&amp;RightArrow;</mo> </mover> <mo>&amp;times;</mo> <mover> <mi>W</mi> <mo>&amp;RightArrow;</mo> </mover> </mrow>
<mrow> <mover> <mi>W</mi> <mo>&amp;RightArrow;</mo> </mover> <mo>=</mo> <mrow> <mo>(</mo> <mi>W</mi> <mo>(</mo> <mn>0</mn> <mo>)</mo> <mo>,</mo> <mi>W</mi> <mo>(</mo> <mn>1</mn> <mo>)</mo> <mo>,</mo> <mi>W</mi> <mo>(</mo> <mn>2</mn> <mo>)</mo> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>W</mi> <mo>(</mo> <mi>m</mi> <mo>)</mo> <mo>)</mo> </mrow> </mrow>
其中,C(i,j)表示该阴影区所对应的非阴影区整体的加权均值和方差,为非阴影区每个分割块的均值和方差组成的向量,CB(m)是每个分割块的校正参数,m是阴影区内对应的分割块数,为非阴影区不同分割块与对应阴影分割块的距离组成的权向量;
步骤7.2,按照上一步的结果,用整体阴影和非阴影区的均值和方差进行计算,按照下式进行分波段逐像素的校正,
其中,Is′(i,j)与Is(i,j)分别为阴影区校正前后各像素的值,μn与μs分别为非阴影区与阴影区整体的均值,σn-s为阴影区与非阴影区的协方差,与σs为阴影区的方差。
3.如权利要求2所述的一种基于高分辨率遥感影像的城市阴影检测与去除方法,其特征在于:步骤8的具体实现方式如下,
首先在与该阴影边缘临近的同质分割块中寻找阴影和非阴影两个分割块作为参考,然后计算两参考块的均值与方差的平均值,并用按照下式进行校正,
其中Ie′(i,j)与Ie(i,j)分别为阴影边缘校正前后各像素的值,μavg和σavg分别为阴影边缘两侧的阴影分割块与对应非阴影分割块平均的均值和方差,μe和σe分别为阴影边缘的均值和方差。
4.如权利要求3所述的一种基于高分辨率遥感影像的城市阴影检测与去除方法,其特征在于:步骤1的具体实现方式如下,
步骤1.1,对影像进行分波段的双边滤波,并利用四波段数据计算得到的B’、I、Q、A四种特征值,分别用Otsu阈值法进行处理,得到的四个结果的交集,即为初步阴影检测结果,其中B′=B/(R+G+B),R、G、B分别表示红、绿、蓝三波段数据,B′为蓝波段的归一化值,是阴影的主特征值,I为HSI色彩空间下的亮度值,Q为复合特征Q=B′-I,A为复合特征,
其中NDVI为归一化植被指数,NDVI=(NR-R)/(NR+R),t为判定植被的阈值;
步骤1.2,用图割法对阴影进行区域增长,认为待增长区是前景,非阴影区是可能背景,图割后得到的前景与可能前景都视为阴影区,背景与可能背景都视为非阴影区。
5.如权利要求4所述的一种基于高分辨率遥感影像的城市阴影检测与去除方法,其特征在于:步骤2中对IHE进行双边滤波的实现方式如下,
用双边滤波对IHE图进行处理得到滤波影像,双边滤波器中:
<mrow> <msub> <mi>I</mi> <mrow> <mi>B</mi> <mi>F</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <msub> <mi>I</mi> <mrow> <mi>H</mi> <mi>E</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mi>&amp;omega;</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mi>&amp;omega;</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>
<mrow> <msub> <mi>W</mi> <mi>r</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>e</mi> <mfrac> <mrow> <mo>-</mo> <mo>|</mo> <mo>|</mo> <msub> <mi>I</mi> <mrow> <mi>H</mi> <mi>E</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>I</mi> <mrow> <mi>H</mi> <mi>E</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>i</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msup> <mi>j</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <msub> <mi>&amp;sigma;</mi> <mi>r</mi> </msub> <mn>2</mn> </msup> </mrow> </mfrac> </msup> </mrow>
式中IBF是双边滤波后得到的影像,i,j分别为每个像素的横纵坐标,权重系数ω(i,j)=Ws(i,j)×Wr(i,j),其中Ws(i,j)和Wr(i,j)分别为IHE图的相似性因子和空间距离判定因子,σs和σr表示高斯函数标准偏差,(i,j)表示待滤波中心像素坐标,(i’,j’)表示窗口其他像素坐标,IHE(i’,j’)表示双边滤波后中心像素的亮度值。
6.如权利要求5所述的一种基于高分辨率遥感影像的城市阴影检测与去除方法,其特征在于:步骤4中通过对城市非水域阴影掩膜影像各进行一次腐蚀和一次膨胀,然后将两次的结果相减,得到的结果视为阴影边缘。
7.如权利要求6所述的一种基于高分辨率遥感影像的城市阴影检测与去除方法,其特征在于:步骤5的具体实现方式如下,
步骤5.1,用Hog算子提取影像的弱纹理强度值,并用Gabor变化提取不同尺度和方向的纹理向量式中σklkl分别是每个分割区块的小波变换系数的方差与均值,k、l分别表示小波变换的尺度和方向;
步骤5.2,将原始影像从RGB颜色空间转换到CIE-Lab颜色空间,对应每个像素的(L,a,b)颜色值和(x,y)坐标组成一个5维向量V[L,a,b,x,y];
步骤5.3,按照城市非水域阴影掩膜将影像分为阴影和非阴影两张影像,并分成K个种子点,根据步骤5.2中计算所得像素间V向量的距离大小,对每个种子点的周围空间进行K-means聚类,得到最终SLIC分割结果。
8.如权利要求7所述的一种基于高分辨率遥感影像的城市阴影检测与去除方法,其特征在于:步骤6中对每个阴影区内的分割块,在周围临近的非阴影区中计算出特征最相近的分割块,实现方式如下,
针对每个阴影区内的分割块,计算其周围非阴影区各分割块的Gabor向量间的距离dG,结合L光谱分量、Hog算子,计算相似性,并用距离作为权值,找到与该阴影分割块地物类型最接近的非阴影区分割块,
式中Dist(s,n)表示非阴影区与阴影区内分割块特征值的欧氏距离,距离最小的表示非阴影区地物类型最接近阴影区的分割块,s表示阴影区分割块,n表示非阴影区分割块,为表示两个分割块距离的权值,(Xs,Ys)、(Xn,Yn)分别为阴影分割块与非阴影分割块的中心像素坐标。
CN201711185046.5A 2017-11-23 2017-11-23 一种基于高分辨率遥感影像的城市阴影检测与去除方法 Expired - Fee Related CN107862667B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711185046.5A CN107862667B (zh) 2017-11-23 2017-11-23 一种基于高分辨率遥感影像的城市阴影检测与去除方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711185046.5A CN107862667B (zh) 2017-11-23 2017-11-23 一种基于高分辨率遥感影像的城市阴影检测与去除方法

Publications (2)

Publication Number Publication Date
CN107862667A true CN107862667A (zh) 2018-03-30
CN107862667B CN107862667B (zh) 2019-12-24

Family

ID=61703584

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711185046.5A Expired - Fee Related CN107862667B (zh) 2017-11-23 2017-11-23 一种基于高分辨率遥感影像的城市阴影检测与去除方法

Country Status (1)

Country Link
CN (1) CN107862667B (zh)

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108594255A (zh) * 2018-04-20 2018-09-28 武汉大学 一种激光测距辅助光学影像联合平差方法及系统
CN108830844A (zh) * 2018-06-11 2018-11-16 北华航天工业学院 一种基于多时相高分辨率遥感影像的设施蔬菜提取方法
CN108846806A (zh) * 2018-05-14 2018-11-20 北京洛斯达数字遥感技术有限公司 图像处理方法、图像处理装置以及记录介质
CN108846402A (zh) * 2018-05-25 2018-11-20 南京师范大学 基于多源数据的梯田田坎自动化提取方法
CN108875292A (zh) * 2018-05-16 2018-11-23 中国水利水电科学研究院 基于遥感的流域水文的仿真系统及方法
CN110163141A (zh) * 2019-05-16 2019-08-23 西安电子科技大学 基于遗传算法的卫星图像预处理方法
CN111191628A (zh) * 2020-01-06 2020-05-22 河海大学 基于决策树与特征优化的遥感影像震害建筑物识别方法
CN112559786A (zh) * 2020-12-08 2021-03-26 中国联合网络通信集团有限公司 一种光学遥感影像成像时间的确定方法及装置
CN112819720A (zh) * 2021-02-02 2021-05-18 Oppo广东移动通信有限公司 图像处理方法、装置、电子设备及存储介质
CN113012059A (zh) * 2021-02-08 2021-06-22 瞬联软件科技(北京)有限公司 一种文字图像的阴影消除方法、装置及电子设备
CN113177473A (zh) * 2021-04-29 2021-07-27 生态环境部卫星环境应用中心 遥感影像水体自动化提取方法和装置
CN113487502A (zh) * 2021-06-30 2021-10-08 中南大学 一种坑洼图像的阴影去除方法
CN117252789A (zh) * 2023-11-10 2023-12-19 中国科学院空天信息创新研究院 高分辨率遥感影像阴影重建方法、装置及电子设备
CN117541935A (zh) * 2023-11-28 2024-02-09 自然资源部国土卫星遥感应用中心 一种复杂城市环境下的资源三号遥感影像绿地提取方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101114385A (zh) * 2007-08-07 2008-01-30 深圳先进技术研究院 一种数字城市全自动生成的方法
CN103295013A (zh) * 2013-05-13 2013-09-11 天津大学 一种基于成对区域的单幅图像阴影检测方法
CN105590316A (zh) * 2015-12-11 2016-05-18 中国测绘科学研究院 面向对象的高分辨率遥感影像阴影提取方法
US9430715B1 (en) * 2015-05-01 2016-08-30 Adobe Systems Incorporated Identifying and modifying cast shadows in an image
CN105913441A (zh) * 2016-04-27 2016-08-31 四川大学 一种用于改善视频中目标检测性能的阴影去除方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101114385A (zh) * 2007-08-07 2008-01-30 深圳先进技术研究院 一种数字城市全自动生成的方法
CN103295013A (zh) * 2013-05-13 2013-09-11 天津大学 一种基于成对区域的单幅图像阴影检测方法
US9430715B1 (en) * 2015-05-01 2016-08-30 Adobe Systems Incorporated Identifying and modifying cast shadows in an image
CN105590316A (zh) * 2015-12-11 2016-05-18 中国测绘科学研究院 面向对象的高分辨率遥感影像阴影提取方法
CN105913441A (zh) * 2016-04-27 2016-08-31 四川大学 一种用于改善视频中目标检测性能的阴影去除方法

Cited By (25)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108594255B (zh) * 2018-04-20 2021-09-03 武汉大学 一种激光测距辅助光学影像联合平差方法及系统
CN108594255A (zh) * 2018-04-20 2018-09-28 武汉大学 一种激光测距辅助光学影像联合平差方法及系统
CN108846806A (zh) * 2018-05-14 2018-11-20 北京洛斯达数字遥感技术有限公司 图像处理方法、图像处理装置以及记录介质
CN108846806B (zh) * 2018-05-14 2020-11-10 北京洛斯达科技发展有限公司 图像处理方法、图像处理装置以及记录介质
CN108875292A (zh) * 2018-05-16 2018-11-23 中国水利水电科学研究院 基于遥感的流域水文的仿真系统及方法
CN108846402A (zh) * 2018-05-25 2018-11-20 南京师范大学 基于多源数据的梯田田坎自动化提取方法
CN108846402B (zh) * 2018-05-25 2022-02-11 南京师范大学 基于多源数据的梯田田坎自动化提取方法
CN108830844A (zh) * 2018-06-11 2018-11-16 北华航天工业学院 一种基于多时相高分辨率遥感影像的设施蔬菜提取方法
CN108830844B (zh) * 2018-06-11 2021-09-10 北华航天工业学院 一种基于多时相高分辨率遥感影像的设施蔬菜提取方法
CN110163141B (zh) * 2019-05-16 2023-04-07 西安电子科技大学 基于遗传算法的卫星图像预处理方法
CN110163141A (zh) * 2019-05-16 2019-08-23 西安电子科技大学 基于遗传算法的卫星图像预处理方法
CN111191628A (zh) * 2020-01-06 2020-05-22 河海大学 基于决策树与特征优化的遥感影像震害建筑物识别方法
CN112559786A (zh) * 2020-12-08 2021-03-26 中国联合网络通信集团有限公司 一种光学遥感影像成像时间的确定方法及装置
CN112559786B (zh) * 2020-12-08 2024-03-15 中国联合网络通信集团有限公司 一种光学遥感影像成像时间的确定方法及装置
CN112819720A (zh) * 2021-02-02 2021-05-18 Oppo广东移动通信有限公司 图像处理方法、装置、电子设备及存储介质
CN112819720B (zh) * 2021-02-02 2023-10-03 Oppo广东移动通信有限公司 图像处理方法、装置、电子设备及存储介质
CN113012059A (zh) * 2021-02-08 2021-06-22 瞬联软件科技(北京)有限公司 一种文字图像的阴影消除方法、装置及电子设备
CN113177473B (zh) * 2021-04-29 2021-11-16 生态环境部卫星环境应用中心 遥感影像水体自动化提取方法和装置
CN113177473A (zh) * 2021-04-29 2021-07-27 生态环境部卫星环境应用中心 遥感影像水体自动化提取方法和装置
CN113487502A (zh) * 2021-06-30 2021-10-08 中南大学 一种坑洼图像的阴影去除方法
CN113487502B (zh) * 2021-06-30 2022-05-03 中南大学 一种坑洼图像的阴影去除方法
CN117252789A (zh) * 2023-11-10 2023-12-19 中国科学院空天信息创新研究院 高分辨率遥感影像阴影重建方法、装置及电子设备
CN117252789B (zh) * 2023-11-10 2024-02-02 中国科学院空天信息创新研究院 高分辨率遥感影像阴影重建方法、装置及电子设备
CN117541935A (zh) * 2023-11-28 2024-02-09 自然资源部国土卫星遥感应用中心 一种复杂城市环境下的资源三号遥感影像绿地提取方法
CN117541935B (zh) * 2023-11-28 2024-04-30 自然资源部国土卫星遥感应用中心 一种复杂城市环境下的资源三号遥感影像绿地提取方法

Also Published As

Publication number Publication date
CN107862667B (zh) 2019-12-24

Similar Documents

Publication Publication Date Title
CN107862667A (zh) 一种基于高分辨率遥感影像的城市阴影检测与去除方法
CN104851113B (zh) 多分辨率遥感影像的城市植被自动提取方法
CN106651872B (zh) 基于Prewitt算子的路面裂缝识别方法及系统
CN102722891B (zh) 一种图像显著度检测的方法
CN103559500B (zh) 一种基于光谱与纹理特征的多光谱遥感图像地物分类方法
CN109191432B (zh) 基于域变换滤波多尺度分解的遥感图像云检测方法
CN111915704A (zh) 一种基于深度学习的苹果分级识别方法
CN107358585B (zh) 基于分数阶微分及暗原色先验的雾天图像增强方法
CN107240084B (zh) 一种单幅图像去雨方法及装置
US20100008576A1 (en) System and method for segmentation of an image into tuned multi-scaled regions
CN105761266A (zh) 从遥感图像中提取矩形建筑物的方法
CN104217440B (zh) 一种从遥感图像中提取建成区的方法
CN112488046B (zh) 一种基于无人机高分辨率影像的车道线提取方法
CN108765347A (zh) 一种适合遥感影像的色彩增强方法
CN106294705A (zh) 一种批量遥感影像预处理方法
CN108875747A (zh) 一种基于机器视觉的小麦不完善粒识别方法
CN111275652B (zh) 一种去除城市遥感图像中雾霾的方法
CN105913421A (zh) 基于自适应形状暗通道的遥感图像云检测方法
CN110175556B (zh) 基于Sobel算子的遥感图像云检测方法
CN113077486B (zh) 一种山区植被覆盖率监测方法及系统
CN107818303A (zh) 无人机油气管线影像自动对比分析方法、系统及软件存储器
CN107992856A (zh) 城市场景下的高分遥感建筑物阴影检测方法
CN115100077B (zh) 一种图像增强方法与装置
CN106650663A (zh) 建筑物真伪变化的判定方法及含此方法的伪变化去除方法
CN106875407A (zh) 一种结合形态学和标记控制的无人机影像林冠分割方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20191224

Termination date: 20211123