基于Google Earth遥感影像的城镇信息提取方法

摘要

本发明公开一种基于Google Earth遥感影像的城镇信息提取方法，首先对遥感影像进行中值滤波等预处理，再根据sobel边缘检测算子生成影像的边缘图像，然后构造高斯空间域加权边缘密度生成算子，生成边缘密度图像，再对边缘密度图像进行OSTU最优阈值分割，最后对分割图进行后处理，最终获取遥感影像中的城镇信息。本发明的提取方法效果好、速度快、准确率高，适用于城镇尤其是城郊结合处的城镇信息的提取。

说明书

技术领域

本发明涉及一种城镇信息提取方法，具体涉及一种基于Google Earth遥感影像的城镇信息提取方法。 

背景技术

城镇作为独立于其周边农用地、植被、裸地或者水体的地物，对于土地规划、土地覆盖、利用制图都有着重要意义。二十一世纪以来，随着空间技术的发展和各类图像处理、模式识别算法的出现，基于遥感手段的城镇信息提取算法越来越丰富，在各项制图工作和研究中取得了广泛应用。目前，常用的城镇信息提取算法都是基于多光谱遥感影像如SPOT和TM遥感影像，利用特定波段的光谱值生成归一化建筑物指数(NDBI)或其改进形式，再基于监督分类或者阈值分割方法进行城镇信息提取，但是鲜见将Google Earth遥感影像作为数据源或使用空间结构信息的提取方法。 

本发明针对Google Earth遥感影像中城镇区域和其他区域在光谱空间排列结构上的差异而设计了本算法，首先对遥感影像进行预处理，再通过sobel边缘检测算子生成边缘图像，然后根据空间域高斯边缘密度生成算子生成边缘密度图像，最后通过0STU阈值分割方法完成城镇信息提取。 

关于城镇信息提取问题，鲜见基于空间结构化信息的方法，一般仅利用了光谱信息，未能考虑到空间地物丰富的结构信息，故未能够最大化利用遥感影像中丰富的信息；其次，Google Earth遥感影像获取容易，但针对这一数据源开展的研究并不多，因此有必要针对该数据源做相应的研究工作。 

综上，提出基于Google Earth遥感影像结构信息的城镇信息提取算法非常有意义，也将为相近研究提供一种新思路。 

发明内容

发明目的：本发明的目的是针对现有城镇信息提取技术中存在未利用空间信息和Google Earth遥感影像在该领域利用较少的情况，提供一种基于GoogleEarth遥感影像的城镇信息提取方法。 

技术方案：本发明是一种基于Google Earth遥感影像的城镇信息提取方法， 包括以下步骤： 

(1)将原始RGB图像转化为灰度图像，利用中值滤波算法进行图像去噪，中值滤波模板大小取为3×3，以每个像元为中心，取其3×3邻域内的9个像元灰度值的中值作为该像元的滤波结果； 

(2)选择sobel边缘检测算子，使用横向和纵向两个模板对灰度图像进行卷积得到横向和纵向两个方向的梯度图像，累加两幅图像得到原始灰度图的近似梯度图像，通过下述公式梯度阈值g_thrd：g_thrd=k*mean(im_gradient)，其中im_gradient为梯度图像，mean(*)为梯度均值运算操作，k为阈值乘系数，乘系数一般取值范围为[0.5，2]，例如可以选取k=1，大于该阈值的认为是图像边缘，以此生成边缘特征图像； 

(3)在上述步骤(2)中所得的边缘特征图像中，根据空间地物在空间范围内的相关性，越靠近某像元的边缘对该像元的边缘密度值贡献越大，由于二维高斯函数的旋转对称不变形和单值性，因此二维高斯函数可以反应此相关性，生成高斯空间域加权边缘密度生成算子B_m×m，其中B_m×m中元素B_ij表示如下： $>B_{\mathrm{ij}}=\frac{\exp (-\sqrt{{(i-\frac{m+1}{2})}^2+{(j-\frac{m+1}{2})}^2}/{2\sigma }^2)}{\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{m}{\Sigma }}\exp (-\sqrt{{(i-\frac{m+1}{2})}^2+{(j-\frac{m+1}{2})}^2}/{2\sigma }^2)},$>根据上述算子B_m×m生成边缘密度图像，B_m×m的行数以及列数均为m； 

(4)基于步骤(3)中的边缘密度图像，通过OSTU方法选定阈值，对边缘密度图像进行阈值分割，选取大于给定面积阈值的图像斑块，最后进行闭运算以获取最终城镇信息。 

进一步的，所述步骤(3)中边缘密度生成算子的生成步骤如下： 

1)确定边缘密度生成算子B的大小m和参数σ²：m代表算子空间作用范围，σ²用于度量算子中各元素随着该元素至中心像元距离递增时的递减能力，σ²越小则算子的递减能力越强，越能体现中心像元的重要性，σ²越大则算子的递减能力越弱，越能体现邻域像元的重要性，σ²的选择能体现中心像元和其邻域像 元在空间范围内的重要程度，m和σ²可以通过多次试验根据最终分割结果与目视解译结果进行比较择优确定，例如，m=7，σ²=64其中 

$>B=\left(\begin{array}{ccccc}{b}_{11}& ·& ·& ·& b_{1m}\\ b_{21}& ·& ·& ·& b_{2m}\\ ·& & & & ·\\ ·& & & & ·\\ ·& & & & ·\\ b_{\mathrm{ml}}& ·& ·& ·& b_{\mathrm{mm}}\end{array}\right)(i=\mathrm{1,2},···,m;j=\mathrm{1,2},···,m)$>

上述i和j分别代表算子B中各元素对应的行号和列号； 

2)计算归一化参数c， 

$>c=\frac{1}{\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{m}{\Sigma }}\exp (-\sqrt{{(i-\frac{m+1}{2})}^2+{(j-\frac{m+1}{2})}^2}/{2\sigma }^2)};$>

3)根据归一化参数计算各个元素b(i，j)，组成边缘密度生成算子B， 

$>b(i,j)=\mathrm{cexp}(-\sqrt{{(i-\frac{m+1}{2})}^2+{(j-\frac{m+1}{2})}^2}/{2\sigma }^2);$>

4)生成大小为M×N的初始边缘密度图像D_M×N，计算各像元的边缘密度值， 

$>d(\mathrm{col},\mathrm{row})=\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{m}{\Sigma }}b(i,j)\mathrm{im}(\mathrm{row}-\frac{m-1}{2}+i-1,\mathrm{col}-\frac{m-1}{2}+j-1)$>

其中，im为原始遥感影像灰度图，col代表像元在原始遥感影像灰度图中的列数，row代表像元在原始遥感影像灰度图中的行数，所得边缘密度图像为： 

$>D=\left(\begin{array}{ccccc}{d}_{11}& ·& ·& ·& d_{1N}\\ d_{21}& ·& ·& ·& d_{2N}\\ ·& & & & ·\\ ·& & & & ·\\ ·& & & & ·\\ d_{M1}& ·& ·& ·& d_{\mathrm{MN}}\end{array}\right)(\mathrm{row}=\frac{m-1}{2}+\mathrm{1,2},···,M-\frac{m-1}{2},\mathrm{col}=\frac{m-1}{2}+\mathrm{1,2},···,N-\frac{m-1}{2})$>

因为边缘密度生成算子大小为m×m，因此对于图像靠近边界的行或列是未进行边缘密度计算的，对于每一个未参与计算的像元，在参与计算的像元中，寻找和该像元行列坐标距离最为接近的像元，并将此像元的边缘密度值赋给该像元，以此填补所有未计算像元的边缘密度值，最终得到的边缘密度图像为： 

$>D=\left(\begin{array}{ccccc}{d}_{11}& ·& ·& ·& d_{1N}\\ d_{21}& ·& ·& ·& d_{2N}\\ ·& & & & ·\\ ·& & & & ·\\ ·& & & & ·\\ d_{M1}& ·& ·& ·& d_{\mathrm{MN}}\end{array}\right)(\mathrm{row}=\mathrm{1,2},···,M,\mathrm{col}=\mathrm{1,2},···,N)$>

5)遍历边缘密度图像中各像元，使用冒泡法搜索边缘密度图像中各边缘密度值的最大值和最小值，首先初始化边缘密度最大值和最小值：d_max=0，d_min=100000；然后依次将i行j列像元的边缘密度值d_ij分别与d_max及d_min进行比较，若d_max<d_ij，则将d_ij赋给d_max，若d_min＞d_ij，则将d_ij赋给d_min，以此获取边缘密度最大值和最小值，对各像元边缘密度值线性拉伸至0-255，并依据四舍五入法取整： 

d_ij=round(255((d_ij-d_min)/(d_max-d_min))) 

上述公式中，d_ij为第i行第j列像元的边缘密度值，d_max为边缘密度最大值，d_min为边缘密度最小值，round(*)为四舍五入取整算子。 

进一步的，所述步骤(4)具体包括以下步骤： 

i)将边缘密度图像中的各个像元的边缘密度值按照一定间隔分为若干等级，即边缘密度图像的灰度级，然后进行排列，可以使用集合P={p₁，p₂，…，p_L}来表示，其中L为边缘密度图像的灰度等级，由于边缘密度图像已经进行线性拉伸至[0，255]，则本发明中，L=256，边缘密度值为p_i的像元数量为r_i，则总的像元数量为其中若L<256，则相当于对边缘密度图像进行灰度级压缩； 

ii)以第k个边缘密度值p_k为分界处，将集合p分为两类，分别为P₀={p₁，p₂，…，p_k}和P₁={p_k+1，p_k+2，…，p_L}，(k=1，2，…，L)，依次计算两个类别之间的类间方差，类间方差计算公式如下所示： 

$>{\sigma }_R^2\left(k\right)=\frac{{[{\mu }_T\omega \left(k\right)-\mu \left(k\right)]}^2}{\omega \left(k\right)[1-\omega \left(k\right)]}$>

其中，$>{\mu }_T=\underset{i=1}{\overset{L}{\Sigma }}{\mathrm{ip}}_i,\mu \left(k\right)=\underset{i=1}{\overset{k}{\Sigma }}{\mathrm{ip}}_i,\omega \left(k\right)=\underset{i=1}{\overset{k}{\Sigma }}{p}_i,$>P_i=r_i/R； 

iii)当类间方差取得最大值时：则此时的k对应着最优分割阈值，即为自适应最优分割阈值。 

有益效果：本发明提出一种基于Google Earth遥感影像的城镇信息提取算方法，能够方便快速地提取出城镇信息，用于大尺度城市制图和城市扩张研究等。本发明可以充分利用遥感影像的空间结构信息，本发明引入的空间域高斯加权边缘密度加权算子充分考虑了地物之间的空间相关性，使得各地物的边缘密度生成更为合理，最后通过OSTU阈值分割方法方便快速地提取城镇信息。 

附图说明

图1是本发明实施例中采集的城镇信息原始遥感影像的灰度图像； 

图2是图1转化为灰度图的梯度图像； 

图3是由梯度图像转化得到的边缘图像； 

图4是图3的边缘密度图像； 

图5是图4的OSTU分割结果示意图； 

图6是图5中的OSTU分割结果叠加于边缘图像获得的结果的示意图； 

图7是本发明的城镇信息提取流程图。 

具体实施方式

下面对本发明技术方案结合附图进行详细说明。 

本发明的一种基于Google Earth遥感影像的城镇信息提取方法，包括以下步骤： 

(1)将原图像转化为灰度图像，利用中值滤波算法去除原始遥感影像中存在的噪声； 

(2)使用sobel边缘检测算子生成边缘特征图像； 

(3)根据空间域高斯函数，生成一定大小的边缘密度生成算子，然后生成边缘密度图像； 

(4)基于步骤(3)中的边缘密度图像，通过OSTU方法选定阈值，对边缘密度图像进行阈值分割，选取大于给定面积阈值的图像斑块，最后进行闭运算以获取最终城镇信息。 

进一步的，所述步骤(3)中边缘密度生成算子的生成步骤如下： 

1)确定边缘密度生成算子B的大小m和参数σ²，其中 

$>B=\left(\begin{array}{ccccc}{b}_{11}& ·& ·& ·& b_{1m}\\ b_{21}& ·& ·& ·& b_{2m}\\ ·& & & & ·\\ ·& & & & ·\\ ·& & & & ·\\ b_{m1}& ·& ·& ·& b_{\mathrm{mm}}\end{array}\right)(i=\mathrm{1,2},···,m;j=\mathrm{1,2},···,m)$>

2)计算归一化参数c， 

$>c=\frac{1}{\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{m}{\Sigma }}\exp (-\sqrt{{(i-\frac{m+1}{2})}^2+{(j-\frac{m+1}{2})}^2}/{2\sigma }^2)};$>

3)根据归一化参数计算各个元素b(i，j)，组成边缘密度生成算子B， 

$>b(i,j)=\mathrm{cexp}(-\sqrt{{(i-\frac{m+1}{2})}^2+{(j-\frac{m+1}{2})}^2}/{2\sigma }^2);$>

4)生成大小为M×N的初始边缘密度图像D_M×N，计算各像元的边缘密度值， 

$>d(\mathrm{row},\mathrm{col})=\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{m}{\Sigma }}b(i,j)\mathrm{im}(\mathrm{row}-\frac{m-1}{2}+i-1,\mathrm{col}-\frac{m-1}{2}+j-1)$>

其中，im为原始遥感影像灰度图，col代表像元在原始遥感影像灰度图中的列数，row代表像元在原始遥感影像灰度图中的行数，所得边缘密度图像为： 

$>D=\left(\begin{array}{ccccc}{d}_{11}& ·& ·& ·& d_{1N}\\ d_{21}& ·& ·& ·& d_{2N}\\ ·& & & & ·\\ ·& & & & ·\\ ·& & & & ·\\ d_{M1}& ·& ·& ·& d_{\mathrm{MN}}\end{array}\right)(\mathrm{row}=\frac{m-1}{2}+\mathrm{1,2},···,M-\frac{m-1}{2},\mathrm{col}=\frac{m-1}{2}+\mathrm{1,2},···,N-\frac{m-1}{2})$>

因为边缘密度生成算子大小为m×m，因此对于图像靠近边界的行或列是未进行边缘密度计算的，对于每一个未参与计算的像元，在参与计算的像元中，寻找和该像元行列坐标距离最为接近的像元，并将此像元的边缘密度值赋给该像元，以此填补所有未计算像元的边缘密度值，最终得到的边缘密度图像为： 

$>D=\left(\begin{array}{ccccc}{d}_{11}& ·& ·& ·& d_{1N}\\ d_{21}& ·& ·& ·& d_{2N}\\ ·& & & & ·\\ ·& & & & ·\\ ·& & & & ·\\ d_{M1}& ·& ·& ·& d_{\mathrm{MN}}\end{array}\right)(\mathrm{row}=\mathrm{1,2},···,M,\mathrm{col}=\mathrm{1,2},···,N)$>

5)遍历所有像素点，搜索边缘密度图像的最大值和最小值，对各个像元进行线性拉伸至0-255，并依据四舍五入法取整： 

d_ij=round(255((d_ij-d_min)/(d_max-d_min))) 

进一步的，所述步骤(4)具体包括以下步骤： 

i)将边缘密度图像中的各个像元的边缘密度值按照一定间隔(例如：间隔为1)分为若干等级，即边缘密度图像的灰度级，然后进行排列，可以使用集合P={p₁，p₂，…，p_L}来表示，其中L为边缘密度图像的灰度等级，由于边缘密度图像已经进行线性拉伸至[0，255]，本实施例中，L=256，边缘密度值为p_i的像元数量为r_i，则总的像元数量为

ii)以第k个边缘密度值p_k为分界处，将集合p分为两类，分别为P₀={p₁，p₂，…，p_k}和P₁={p_k+1，p_k+2，…，p_L}，(k=l，2，…，L)，依次计算两个类别之间的类间方差类间方差计算公式如下所示： 

$>{\sigma }_R^2\left(k\right)=\frac{{[{\mu }_T\omega \left(k\right)-\mu \left(k\right)]}^2}{\omega \left(k\right)[1-\omega \left(k\right)]}$>

其中，$>{\mu }_T=\underset{i=1}{\overset{L}{\Sigma }}{\mathrm{ip}}_i,\mu \left(k\right)=\underset{i=1}{\overset{k}{\Sigma }}{\mathrm{ip}}_i,\omega \left(k\right)=\underset{i=1}{\overset{k}{\Sigma }}{p}_i,$>P_i=r_i/R。 

iii)当类间方差取得最大值时：则此时的k对应着最优分割阈值，即为自适应最优分割阈值。 

如上所述，尽管参照特定的优选实施例已经表示和表述了本发明，但其不得解释为对本发明自身的限制。在不脱离所附权利要求定义的本发明的精神和范围前提下，可对其在形式上和细节上作出各种变化。