由合成孔径雷达图像中检测建筑和提取数目信息的方法

摘要

本发明公开了一种由合成孔径雷达图像中检测建筑和提取数目信息的方法，该方法通过高低阈值、形态学开闭操作以及去除小面积区域，保证了建筑检测结果有较小虚警率，同时最大限度的连接了来自同一建筑的相邻亮特征，而且，该方法可以提高区域内建筑数目信息的提取精度并能很好的分离非常接近甚至相互重叠的来自不同建筑的亮特征。而且，本方法采用单幅高分辨率SAR数据，其是最容易获取的数据形式，相比其他形式的数据可以节约数据成本。本发明中的各项子算法非常简单和高效（主要为阈值运算和数学形态学运算），这使得它可以方便应用于大数据量和大面积区域情况。

说明书

技术领域

本发明涉及遥感信息提取与应用领域，更具体地涉及一种从高分 辨率合成孔径雷达图像中检测建筑和提取建筑数目信息的方法。

背景技术

合成孔径雷达（Synthetic Aperture Radar,SAR）以其高分辨率和 全天候，全天时，大面积的数据获取能力已成为世界各国普遍重视的 对地观测遥感技术。先进的星载雷达系统，如德国雷达卫星 TerraSAR-X，加拿大雷达卫星RADARSAT-2和意大利的 COSMO/SkyMet，已经可以提供米级分辨率的图像数据，而先进的机 载SAR系统已经能够达到分米级水平。在这样的分辨率下，城市建 筑的几何和细节信息已经清晰可见，城市建筑信息的提取完全可以基 于高分辨率SAR数据源进行。一个或多个空间区域内的建筑数目以 及此数目的变化，在很多方面如城市规划，土地管理，环境保护资源 分配，居住区环境评估及相关政策制定具有重要指导性作用。

建筑检测是进一步提取建筑数目信息的关键步骤。理想情况下， 如果检测结果中每个独立检出区域对应一栋独立建筑，则指定区域内 的建筑数目信息可以通过简单的对检出区域计数得到。由于SAR特 殊的侧视相干成像原理，单个建筑在高分辨率SAR图像上常表现为 高亮的点簇或者线簇（由掩叠和二次散射现象造成）特征。这个特征 明显区别于大部分自然地物，是目前SAR图像建筑检测的主要基础。 为了得到一个区域内建筑的准确数目，建筑检测需要连接同一建筑的 离散亮特征，同时分离非常接近甚至相互重叠的来自不同建筑的亮特 征。

到目前为止，人们进行了许多关于高分辨率SAR建筑检测方面 的研究工作。

例如，Thiele等基于多视角SAR数据利用边缘检测算子提取建 筑边界，然后结合InSAR得到高度信息判断建筑是否存在。Sportouche 等基于高分辨率光学图像和SAR数据,利用边缘检测算子提取亮的线 状特征，利用阈值检测建筑阴影，将满足固定位置和近邻关系的亮特 征与阴影组合判断为建筑候选，最后融合光学图像提取的信息进行综 合判断和三维信息提取。Simeneto等基于立体像对SAR数据，只针 对宽大的方形建筑，利用像素分类以及边缘检测得到SAR图像建筑 亮特征的候选区域，然后利用houph变换检测L型，T型以及X型 的亮特征，最后融合像对的不同图像上的亮特征检测结果来进一步减 小建筑检测的虚警率。Xu等基于多视SAR数据使用Wishart分布描 述高分辨率SAR中的建筑的统计分布，然后基于恒虚警边缘检测和 houph变换检测建筑的平行结构以定位建筑。Michaelsen等基于小于 1m空间分辨率精度的SAR影像，首先检测明显的亮点和线特征，然 后利用GESTALT系统进行分组和优化，最终检测建筑。Wegner等 基于InSAR数据和光学正射影像，在条件随机场（Conditional Random  Field,CRF）框架下利用建筑的光学图像特征和SAR数据特征（主要 是由二次散射造成的高亮线特征）进行建筑和非建筑的分类。

这些工作对于建筑数目信息提取来说存在以下几个问题：

（1）这些工作都是以进一步提取建筑的三维几何信息为目的。 为了能够提取出精确建筑几何信息，针对尽可能高分辨率以及尽可 能多的数据，例如空间分辨率在1m及以下的SAR数据，多视角 SAR数据，SAR立体像对数据，InSAR数据，光学与SAR数据， 光学与InSAR数据。对于建筑数目信息提取来说，如此苛刻的数据 要求并不必要，反而会造成无谓的成本浪费。

（2）只考虑了各个建筑间相对离散的区域，当应用于建筑密集 区域时，建筑数目信息的精度会迅速下降。换句话说，即没有考虑 如何处理接近甚至相互重叠的来自不同建筑的亮特征情况。

（3）多数工作只检测特殊类型的建筑，例如宽大的长方形建筑， 具有L形状亮特征的建筑。

因此，现有的由合成孔径雷达图像中检测建筑和提取数目信息 的方法还是不理想。

发明内容

本发明目的是针对现有高分辨率SAR图像建筑检测方法的不 足，提供一种由合成孔径雷达图像中检测建筑和提取建筑数目信息的 方法，这种方法不仅考虑了各个建筑间相对离散的区域，还考虑了如 何处理接近甚至相互重叠的来自不同建筑的亮特征情况，并且这种方 法可以检测多种类型的建筑，算法也简单。

为实现本发明的目的，本发明的由合成孔径雷达图像中检测建筑 和提取数目信息的方法，包括通过合成孔径雷达传感器得到高分辨率 图像并利用计算机进行处理，其中孔径雷达图像又称SAR图像，在 使用计算机进行处理时，可以在所述SAR图像的正常像素取值范围 内选取特定取值作为阈值，其中，该方法包括以下具体步骤：

步骤1，通过安装在计算机中的高分辨率SAR图像处理程序打 开高分辨率SAR图像；

步骤2，在所述SAR图像的像素正常取值范围内选择高低两个 阈值，基于选择的高低两个阈值，利用计算机分别对打开的SAR图 像进行二值化处理，设定建筑区域为前景；

步骤3，将步骤2中得到的低阈值二值化结果中的前景区域设定 为起点，在高阈值二值化结果中的前景区域内扩展范围，使两个结果 中位置相连接或者重叠的前景区域合并为新的前景区域，去除没有高 阈值二值化结果中的前景区域和与之有重叠或者相连接的低阈值二 值化结果中的前景区域；

步骤4，合并步骤3结果中相邻很近的前景区域，关闭小的间隙；

步骤5，分离步骤4结果中微弱相连的区域，消除小的离散点和 尖峰；

步骤6，去除步骤5结果中面积较小的前景区域；

步骤7，设定五个以上的代表真实建筑区域的前景区域作为建筑 区域样本；

步骤8，利用计算机计算选取的建筑区域样本的图像特征，并获 取该建筑区域内的真实建筑数目；

步骤9，选择下一样本，重复步骤8，直到完成所有样本的计算；

步骤10，利用步骤7至步骤9所计算出的数据集，对样本的图 像特征与真实建筑数目进行相关分析，筛选出与真实建筑数目相关性 较高的图像特征集；

步骤11，以真实建筑数目为因变量，样本图像特征子集为自变 量做多变量线性回归分析，构建线性回归方程；

步骤12，利用公式在计算机中评估线性归回方程的有效性；

步骤13，选择不同的样本图像特征子集，重复步骤11至步骤12， 直至得到最佳性能的线性回归方程；

步骤14，利用线性回归方程初步确定步骤6结果中每个前景区 域内的建筑数目；

步骤15，将步骤6结果中的各前景区域切分为线性回归方程所 估计数目的独立区域，得到最终建筑检测结果；

步骤16，利用步骤15的建筑检测结果，以指定区域内独立前景 区域的个数作为该区域内建筑数目，通过计算机计算出指定区域内的 建筑数目信息。

进一步的，上述的由合成孔径雷达图像中检测建筑和提取数目信 息的方法，其中，所述步骤16还可以采用以下方法：利用步骤6建 筑检测结果和步骤13得到的线性回归方程通过计算机计算指定区域 内的建筑数目信息。具体如下，假设在指定区域内包含m个前景区 域，则整个区域内建筑数目y由下式算出：

$y={\Sigma }_{i=1}^m{y}_i$

y_i=β₀+β₁x_1,i+...+β_jx_j,i+....+β_kx_k,i

其中，x_j，i为由第i个前景区域内算出的第j个图像特征， β_i,i=0，1，...,k为回归方程系数，y_i为利用回归方程计算出的第i个前 景区域内的建筑数目。

更进一步的，上述的由合成孔径雷达图像中检测建筑和提取数目 信息的方法，其中，所述步骤15中将步骤6结果中的各前景区域切 分为线性回归方程所估计数目的独立小区域包括以下子步骤：

所述步骤15中将步骤6结果中的各前景区域切分为线性回归方 程所估计数目的独立小区域包括以下子步骤：

步骤15.1，利用安装在计算机中的SAR图像处理程序复制一个 前景区域到一个包含该区域的小矩形内，构成一幅小的二值图像（前 景区域与背景区域取不同值）；对前景区域反复进行形态学腐蚀操作， 标记极限腐蚀区域和背景区域；

步骤15.2，计算二值图像的距离变换图；

步骤15.3，计算所有标记区域的显著性指标；

步骤15.4，利用计算机软件将各标记区域按照显著性指标从高到 低排序，并按次序选择前N（N为回归方程计算出的该区域内建筑数 目）个标记区域；

步骤15.5，以步骤15.4选择的标记区域为初始溢流区域，在步 骤15.2得到的距离变换图上运行流域变换得到对前景区域的分割结 果；

步骤15.6，以步骤15.5得到的分割结果代替步骤6结果中对应 的前景区域；

步骤15.7，重复15.1至15.6步骤，直至所有步骤6结果中的前 景区域都被处理完成。

更进一步的，上述的由合成孔径雷达图像中检测建筑和提取数目 信息的方法，其中，所述SAR图像处理程序是指可以对SAR图像进 行输入输出或对其中的值进行改变并提取信息操作的计算机软件。

再更进一步的，上述的由合成孔径雷达图像中检测建筑和提取数 目信息的方法，其中，所述SAR图像处理程序是基于IDL语言编写 的软件。

本发明具有以下优点和积极效果：

（1）通过高低阈值、形态学开闭操作以及去除小面积区域，保 证了建筑检测结果有较小虚警率，同时最大限度的连接了来自同一建 筑的相邻亮特征；

（2）利用回归模型建模建筑区域特征与区域内实际建筑栋数间 的复杂关系，提高了区域内建筑数目信息的提取精度；

（3）利用形态学流域变换将大检测出区域切分为符合回归模型 估计数目的小区域，很好的分离了非常接近甚至相互重叠的来自不同 建筑的亮特征。这个结果可进一步用于其他建筑信息（如建筑位置， 形状，高度等信息）的提取工作；

（4）针对单幅高分辨率SAR数据，其是目前最容易获取的数据 形式，相比其他形式的数据可以节约数据成本；

（5）本发明中的各项子算法非常简单和高效（主要为阈值运算 和数学形态学运算）。这使得它可以方便应用于大数据量和大面积区 域情况。

附图说明

图1是本发明的分辨率为3米的SAR城区图像；

图2是本发明的步骤6去除面积较小的前景区域得到的建筑检测 结果图像；

图3是本发明步骤15将一个大的前景区域切分4个独立区域的 过程图像；

图4是本发明前景区域对应的光学图像；

图5是本发明步骤6检测出的对应于图4场景的前景区域图像；

图6是本发明步骤15将图5进行切分的结果图像。

具体实施方式

下面结合附图和实施例对本发明作进一步详细说明。

先看图1、图2，图1是本发明分辨率为3米的SAR城区图像， 其中高亮的区域多数为建筑，图2是步骤6中去除步骤5结果中面积 较小的前景区域得到的结果。再看图3，图3示意了步骤15将一个 大的前景区域切分4个（线性回归方程所估计出的数目）独立区域的 过程，该过程包括a、b、c、d、e五个步骤，其中a为包含多个建筑 的前景区域；b为距离变换图像；c为标记出的极限腐蚀区域和背景 区域（白色部分为标记）；d为基于显著性指标选出的前4个最显著 的独立标记区域和背景区域；e为基于标记图像d在距离变换图像b 上运行流域变换得到的分割结果。再看图4、图5、图6，图4是本 发明前景区域对应的光学图像；图5是本发明步骤6检测出的对应于 图4场景的前景区域；图6是本发明步骤15将图5进行切分的结果。

本发明所提供检测建筑和提取建筑数目信息方法的具体实施方 式，按步骤详述如下：

步骤1，通过安装在计算机中的高分辨率SAR图像处理程序打开 高分辨率SAR图像；其中高分辨率SAR图像是使用合成孔径雷达 （Synthetic Aperture Rader）传感器得到的，传感器可以安装在飞机平 台上，卫星平台上或者飞艇平台上等等，在使用计算机对SAR图像进 行处理时，可以在所述SAR图像的正常像素取值范围内选取特定取值 作为阈值。

步骤2，在所述SAR图像的像素正常取值范围内选择高低两个阈 值，基于选择的高低两个阈值，利用计算机分别对打开的SAR图像进 行二值化处理，设定建筑区域为前景。

在这里，所述二值化处理过程如下：

设一个图像像素的取值为x，这个x有一个取值范围【xmin,xmax】。 阈值就是在x取值范围内选择的一个具体取值。利用一个阈值可以将 一幅像素取值范围为【xmin,xmax】的图像变为像素取值范围为{0,1} 的二值图像，具体是这样实现：将每个图像像素的具体取值与这个阈 值进行比较，如果像素取值大于阈值则将像素取值赋为1（或0），反 之则赋为0（或1）。

具体实施时，高低阈值可以根据经验或者根据样本数据训练得 到，也可以通过直接在图像上交互取阈值的方式取得。直接交互取得 阈值的方式更加符合实际应用的大多数情况。在这里，较高的阈值检 测出部分可信度高的建筑亮区，检测结果具有很小的虚警；较低的阈 值检测出绝大部分建筑亮区，但是检测结果中包含了大量虚警。

步骤3，将步骤2中得到的低阈值二值化结果中的前景区域设定 为起点，在高阈值二值化结果中的前景区域内扩展范围，使两个结果 中位置相连接或者重叠的前景区域合并为新的前景区域，去除没有高 阈值二值化结果中的前景区域和与之有重叠或者相连接的低阈值二 值化结果中的前景区域。

具体实施时，这个步骤可以采用形态学重建方法（以高阈值二值 化结果为标记图，低阈值二值化结果为掩模图）完成。形态学重建方 法为图像领域公知的技术，在此不予赘述。

步骤4，合并步骤3结果中相邻很近的前景区域，关闭小的间隙。

具体实施时，可利用形态学闭算子，根据需求选择结构基元的种 类和大小。例如采用大小为5*5的矩形结构单元进行形态学闭操作， 连接距离2个像素以内的前景区域，填充小于此结构基元的间隙与孔 洞。

步骤5，分离步骤4结果中微弱相连的区域，消除小的离散点和尖 峰。

具体实施时，可利用形态学开算子，根据需求选择结构基元的种 类和大小。例如采用大小为3*3的矩形结构单元进行形态学开操作， 分离仅有单个像素相连接的相邻前景区域，去除小于此结构基元的离 散点与尖峰。

步骤6，去除步骤5结果中面积较小的前景区域。

由于SAR相干成像原理以及地物情况的复杂性，裸地以及草地 林地等地物也会产生部分亮点，这些亮点相对建筑亮区来说相对孤立 并且面积较小。因此通过去除小于某个面积阈值的前景区域可以去除 大量的虚警情况。

步骤7，设定五个以上的代表真实建筑区域的前景区域作为建筑 区域样本。

具体实施时，建筑区域样本的数量通常取较多的数量。所需样本 的最小数目是与步骤11回归方程参数数目相关的，即如果回归方程 的参数数目为10，则样本数目最少为10个，如果样本数目少于10， 那么回归方程的参数就难以有效估计了。实际上，为了使回归方程参 数估计的更准确，建筑区域样本的数量要远大于回归方程的的参数数 目。

步骤8，利用计算机计算选取的建筑区域样本的图像特征，并获 取该建筑区域内的真实建筑数目。

具体实施时，可以首先将样本区域轮廓映射到对应地区的高分辨 率光学影像上，在光学影像上辨认和计数得到该区域内真实建筑数 目。区域样本的图像特征具有非常多的种类，本发明优先使用面积， 致密度，区域轮廓拟合椭圆的长短轴长度等特征。

步骤9，选择下一样本，重复步骤8，直到完成所有样本的计算。

步骤10，利用步骤7至步骤9所计算出的数据集，对样本的图像特 征与真实建筑数目进行相关分析，筛选出与真实建筑数目相关性较高 的图像特征集。

为保证后面步骤构建回归方程的有效性，尽量筛选出与真实建筑 数目具有较高相关性的图像特征。例如可以采用常见的皮尔逊积差相 关系数（即两个变量的协方差除以两个变量的标准差）进行相关性的 评价。

步骤11，以真实建筑数目为因变量，样本图像特征子集为自变量 做多变量线性回归分析，构建线性回归方程。

以y表示真实建筑数目，x_i表示第i个图像特征，则k元线性 回归方程可表示为：

y=β₀+β₁x₁+...+β_ix_i+....+β_kx_k                 公式（1）

其中β₀,β₁,β₂,..,β_i,..,β_k为回归系数，可使用最小二乘方法基于样 本数据集{y,x₁，...，x_i，...，x_k}_{j=1，...，N}估计得到。这个样本数据集是由步骤7至 步骤10计算出的每个建筑样本区域的真实建筑数目和图像特征数据 构成。

步骤12，利用公式在计算机中评估线性归回方程的有效性。

在使用线性回归方程之前，一个重要的步骤是检验模型的拟合优 度以及所估计模型参数和方程的统计显著性。检验模型的拟合优度可 以使用修正的R平方系数：

${\bar{R}}^2=1\frac{\Sigma {(y_i-{\hat{y}}_i)}^2/(n-k-1)}{\Sigma (y_i-\bar{y})/(n-1)}=1-(1-R^2)\frac{n-1}{n-k-1}$公式（2）

这里为因变量y样本数据的总体均值，为回归方程针对第i 个自变量样本估计的因变量取值，y_i是第i个因变量样本取值，n为样 本数据的数目，k是自变量的数目。此系数越接近1则说明模型的拟合 优度越好。

检验所估计的回归方程各参数的统计显著性可采用t-检验，构建 如下统计量：

$T=\frac{b_i-{\beta }_i}{\mathrm{SE}\left(b_i\right)}~t_{\frac{\alpha }{2}}(n-k-1)$公式（3）

其中，b_i是对第i个回归方程参数β_i的估计值,SE(b_i)是b_i的样本 标准差,n为样本数据的数目，k是自变量的数目,α为显著性水平。t- 检验的目的是检测因变量是否被每个自变量所影响。

检测回归方程的统计显著性可采用F-检验，构建如下统计量：

$F=\frac{\Sigma ({\hat{y}}_i-\bar{y})/k}{\Sigma (y_i-{\hat{y}}_i)/(n-k-1)}~F_{\alpha }(k,n-k-1)$公式（4）

这里为因变量y样本数据的总体均值，为回归方程针对第 i个自变量样本估计的因变量取值，n为样本数据的数目，k是自变量 的数目,α为显著性水平。F-检验的目的是检测因变量与自变量的线 性关系。

步骤13，选择不同的样本图像特征子集，重复步骤11至步骤12， 直至得到最佳性能的线性回归方程。

步骤14，利用线性回归方程初步确定步骤6结果中每个前景区域 内的建筑数目。

具体实施时使用下面方程回归方程计算每个前景区域内的建筑 数目：

y=β₀+β₁x₁+...+β_ix_i+....+β_kx_k

这里β₀,β₁,β₂，...,β_i，...,β_k为步骤13得到的回归方程系数，x_i为当前前 景区域的第i个图像特征。

步骤15，将步骤6结果中的各前景区域切分为线性回归方程所估 计数目的独立区域，得到最终建筑检测结果。

步骤16，通过计算机计算指定区域内的建筑数目信息。

这里有两个方案:

方案一，利用步骤6建筑检测结果和步骤13得到的线性回归方程 计算指定区域内的建筑数目信息。具体如下，假设在指定区域内包含 m个前景区域，则整个区域内建筑数目y由下式算出：

$y={\Sigma }_{i=1}^m{y}_i$

y_i=β₀+β₁x_1，i+β₂x_2，i+β₃x_3，i

其中，x_j，i为由第i个前景区域内算出的第j个图像特征， β_i,i＝0，1,2,3为回归方程系数（步骤13中得到的最佳回归方程的系数）， y_i为利用回归方程计算出的第i个前景区域内的建筑数目。

方案二，利用步骤15的建筑检测结果，以指定区域内独立前景区 域的个数作为该区域内建筑数目。

其中方案一为优选方案。

所述步骤15中将步骤6结果中的各前景区域切分为线性回归方 程所估计数目的独立小区域包括以下子步骤：

步骤15.1，利用安装在计算机中的SAR图像处理程序复制一个前 景区域到一个包含该区域的小矩形内，构成一幅小的二值图像（前景 区域与背景区域取不同值）。对前景区域反复进行形态学腐蚀操作， 标记极限腐蚀区域和背景区域；

具体实施时，采用一个半径为r的结构基元B对前景区域A进行 腐蚀操作。由数学形态学理论，区域A被结构基元B腐蚀表示为 腐蚀具有收缩图像区域的作用，对于给定前景区 域A，使用结构基元B反复做腐蚀运算，不断剥掉厚度为r的一层， 随着腐蚀的不断作用，不连通的区域会不断产生，同时某些区域又会 逐渐消失。一个连通成分在消失前的最后一步，称为最终连通成分， 而所有最终连通成分的并，称为相对半径r的极限腐蚀区域。本步 骤要标记小矩形内的所有最终连通成分以及背景区域。

步骤15.2，计算二值图像的距离变换图；

本步骤中距离变换图采用如下方法得到，以I代表当前二值图 像，通过计算二值图像中每个像素到离它最近的背景像素间的距离得 到初始的距离变换图d(I)，本步骤最终使用的距离变换图像为L-d(I)， 其中L为一个常数值。如为了保证距离变换取值不为负值，L需要选择 大于等于max(d(I))的取值。计算距离所采用的距离测度可以选取常用 的欧氏距离，棋盘距离，城区距离以及马氏距离等。

步骤15.3，计算所有标记区域的显著性指标；

将步骤15.2中距离变换图看作是一个三维的地形（像素的x和y坐 标标志了其地理位置，像素的取值代表了地形的高度），则步骤15.1 中标记的不同区域可以看作是这个地形中的不同盆地。评价标记区域 的显著性可以有不同的方式。最简单的评价方式是根据区域所对应盆 地的盆底高度进行，即盆底高度越低则此盆地越显著。此种评价并没 有考虑各盆地间的相对关系，将会造成地势低地区的盆地显著性评价 值都很高，而地势高地区的盆地显著性评价值都很低。本发明中推荐 使用动力学指标来评价盆地的显著性。

盆地M的动力学定义：由盆地M的底出发到达任意更低盆地底的 所有路径中选出路径的最高点高度最低的一条路径，将此路径的最高 点高度减去盆地M的盆底高度，差值即为盆地M的动力学值。从定义 可以看出，动力学不是一个局部的显著性标准，可能与其周围很大范 围内盆地（不仅仅是相邻盆地）的性质有关。本发明推荐在溢流过程 中两个盆地的水流相交时计算盆地的动力学值。首先说明溢流过程， 假设地形的各个标记区域所对应盆地的最低点都是凿穿的并与同一 地下水源相连，地下水位不断上升直至淹没整个地形的过程称为溢流 过程。溢流过程中两个不同盆地的水相交时计算盆地动力学值的方法 如下：认为地势高的盆地被地势低的盆地吞并，同时计算出地势较高 盆地的动力学值，等于当前水位的高度h减去此盆地的盆底高度h(M), 即dy(M)=h-h(M)。这样全局最低盆地将得不到盆地动力学值，可以 令其为一个极大值。此外当两个地势相同的盆地的水相交时，以上算 法也无法判定应当给哪个盆地计算动力学值，可以任选一个盆地计算 其动力学值。

步骤15.4，利用计算机软件将各标记区域按照显著性指标从高到 低排序，并按次序选择前N（N为回归方程计算出的该区域内建筑数 目）个标记区域；

步骤15.5，以步骤15.4选择的标记区域为初始溢流区域，在步骤 15.2得到的距离变换图上运行流域变换得到对前景区域的分割结果；

将步骤15.2中距离变换图看作是一个三维的地形（像素的x和y坐 标标志了其地理位置，像素的取值代表了地形的高度），步骤15.4选 择的标记区域看作是这个地形中的不同盆地，并假设只有这些标记区 域是初始溢流区域（即所对应盆地的最低点是凿穿的并与同一地下水 源相连）。流域变换在针对这个三维地形进行的一个溢流过程基础上 进行，当地下水逐渐上涨时，水会逐渐淹没盆地并向山脊上升。由于 所有盆地的水是同步上升的，则不同盆地的水流交汇处恰好是山脊 点。在每一个交汇处筑起堤坝，阻止不同水流的交汇。这样不断提高 水位，直到水面上只有堤坝。可以看到堤坝是闭合的并且将整个地形 区域切分为了不同的小块区域。切分区域的数目等于最初标记区域的 数目。这些切分开的区域就是本步骤所得的分割结果；

步骤15.6，以步骤15.5得到的分割结果代替步骤6结果中对应的前 景区域；

步骤15.7，重复15.1至15.6步骤，直至所有步骤6结果中的前景区 域都被处理完成。

结合本发明的内容进一步提供以下实施例：

本实施例将本发明所提供技术方案应用在基于1.25m*1.25m空间 分辨率的TerraSAR-X强度数据提取湖北省武汉市部分区域内的建筑 数目信息。整个过程如下：

步骤（1）打开一幅大小为20000×20000像素的高分辨率SAR图像。 该图像是TerraSAR-X拍摄的中国武汉市某处，图像空间分辨率为 1.25m*1.25m（是由3m*3m空间分辨率上过采样得到）。

步骤（2），利用在图像上交互选取阈值和观察的方式确定高低两 个阈值分别为700和400，对SAR图像分别进行二值化处理，大于阈值 区域为前景区域；

步骤（3），利用形态学重构方法融合两个阈值化结果，即以低阈 值（400）二值化结果图为掩模图像，高阈值（700）二值化结果图为 标记图像，进行形态学膨胀重构。最终，两个结果中位置相连接或者 重叠的前景区域合并为新的前景区域，没有高阈值二值化结果中的前 景区域与之有重叠或者相连接的低阈值二值化结果中的前景区域被 去除；

步骤（4），采用5*5矩形结构基元对步骤（3）结果进行形态学闭 操作，合并相邻很近的前景区域，并关闭小的间隙。

步骤（5），采用3*3矩形结构基元对步骤（4）结果进行形态学开 操作，分离微弱相连的区域，并消除小的离散点和尖峰。

步骤（6），去除步骤（5）结果中面积小于25个像素的前景区域；

步骤（7），对照光学参考图，在步骤（6）结果中选择48个独立 前景区域，作为建筑区域样本；

步骤（8），计算选取的建筑区域样本的5个图像特征（面积，质 密度，拟合椭圆的短轴长度，拟合椭圆的长轴长度，周长）；将该样 本区域映射到光学参考图的对应位置，通过人工辨认和计数得到该区 域内的真实建筑数目；

步骤（9），选择下一样本，重复步骤（8），直到完成所有48个样 本的计算。将所有数据结果存为数据集{y,x₁，...，x_i，...，x₅}_{j＝1，...，48}，其中y为 建筑区域样本内的真实建筑栋数，x_i为第i个图像特征。

步骤（10），利用步骤（9）数据集，分别计算每个图像特征与真 实建筑数目的皮尔逊积差相关系数（即两个变量的协方差除以两个变 量的标准差）。发现3个特征（面积，致密度，拟合椭圆短轴长度）与 真实建筑数目间的相关系数绝对值大于0.8，其他特征与真实建筑数 目间的相关系数绝对值小于0.6。因此挑选特征“面积”，“致密度”，“拟 合椭圆短轴长度”作为与真实建筑数目相关性较高的图像特征集，样 本数据集更新为{y,x₁,x₂,x₃}_{j＝1，...，48}。

步骤（11），基于步骤（10）的样本数据集，构建以下3变量线性 回归方程：

y=β₀+β₁x₁+β₂x₂+β₃x₃

其中，y为某区域内建筑栋数，x₁,x₂,x₃分别为区域的“面积”，“致 密度”，“拟合椭圆短轴长度”特征。利用最小二乘算法估计回归方程 系数为β₀=0.4276424，β₁=0.0004744，β₂=0.0067597，β₃=0.0197668。

步骤（12），根据公式（4）计算F统计量等于71.47016，查F分 布表可得F_0.05(3,44)=2.82。由于71.47016远远大于F_0.05(3,44)所以可 以判定整体的线性关系有效且显著；根据公式（3）分别计算回归方 程3个系数的相关t统计量分别为t1=4.148，t2=3.022，t3=2.466, 查t分布表可得t_0.025(44)=2.01。由于各个系数的t统计量都大于 t_0.025(44)，则可判定各个自变量与因变量间具有显著的线性关系；根 据公式（2）计算修正的R平方系数为0.8297,比较接近1，说明构建 的自回归方程可以很大程度的拟合自变量与因变量间的关系。

步骤（13），接着分别构建y=β₀+β₁x₁+β₂x₂，y=β₀+β₂x₂+β₃x₃， y=β₀+β₁x₁+β₃x₃，y=β₀+β₁x₁，y=β₀+β₂x₂，y=β₀+β₃x₃六种情况下 的自回归方程，并且评价其性能。比较七种（加上前面步骤构建的一 种情况）不同情况下的自回归方程性能发现，步骤（11）所构建的回 归方程性能最佳，因此选择这个回归方程作后续应用。

步骤（14），对步骤（6）结果中每个独立的前景区域，计算其“面 积”，“致密度”，“拟合椭圆短轴长度”三个特征，代入步骤（13）得 到的线性回归方程得到该前景区域内的估计建筑数目；

步骤（15），对步骤（6）结果中每个独立的前景区域进行如下操 作：

步骤15.1），复制该前景区域到一个包含该区域的小矩形内，构 成一幅小的二值图像（前景区域取1背景区域取0）。利用3*3矩形结构 基元对前景区域反复进行形态学腐蚀操作，并标记极限腐蚀区域和背 景区域；

步骤15.2），计算步骤15.1）中二值图像的距离变换图，其中距离 测度采用欧式距离；

步骤15.3），计算步骤15.1）所有标记区域的盆地动力学指标；

步骤15.4），将各标记区域按照盆地动力学指标从高到低排序， 并按次序选择前N（N为回归方程在步骤（14）中计算出的该区域内 建筑数目的最邻近整数值）个标记区域；

步骤15.5），以步骤15.4）选择的标记区域为初始溢流区域，在步 骤15.2）得到的距离变换图上运行流域变换得到对前景区域的分割结 果；

步骤15.6），以步骤15.5）得到的分割结果代替步骤6结果中对应 的前景区域；

步骤（16），在打开的高分辨率SAR图像（TerraSAR-X拍摄的中 国武汉市某处）上圈选8片不规则形状的区域，采用两个不同方案计 算指定区域内的建筑数目信息。

方案一，利用步骤（6）建筑检测结果和步骤（13）得到的线性 回归方程计算指定区域内的建筑数目信息。具体如下，假设在指定区 域内包含m个前景区域，则整个区域内建筑数目y由下式算出：

$y={\Sigma }_{i=1}^m{y}_i$

y_i=β₀+β₁x_1，i+β₂x_2，i+β₃x_3，i

其中，x_j，i为由第i个前景区域内算出的第j个图像特征， β_i,i＝0，1,2,3为回归方程系数（步骤（13）中得到的最佳回归方程的系 数），y_i为利用回归方程计算出的第i个前景区域内的建筑数目。

方案二，利用步骤（15）的建筑检测结果，以指定区域内独立前 景区域的个数作为该区域内建筑数目。

其中方案一为优选方案。

表1显示了圈选的8片不规则区域内建筑数据信息提取的结果。采 用步骤（16）中方案一提取8片区域内建筑数目信息的绝对误差的平 均值为0.925。当采用步骤（16）中方案二提取区域内建筑数目信息 时，由于需要对回归方程估计出的建筑数目取整，绝对误差性能指标 有所不同，绝对误差的平均值为1.125。

表1：8片测试区域内建筑数目信息检测结果

当然，以上只是发明的具体应用范例，本发明还有其他的实施方式， 凡采用等同替换或等效变换形成的技术方案，均落在本发明所要求的 保护范围之内。