CN103761717B - 一种基于全色遥感影像的城市水体提取方法 - Google Patents
一种基于全色遥感影像的城市水体提取方法 Download PDFInfo
- Publication number
- CN103761717B CN103761717B CN201410037069.1A CN201410037069A CN103761717B CN 103761717 B CN103761717 B CN 103761717B CN 201410037069 A CN201410037069 A CN 201410037069A CN 103761717 B CN103761717 B CN 103761717B
- Authority
- CN
- China
- Prior art keywords
- scale
- row
- sigma
- col
- window
- 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.)
- Expired - Fee Related
Links
Landscapes
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
本发明公开一种基于全色遥感影像的城市水体提取方法,通过遥感影像预处理、最优尺度选取、均变纹理生成和最优阈值分割四个步骤最终获得所需提取的城市地区的水体信息。本发明计算方法简单,运算量较小,并且最终提取结果较为精准,能够应用于城市规划、环境科学、地理信息制图等多个领域。
Description
技术领域
本发明涉及一种图像信息提取方法,具体涉及一种基于全色遥感影像的城市水体提取方法。
背景技术
城市水体作为一种重要的空间地物和自然资源,对于城市规划、环境科学和城市遥感制图都发挥着不可或缺的作用。近年来,随着卫星遥感技术的发展和各种空间地物提取算法的出现,水体提取方法不断进步,在各项研究中得到了广泛应用。目前,常用的水体提取算法都是基于多光谱遥感影像如SPOT和TM遥感影像,使用光谱变换生成各种水体指数,然后通过监督分类或者阈值分割进行水体提取,但是少有基于全色遥感影像和针对城市地区水体提取的算法。
针对水体提取问题,目前少有针对全色遥感影像和城市区域的提取算法,并且一般仅利用了光谱信息,没有考虑到纹理特征,没能够充分利用遥感影像的结构化特征,未能够最大化利用遥感中丰富的信息。
因此,提出基于全色遥感影像并且针对城市区域的水体提取算法非常有意义,也会对同类研究起到积极作用。
发明内容
发明目的:为解决现有技术中存在的不足,本发明提供了一种基于全色遥感影像的城市水体提取方法。
技术方案:本发明的一种基于全色遥感影像的城市水体提取方法,包括以下步骤:
(1)遥感影像预处理:利用中值滤波抑制遥感影像中的噪声;
(2)最优尺度选取:在经步骤(1)预处理的图像中,选定最小尺度、最大尺度以及尺度变化步长,依次计算全局平均方差,求得最优空间尺度;
(3)均变纹理生成:依据步骤(2)中所得的最优空间尺度,遍历遥感影像中的各个像素,分别计算均变方差,进而生成遥感影像的均变方差纹理图;
(4)最优阈值分割:基于步骤(3)中所得遥感影像的均变方差纹理图,人工选定分割阈值,对方差纹理图进行分割,选取大于给定面积阈值的斑块,进行闭运算,获得最终水体信息。
进一步的,所述步骤(2)中确定最优尺度的详细步骤如下:
(1)确定最小尺度scale_min,最大尺度scale_max以及尺度变化步长scale_step,考虑到尺度过大之后会造成图像中出现明显的边缘效应,因此过大的尺度不予考虑,例如本发明中各项参数选取如下:scale_min=1,scale_max=11,scale_step=2;
(2)对于第k个空间尺度scalek=scale_min+scale_step*(k-1),其中,k=1,2,…,(scale_max-scale_min)/scale_step+1,在原始遥感影像上设置一个大小为(2*scalek+1)×(2*scalek+1)的窗口,不断移动该窗口,对原始图像im中第i行第j列的像元im(i,j),窗口中的所有像元光谱值记为im(i-scalek:i-scalek,j-scalek:j-scalek),然后计算窗口中所有像元光谱值的方差imvar(i,j,k)
上述公式中,ii代表某个像元在窗口中的行号,jj代表某个像元在窗口中的列号,将imvar(i,j,k)作为该窗口中心像元im(i,j)在空间尺度scalek下的方差,并按照上述方法计算处全图中所有像元的方差,并求其平均值imvar0
(3)比较k个空间尺度下的imvar0,选择imvar0最大时所对应的空间尺度scalebest作为遥感影像的最优空间尺度;
进一步的,上述步骤(3)中遥感影像均变方差纹理图的具体生成方法如下:
(1)依据上述步骤得到的最优空间尺度scalebest,在原始遥感影像上设置一个大小为(2*scalebest+1)×(2*scalebest+1)的窗口;
(2)不断移动该窗口,对于像元im(i,j),窗口中的所有像元光谱值记为im(i-scalebest:i-scalebest,j-scalebest:j-scalebest),每个窗口中的像元总个数记作N=(2*scalebest+1)*(2*scalebest+1),将其按照列方向排列记作imN,将窗口中所有像元的行坐标和列坐标依次记作rowN和colN;
(3)依此对每个窗口中的所有像元的光谱值进行最小二乘平面拟合,拟合方程为:Ax+By+Cz+1=0,将待拟合的N个点表示成以下矩阵形式:
依据最小二乘法则得到平面拟合系数为:
其中,h代表N个像元中的第h个像元;
(4)计算窗口中所有像元光谱值至拟合平面的距离:
将dvar作为该窗口中心像元的均变方差纹理,遍历所有点,即生成均变方差纹理图imvar_texture;
遍历所有像素点,搜索均变方差纹理图imvar_texture的最大值和最小值,分别记作imvar_texture_max和imvar_texture_min,对每个像元进行如下操作用于将纹理值拉伸至[0,255]:
进一步的,上述步骤(4)中最优阈值分割的具体步骤如下:
(1)在所得的方差纹理图中,确定一个分割阈值;
(2)在步骤(1)中确定的分割阈值中,根据实际情况确定一个面积阈值用于抑制微小斑块,依次遍历每一个斑块,计算其面积,如果该斑块的面积小于该阈值,则去除该斑块,否则保留该斑块;
(3)基于上述结果,进行闭运算,用于填充水体中的微小空洞,至此,已经得到最终的水体信息,如果有必要需要显示面积较大的水体,可以再设置一个面积阈值,保留面积大于给定阈值的水体区域。
有益效果:本发明的一种基于全色遥感影像的水体提取方法,基于水体因其流动性而造成的水体区域光谱响应局部变化率较小的原理,并提出均变方差纹理的概念用以度量光谱值局部变化率的均一性,充分利用纹理信息用以消除和水体相似光谱地物诸如阴影和居民地等的影响,使得水体部分的纹理特征和其他区域的纹理值之间差异更大,最后利用最优阈值分割方法高效快速地提取水体信息,为城市水体信息的提取提出了一种高效、快速的新方法。本发明计算方法简单,运算复杂度较低,并且最终提取结果较为精准,能够应用于城市规划、环境科学、地理信息制图等多个领域。
附图说明
图1为本发明的流程图;
图2为本发明中提取的水体的原始遥感影像;
图3为图1的均变方差纹理图;
图4为图3的分割结果示意图;
图5为图4去除微小噪声后的结果示意图;
图6为对图5结果进行闭运算得到的结果示意图;
图7为对图6中较大面积水体的显示结果示意图。
具体实施方式
下面对本发明技术方案进行详细说明,但是本发明的保护范围不局限于所述实施例。
下面对本发明技术方案结合附图进行详细说明。
本发明的一种基于全色遥感影像的城市水体提取方法,包括以下步骤:
(1)遥感影像预处理:利用中值滤波抑制遥感影像中的噪声;
(2)最优尺度选取:在经步骤(1)预处理的图像中,选定最小尺度、最大尺度以及尺度变化步长,依次计算全局平均方差,求得最优空间尺度;
(3)均变纹理生成:依据步骤(2)中所得的最优空间尺度,遍历遥感影像中的各个像素,分别计算均变方差,进而生成遥感影像的均变方差纹理图;
(4)最优阈值分割:基于步骤(3)中所得遥感影像的均变方差纹理图,人工选定分割阈值,对方差纹理图进行分割,选取大于给定面积阈值的斑块,进行闭运算,获得最终水体信息。
本发明中步骤(2)中确定最优尺度的详细算法如下:
(1)确定最小尺度scale_min,最大尺度scale_max以及尺度变化步长scale_step;
(2)对于第k个空间尺度scalek=scale_min+scale_step*(k-1),其中,k=1,2,…,(scale_max-scale_min)/scale_step+1,在原始遥感影像上设置一个大小为(2*scalek+1)×(2*scalek+1)的窗口,不断移动该窗口,对于像元im(i,j),窗口中的所有像元光谱值记为im(i-scalek:i-scalek,j-scalek:j-scalek),然后计算窗口中所有像元光谱值的方差imvar(i,j,k)
将imvar(i,j,k)作为该窗口中心像元im(i,j)在空间尺度scalek下的方差,并按照上述方法计算处全图中所有像元的方差,并求其平均值imvar0
(3)比较k个空间尺度下的imvar0,选择imvar0最大时所对应的空间尺度scalebest作为遥感影像的最优空间尺度;
本发明中步骤(3)遥感影像均变方差纹理图的具体生成方法如下:
(1)依据上述步骤得到的最优空间尺度scalebest,在原始遥感影像上设置一个大小为(2*scalebest+1)×(2*scalebest+1)的窗口;
(2)不断移动该窗口,对于像元im(i,j),窗口中的所有像元光谱值记为im(i-scalebest:i-scalebest,j-scalebest:j-scalebest),每个窗口中的像元总个数记作N=(2*scalebest+1)*(2*scalebest+1),将其按照列方向排列记作imN,将窗口中所有像元的行坐标和列坐标依次记作rowN和colN;
(3)依此对每个窗口中的所有像元的光谱值进行最小二乘平面拟合,拟合方程为:Ax+By+Cz+1=0,将待拟合的N个点表示成以下矩阵形式:
依据最小二乘法则得到平面拟合系数为:
计算窗口中所有像元光谱值至拟合平面的距离:
对该窗口中的所有距离求得方差dvar
将dvar作为该窗口中心像元的均变方差纹理,遍历所有点,即生成均变方差纹理图imvar_texture;
遍历所有像素点,搜索均变方差纹理图imvar_texture的最大值和最小值,分别记作imvar_texture_max和imvar_texture_min,对每个像素点作如下操作:
本发明中步骤(4)中最优阈值分割的具体步骤如下:
(1)在所得的方差纹理图中,确定一个分割阈值;
(2)在步骤(1)中确定的分割阈值中,根据实际情况确定一个面积阈值用于抑制微小斑块,依次遍历每一个斑块,计算其面积,如果该斑块的面积小于该阈值,则去除该斑块,否则保留该斑块;
(3)基于上述结果,进行闭运算,用于填充水体中的微小空洞,至此,已经得到最终的水体信息,如果有必要需要显示面积较大的水体,可以再设置一个面积阈值,保留面积大于给定阈值的水体区域。
Claims (3)
1.一种基于全色遥感影像的城市水体提取方法,其特征在于包括以下步骤:
(1)遥感影像预处理:利用中值滤波抑制遥感影像中的噪声;
(2)最优尺度选取:在经步骤(1)预处理的图像中,选定最小尺度、最大尺度以及尺度变化步长,依次计算全局平均方差,求得最优空间尺度;
(3)均变纹理生成:依据步骤(2)中所得的最优空间尺度,遍历遥感影像中的各个像素,分别计算均变方差,进而生成遥感影像的均变方差纹理图,具体生成方法如下:
(31)依据所得的最优空间尺度scalebest,在原始遥感影像上设置一个大小为(2*scalebest+1)×(2*scalebest+1)的窗口;
(32)不断移动该窗口,对于像元im(i,j),窗口中的所有像元光谱值记为im(i-scalek:i+scalek,j-scalek:j+scalek),每个窗口中的像元总个数记作N=(2*scalebest+1)*(2*scalebest+1);
然后将其按照列方向排列记作im1,im2,…,imN,将窗口中所有像元的行坐标和列坐标依次记作rowN和colN;
(33)依次对每个窗口中的所有像元的光谱值进行最小二乘平面拟合,拟合方程为:Ax+By+Cz+1=0,将待拟合的N个像元表示成以下矩阵形式:
依据最小二乘法则得到平面拟合系数为:
上述公式中,h代表N个像元中的第h个像元;
(34)计算窗口中所有像元光谱值至拟合平面的距离:
计算该窗口中的各个像元所对应的距离,求得方差dvar
将dvar作为该窗口中心像元的均变方差纹理,遍历所有点,即生成均变方差纹理图imvar_texture;
遍历所有像素点,搜索均变方差纹理图imvar_texture的最大值和最小值,分别记作imvar_texture_max和imvar_texture_min,对每个像素点作如下操作用于将纹理值拉伸至[0,255]:
(4)最优阈值分割:基于步骤(3)中所得遥感影像的均变方差纹理图,人工选定分割阈值,对方差纹理图进行分割,选取大于给定面积阈值的斑块,进行闭运算,获得最终水体信息。
2.根据权利要求1所述的基于全色遥感影像的城市水体提取方法,其特征在于:所述步骤(2)中确定最优尺度的详细步骤如下:
(21)确定最小尺度scale_min,最大尺度scale_max以及尺度变化步长scale_step;
(22)对于第k个空间尺度scalek=scale_min+scale_step*(k-1),其中,k=1,2,…,(scale_max-scale_min)/scale_step+1,在原始遥感影像上设置一个大小为(2*scalek+1)×(2*scalek+1)的窗口,不断移动该窗口,对原始图像im中第i行第j列的像元im(i,j),窗口中的所有像元光谱值记为im(i-scalek:i+scalek,j-scalek:j+scalek),然后计算窗口中所有像元光谱值的方差imvar(i,j,k):
上述公式中,ii代表某个像元在窗口中的行号,jj代表某个像元在窗口中的列号;将imvar(i,j,k)作为该窗口中心像元im(i,j)在空间尺度scalek下的方差,并按照上述方法计算出全图中所有像元的方差,并求其平均值imvar0;
(23)比较k个空间尺度下的imvar0,选择imvar0最大时所对应的空间尺度scalebest作为遥感影像的最优空间尺度,scalebest=1,则移动窗口大小为3×3。
3.根据权利要求1所述的基于全色遥感影像的城市水体提取方法,其特征在于:所述步骤(4)中最优阈值分割的具体步骤如下:
(41)将方差纹理图中,由均变方差纹理图的灰度直方图确定影像中的水体和非水体的纹理值分布,从而获得分割阈值;
(42)在步骤(1)中确定的分割阈值中,根据实际情况确定一个面积阈值用于抑制微小斑块,依次遍历每一个斑块,计算其面积,如果该斑块的面积小于该阈值,则去除该斑块,否则保留该斑块;
(43)基于上述结果,进行闭运算,用于填充水体中的微小空洞,至此,已经得到最终的水体信息;通过另设一个面积阈值,保留面积大于给定阈值的水体区域,进而显示面积较大的水体。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410037069.1A CN103761717B (zh) | 2014-01-26 | 2014-01-26 | 一种基于全色遥感影像的城市水体提取方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410037069.1A CN103761717B (zh) | 2014-01-26 | 2014-01-26 | 一种基于全色遥感影像的城市水体提取方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN103761717A CN103761717A (zh) | 2014-04-30 |
CN103761717B true CN103761717B (zh) | 2016-06-22 |
Family
ID=50528950
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410037069.1A Expired - Fee Related CN103761717B (zh) | 2014-01-26 | 2014-01-26 | 一种基于全色遥感影像的城市水体提取方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN103761717B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104318051B (zh) * | 2014-09-15 | 2017-06-16 | 中国水利水电科学研究院 | 基于规则的大范围水体信息遥感自动提取系统及方法 |
CN106023133B (zh) * | 2016-04-26 | 2019-03-19 | 武汉大学 | 一种基于多特征联合处理的高分辨率遥感影像水体提取方法 |
CN107507200B (zh) * | 2017-08-31 | 2021-03-30 | 电子科技大学 | 一种基于连通检测与噪声抑制的sar影像高精度大范围水域提取方法 |
CN107527331B (zh) * | 2017-09-26 | 2020-10-02 | 南京大学 | 基于双冒泡法的极地冰山遥感识别方法 |
CN108363949B (zh) * | 2017-12-27 | 2020-08-28 | 二十一世纪空间技术应用股份有限公司 | 一种基于物候分析的棉花遥感监测方法 |
CN109141403B (zh) * | 2018-08-01 | 2021-02-02 | 上海航天控制技术研究所 | 一种星敏感器小窗口访问的图像处理系统及其方法 |
CN114220017A (zh) * | 2022-02-23 | 2022-03-22 | 广东省科学院广州地理研究所 | 遥感数据尺度自适应调整方法、装置、存储介质及设备 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102284827A (zh) * | 2011-07-19 | 2011-12-21 | 大连理工大学 | 基于定位面拟合的平行度修整法 |
CN102855494B (zh) * | 2012-08-20 | 2015-08-12 | 中国测绘科学研究院 | 一种卫星遥感影像的水体提取方法及装置 |
CN103400151B (zh) * | 2013-08-16 | 2016-07-27 | 武汉大学 | 一体化的光学遥感影像与gis自动配准与水体提取方法 |
-
2014
- 2014-01-26 CN CN201410037069.1A patent/CN103761717B/zh not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
CN103761717A (zh) | 2014-04-30 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN103761717B (zh) | 一种基于全色遥感影像的城市水体提取方法 | |
CN103871076B (zh) | 基于光流法和超像素分割的运动目标提取方法 | |
CN103400151B (zh) | 一体化的光学遥感影像与gis自动配准与水体提取方法 | |
CN110119728A (zh) | 基于多尺度融合语义分割网络的遥感图像云检测方法 | |
CN107220949A (zh) | 公路监控视频中运动车辆阴影的自适应消除方法 | |
CN103871062B (zh) | 一种基于超像素描述的月面岩石检测方法 | |
CN106530247B (zh) | 一种基于结构信息的多尺度图像修复方法 | |
CN103455991A (zh) | 一种多聚焦图像融合方法 | |
CN105869178A (zh) | 一种基于多尺度组合特征凸优化的复杂目标动态场景无监督分割方法 | |
CN106952288A (zh) | 基于卷积特征和全局搜索检测的长时遮挡鲁棒跟踪方法 | |
CN106295637A (zh) | 一种基于深度学习与强化学习的车辆识别方法 | |
CN103745453B (zh) | 基于Google Earth遥感影像的城镇信息提取方法 | |
CN101447076A (zh) | 一种web图像中感兴趣区域的分割方法 | |
CN103106658A (zh) | 一种海岛、礁岸线快速提取方法 | |
CN103095996A (zh) | 基于时空显著性检测的多传感器视频融合方法 | |
CN103400383A (zh) | 基于nsct和压缩投影的sar图像变化检测方法 | |
CN104615990A (zh) | 一种基于怀柔全日面单色像的太阳黑子自动识别方法 | |
CN104318100A (zh) | 一种基于特征敏感投影算子的厚度点云薄化方法 | |
CN104091339A (zh) | 一种图像快速立体匹配方法及装置 | |
CN103955945A (zh) | 基于双目视差和活动轮廓的自适应彩色图像分割方法 | |
CN104751111A (zh) | 识别视频中人体行为的方法和系统 | |
CN105138992A (zh) | 一种基于区域主动轮廓模型的海岸线检测方法 | |
CN103971354A (zh) | 低分辨率红外图像重建高分辨率红外图像的方法 | |
CN105469428B (zh) | 一种基于形态学滤波和svd的弱小目标检测方法 | |
CN105184802A (zh) | 一种图像处理的方法及装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20160622 Termination date: 20190126 |
|
CF01 | Termination of patent right due to non-payment of annual fee |