一种基于全色遥感影像的城市水体提取方法

摘要

本发明公开一种基于全色遥感影像的城市水体提取方法，通过遥感影像预处理、最优尺度选取、均变纹理生成和最优阈值分割四个步骤最终获得所需提取的城市地区的水体信息。本发明计算方法简单，运算量较小，并且最终提取结果较为精准，能够应用于城市规划、环境科学、地理信息制图等多个领域。

说明书

技术领域

本发明涉及一种图像信息提取方法，具体涉及一种基于全色遥感影像的城市水体提取方法。

背景技术

城市水体作为一种重要的空间地物和自然资源，对于城市规划、环境科学和城市遥感制图都发挥着不可或缺的作用。近年来，随着卫星遥感技术的发展和各种空间地物提取算法的出现，水体提取方法不断进步，在各项研究中得到了广泛应用。目前，常用的水体提取算法都是基于多光谱遥感影像如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个空间尺度scale_k=scale_min+scale_step*(k-1)，其中，k=1,2,…,(scale_max-scale_min)/scale_step+1，在原始遥感影像上设置一个大小为(2*scale_k+1)×(2*scale_k+1)的窗口，不断移动该窗口，对原始图像im中第i行第j列的像元im(i,j)，窗口中的所有像元光谱值记为im(i-scale_k:i-scale_k,j-scale_k:j-scale_k)，然后计算窗口中所有像元光谱值的方差im_var(i,j,k)

$>{\mathrm{im}}_{\mathrm{mean}}(i,j,k)=\frac{\underset{\mathrm{ii}=1}{\overset{2*{\mathrm{scale}}_k+1}{\Sigma }}\underset{\mathrm{jj}=1}{\overset{2*{\mathrm{scale}}_k+1}{\Sigma }}\mathrm{im}(i-{\mathrm{scale}}_k+\mathrm{ii},j-{\mathrm{scale}}_k+\mathrm{jj})}{(2*{\mathrm{scale}}_k+1)*(2*{\mathrm{scale}}_k+1)}$>

$>{\mathrm{im}}_{\mathrm{var}}(i,j,k)=\frac{\underset{\mathrm{ii}=1}{\overset{2*{\mathrm{scale}}_k+1}{\Sigma }}\underset{\mathrm{jj}=1}{\overset{2*{\mathrm{scale}}_k+1}{\Sigma }}{(\mathrm{im}(i-{\mathrm{scale}}_k+\mathrm{ii},j-{\mathrm{scale}}_k+\mathrm{jj})-{\mathrm{im}}_{\mathrm{mean}})}^2}{(2*{\mathrm{scale}}_k+1)*(2*{\mathrm{scale}}_k+1)}$>

上述公式中，ii代表某个像元在窗口中的行号，jj代表某个像元在窗口中的列号，将im_var(i,j,k)作为该窗口中心像元im(i,j)在空间尺度scale_k下的方差，并按照上述方法计算处全图中所有像元的方差，并求其平均值im_var0

（3）比较k个空间尺度下的im_var0，选择im_var0最大时所对应的空间尺度scale_best作为遥感影像的最优空间尺度；

进一步的，上述步骤（3）中遥感影像均变方差纹理图的具体生成方法如下：

（1）依据上述步骤得到的最优空间尺度scale_best，在原始遥感影像上设置一个大小为(2*scale_best+1)×(2*scale_best+1)的窗口；

（2）不断移动该窗口，对于像元im(i,j)，窗口中的所有像元光谱值记为im(i-scale_best:i-scale_best,j-scale_best:j-scale_best)，每个窗口中的像元总个数记作N=(2*scale_best+1)*(2*scale_best+1)，将其按照列方向排列记作im_N，将窗口中所有像元的行坐标和列坐标依次记作row_N和col_N；

（3）依此对每个窗口中的所有像元的光谱值进行最小二乘平面拟合，拟合方程为：Ax+By+Cz+1=0，将待拟合的N个点表示成以下矩阵形式：

$>\left(\begin{array}{ccc}{\mathrm{row}}_1& {\mathrm{col}}_1& {\mathrm{im}}_1\\ {\mathrm{row}}_2& {\mathrm{col}}_2& {\mathrm{im}}_2\\ ·& ·& ·\\ ·& ·& ·\\ ·& ·& ·\\ {\mathrm{row}}_N& {\mathrm{col}}_N& {\mathrm{im}}_N\end{array}\right)\left(\begin{array}{c}A\\ B\\ C\end{array}\right)=\left(\begin{array}{c}-1\\ -1\\ -1\end{array}\right)$>

依据最小二乘法则得到平面拟合系数为：

$>\left(\begin{array}{c}A\\ B\\ C\end{array}\right)={\left(\begin{array}{ccc}{\mathrm{\Sigma row}}_h^2& \Sigma {\mathrm{row}}_h{\mathrm{col}}_h& \Sigma {\mathrm{row}}_h\mathrm{in}({\mathrm{wor}}_h,{\mathrm{col}}_h)\\ \Sigma {\mathrm{col}}_h{\mathrm{row}}_h& {\mathrm{\Sigma col}}_h^2& \Sigma {\mathrm{col}}_h\mathrm{im}({\mathrm{row}}_h,{\mathrm{col}}_h)\\ \mathrm{\Sigma im}({\mathrm{row}}_h,{\mathrm{col}}_h){\mathrm{row}}_h& \mathrm{\Sigma im}({\mathrm{row}}_h,{\mathrm{col}}_h){\mathrm{col}}_h& \mathrm{\Sigma im}{({\mathrm{row}}_h,{\mathrm{col}}_h)}^2\end{array}\right)}^{-1}\left(\begin{array}{c}\Sigma {\mathrm{row}}_h\\ \Sigma {\mathrm{col}}_h\\ \mathrm{\Sigma im}({\mathrm{row}}_h,{\mathrm{col}}_h)\end{array}\right)$>

其中，h代表N个像元中的第h个像元；

（4）计算窗口中所有像元光谱值至拟合平面的距离：$>d_k=\frac{\mathrm{abs}({\mathrm{Arow}}_h+{\mathrm{Bcol}}_h+\mathrm{Cim}({\mathrm{row}}_h,{\mathrm{col}}_h)+1)}{\sqrt{A^2+B^2+C^2}},h=\mathrm{1,2},···,N,$>其中abs(*)是绝对值算子，计算该窗口中的各个像元所对应的距离，求得方差d_var

$>d_{\mathrm{mean}}=\frac{\underset{h=1}{\overset{N}{\Sigma }}{d}_n}{N}$>

$>d_{\mathrm{var}}=\frac{\underset{h=1}{\overset{N}{\Sigma }}{(d_h-d_{\mathrm{mean}})}^2}{N}$>

将d_var作为该窗口中心像元的均变方差纹理，遍历所有点，即生成均变方差纹理图im_{var_texture}；

遍历所有像素点，搜索均变方差纹理图im_{var_texture}的最大值和最小值，分别记作im_{var_texture}_max和im_{var_texture}_min，对每个像元进行如下操作用于将纹理值拉伸至[0,255]：

$>{\mathrm{im}}_{\mathrm{var}\_\mathrm{texture}}(i,j)=255*\frac{({\mathrm{im}}_{\mathrm{var}\_\mathrm{texture}}(i,j)-{\mathrm{im}}_{\mathrm{var}\_\mathrm{texture}}-\min )}{({\mathrm{im}}_{\mathrm{var}\_\mathrm{texture}}-\max -{\mathrm{im}}_{\mathrm{var}\_\mathrm{texture}}-\min )}.$>

进一步的，上述步骤（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个空间尺度scale_k=scale_min+scale_step*(k-1)，其中，k=1,2,…,(scale_max-scale_min)/scale_step+1，在原始遥感影像上设置一个大小为(2*scale_k+1)×(2*scale_k+1)的窗口，不断移动该窗口，对于像元im(i,j)，窗口中的所有像元光谱值记为im(i-scale_k:i-scale_k,j-scale_k:j-scale_k)，然后计算窗口中所有像元光谱值的方差im_var(i,j,k)

$>{\mathrm{im}}_{\mathrm{mean}}(i,j,k)=\frac{\underset{\mathrm{ii}=1}{\overset{2*{\mathrm{scale}}_k+1}{\Sigma }}\underset{\mathrm{jj}=1}{\overset{2*{\mathrm{scale}}_k+1}{\Sigma }}\mathrm{im}(i-{\mathrm{scale}}_k+\mathrm{ii},j-{\mathrm{scale}}_k+\mathrm{jj})}{(2*{\mathrm{scale}}_k+1)*(2*{\mathrm{scale}}_k+1)}$>

$>{\mathrm{im}}_{\mathrm{var}}(i,j,k)=\frac{\underset{\mathrm{ii}=1}{\overset{2*{\mathrm{scale}}_k+1}{\Sigma }}\underset{\mathrm{jj}=1}{\overset{2*{\mathrm{scale}}_k+1}{\Sigma }}{(\mathrm{im}(i-{\mathrm{scale}}_k+\mathrm{ii},j-{\mathrm{scale}}_k+\mathrm{jj})-{\mathrm{im}}_{\mathrm{mean}})}^2}{(2*{\mathrm{scale}}_k+1)*(2*{\mathrm{scale}}_k+1)}$>

将im_var(i,j,k)作为该窗口中心像元im(i,j)在空间尺度scale_k下的方差，并按照上述方法计算处全图中所有像元的方差，并求其平均值im_var0

（3）比较k个空间尺度下的im_var0，选择im_var0最大时所对应的空间尺度scale_best作为遥感影像的最优空间尺度；

本发明中步骤（3）遥感影像均变方差纹理图的具体生成方法如下：

（1）依据上述步骤得到的最优空间尺度scale_best，在原始遥感影像上设置一个大小为(2*scale_best+1)×(2*scale_best+1)的窗口；

（2）不断移动该窗口，对于像元im(i,j)，窗口中的所有像元光谱值记为im(i-scale_best:i-scale_best,j-scale_best:j-scale_best)，每个窗口中的像元总个数记作N=(2*scale_best+1)*(2*scale_best+1)，将其按照列方向排列记作im_N，将窗口中所有像元的行坐标和列坐标依次记作row_N和col_N；

（3）依此对每个窗口中的所有像元的光谱值进行最小二乘平面拟合，拟合方程为：Ax+By+Cz+1=0，将待拟合的N个点表示成以下矩阵形式：

$>\left(\begin{array}{ccc}{\mathrm{row}}_1& {\mathrm{col}}_1& {\mathrm{im}}_1\\ {\mathrm{row}}_2& {\mathrm{col}}_2& {\mathrm{im}}_2\\ ·& ·& ·\\ ·& ·& ·\\ ·& ·& ·\\ {\mathrm{row}}_N& {\mathrm{col}}_N& {\mathrm{im}}_N\end{array}\right)\left(\begin{array}{c}A\\ B\\ C\end{array}\right)=\left(\begin{array}{c}-1\\ -1\\ -1\end{array}\right)$>

依据最小二乘法则得到平面拟合系数为：

$>\left(\begin{array}{c}A\\ B\\ C\end{array}\right)={\left(\begin{array}{ccc}{\mathrm{\Sigma row}}_h^2& \Sigma {\mathrm{row}}_h{\mathrm{col}}_h& \Sigma {\mathrm{row}}_h\mathrm{in}({\mathrm{wor}}_h,{\mathrm{col}}_h)\\ \Sigma {\mathrm{col}}_h{\mathrm{row}}_h& {\mathrm{\Sigma col}}_h^2& \Sigma {\mathrm{col}}_h\mathrm{im}({\mathrm{row}}_h,{\mathrm{col}}_h)\\ \mathrm{\Sigma im}({\mathrm{row}}_h,{\mathrm{col}}_h){\mathrm{row}}_h& \mathrm{\Sigma im}({\mathrm{row}}_h,{\mathrm{col}}_h){\mathrm{col}}_h& \mathrm{\Sigma im}{({\mathrm{row}}_h,{\mathrm{col}}_h)}^2\end{array}\right)}^{-1}\left(\begin{array}{c}\Sigma {\mathrm{row}}_h\\ \Sigma {\mathrm{col}}_h\\ \mathrm{\Sigma im}({\mathrm{row}}_h,{\mathrm{col}}_h)\end{array}\right)$>

计算窗口中所有像元光谱值至拟合平面的距离：$>d_k=\frac{\mathrm{abs}({\mathrm{Arow}}_h+{\mathrm{Bcol}}_h+\mathrm{Cim}({\mathrm{row}}_h,{\mathrm{col}}_h)+1)}{\sqrt{A^2+B^2+C^2}},h=\mathrm{1,2},···,N,$>

对该窗口中的所有距离求得方差d_var

$>d_{\mathrm{mean}}=\frac{\underset{h=1}{\overset{N}{\Sigma }}{d}_n}{N}$>

$>d_{\mathrm{var}}=\frac{\underset{h=1}{\overset{N}{\Sigma }}{(d_h-d_{\mathrm{mean}})}^2}{N}$>

将d_var作为该窗口中心像元的均变方差纹理，遍历所有点，即生成均变方差纹理图im_{var_texture}；

遍历所有像素点，搜索均变方差纹理图im_{var_texture}的最大值和最小值，分别记作im_{var_texture}_max和im_{var_texture}_min，对每个像素点作如下操作：

$>{\mathrm{im}}_{\mathrm{var}\_\mathrm{texture}}(i,j)=255*\frac{({\mathrm{im}}_{\mathrm{var}\_\mathrm{texture}}(i,j)-{\mathrm{im}}_{\mathrm{var}\_\mathrm{texture}}-\min )}{({\mathrm{im}}_{\mathrm{var}\_\mathrm{texture}}-\max -{\mathrm{im}}_{\mathrm{var}\_\mathrm{texture}}-\min )}.$>

本发明中步骤（4）中最优阈值分割的具体步骤如下：

（1）在所得的方差纹理图中，确定一个分割阈值；

（2）在步骤（1）中确定的分割阈值中，根据实际情况确定一个面积阈值用于抑制微小斑块，依次遍历每一个斑块，计算其面积，如果该斑块的面积小于该阈值，则去除该斑块，否则保留该斑块；

（3）基于上述结果，进行闭运算，用于填充水体中的微小空洞，至此，已经得到最终的水体信息，如果有必要需要显示面积较大的水体，可以再设置一个面积阈值，保留面积大于给定阈值的水体区域。