[发明专利]一种高分辨率遥感影像中的建筑物快速提取方法

权利要求书

1.一种高分辨率遥感影像中的建筑物快速提取方法，其特征在于，包括如下步骤：

步骤1)、读取遥感影像，提取灰度图像：

遥感图像中每个像素点的色彩具有R、G、B三个分量,P(r,g,b)代表一像素点,其中r代表该像素点的R分量的值，g代表该像素点的G分量的值，b代表该像素点的B分量的值；

灰度函数H(P)的值表示像素点P(r,g,b)的灰度值：

H(P)＝0.3r+0.59g+0.11b 式(1)

用式(1)按照从左到右、从上到下的顺序依次扫描遥感图像中的每个像素点，得到每个像素点的灰度值，按照式(1)扫描时同样的顺序组成第一灰度图像；

步骤2)、使用高斯滤波对第一灰度图像进行降噪处理：

2a)、计算大小为(2k+1)×(2k+1)的高斯模板U(x，y)：

其中，k为正整数且k≥1，x，y分别代表高斯模板中元素的横坐标和纵坐标；σ＞0，为第一灰度图像的平滑程度参数；

2b)、取k＝1，大小为3×3的高斯模板U(x，y)中9个元素的模板值分别为m₁、m₂、m₃、m₄、m₅、m₆、m₇、m₈、m₉；

第一灰度图像中，除左、上、右、下四条宽度为一像素的边缘外任一像素点P点的8邻像素，即紧邻像素点P的左上、上、右上、左、右、左下、下、右下8个方向上的像素；像素点P与8邻像素的灰度值按照第一灰度图像中从左至右、从上到下的顺序依次为g₁、g₂、g₃、g₄、g₅、g₆、g₇、g₈、g₉；

用该高斯模板U(x，y)按照从左到右、从上到下的顺序依次扫描第一灰度图像中除左、上、右、下四条宽度为一像素的边缘外的每一个像素P，各像素点P由高斯模板U(x，y)做高斯滤波的方法见式(3)；

所述的U₀表示以高斯模板U(m，n)做高斯滤波后所到的像素点P的灰度值；S为中间变量，表示像素点P和其8邻像素的灰度值的总和；

按照式(3)依从左到右、从上到下的顺序对第一灰度图像中像素点P进行依次扫描时，对当前所扫描的像素点P由高斯模板U(x，y)所确定的邻域内所有像素的加权平均灰度值替代当前像素点P的灰度值，扫描结束后得到第二灰度图像；

步骤3)、使用基于Otsu方法的自适应Canny算法提取图像中的边缘；

31)：用一阶偏导的有限差分来计算梯度的幅值和方向；

采用Canny算法,其中所采用的卷积模板为式(4)；

其中M₁、M₂分别为第二灰度图像x、y轴方向上的一阶偏导数矩阵,依照所述一阶偏导数矩阵M₁、M₂按照从左到右、从上到下的顺序依次扫描第二灰度图像中的每一个像素点，得到每个像素点x、y方向上的梯度幅值G_x、G_y：

G_x＝[f(x+1,y)-f(x,y)+f(x+1,y+1)-f(x,y+1)]/2 式(5)

G_y＝[f(x,y)-f(x,y+1)+f(x+1,y)-f(x+1,y+1)]/2 式(6)

x、y表示所处理像素点的横坐标、纵坐标值，f(x,y)表示第二灰度图像中坐标为(x,y)的像素点的灰度值，该像素点的梯度幅值G(x,y)和方向角θ(x,y)由式(7)和式(8)处理得到；

由式(7)和式(8)对第二灰度图像中每一个像素点按照从左到右、从上到下的顺序依次扫描，得到第二灰度图像中每一个像素点梯度幅值和方向角；

32)：非极大值抑制；

经过上述步骤31)的处理，得到了每个像素的梯度幅值G(x,y)和方向角θ(x,y)，为了得到最真实的边缘，须保留每个像素点在其梯度方向上的极大值，而删掉当前像素点梯度方向上的其他值，即非极大值抑制；对第二灰度图像中每个像素点按照从左到右、从上到下的顺序进行非极大值抑制，即处理得到第三灰度图像，第三灰度图像显示了遥感影像中建筑物影像的边缘；

非极大值抑制的具体处理方法如下；

判断是否是极大值，需要根据当前扫描像素点梯度方向所在的区间，在x轴和y轴方向分别进行插值计算，然后再进行比较,以找到最大值，梯度方向所在的区间指的是梯度方向α所在的角度范围，插值是为了更合理地判断当前像素点的梯度方向：

①当以及时，令G_p＝G₃w+G₂(1-w),G_n＝G₇w+G₈(1-w)；所述的w为中间变量；

②当以及时，令G_p＝G₇w+G₄(1-w),G_n＝G₃w+G₆(1-w)；

③当以及时，令G_p＝G₆w+G₉(1-w),G_n＝G₄w+G₁(1-w)；

④当以及时，令G_p＝G₂w+G₁(1-w),G_n＝G₈w+G₉(1-w)；

其中，G₅为当前处理像素的梯度幅值，G_k(1≤k≤4,6≤k≤9)表示8邻像素的梯度幅值，G_p表示梯度方向上前一个像素的梯度幅值，G_n表示梯度方向上后一个像素的梯度幅值；

当梯度方向α的值在上述相应的区间内时，当且仅当目前所处理像素的梯度幅值G₅大于前述①～④中各插值方法计算出的G_p和G_n时，G₅为极大值，也就表明当前处理像素是边缘点；

33)、计算最优高低阈值

读取第三灰度图像的宽度为W像素、高度为H像素，则第三灰度图像中的像素总数N为：

N＝W×H 式(9)

设第三灰度图像中每一像素的梯度幅值G(x,y)的范围为[0,G]，读取每一像素的梯度幅值G(x,y)，找梯度幅值等于确定值为i的像素，得到梯度幅值等于i的像素个数为N_i个，i∈[0,G]；设T₁为最优低阈值，T₂为最优高阈值；

非边缘点总个数N₁(T₁,T₂)为：

其中，N_i是梯度幅值等于i的像素的个数；

非边缘点所占比例值ω₁(T₁,T₂)为：

潜在边缘点总个数N₂(T₁,T₂)为：

潜在边缘点所占比例值ω₂(T₁,T₂)为：

真边缘点总个数N₃(T₁,T₂)为：

真边缘点所占比例值ω₃(T₁,T₂)为：

非边缘点的平均梯度μ₁(T₁,T₂)、潜在边缘点的平均梯度μ₂(T₁,T₂)和真边缘点的平均梯度μ₃(T₁,T₂)如下：

第三灰度图像总的梯度均值μ_s为：

μ_s＝μ₁(T₁,T₂)ω₁(T₁,T₂)+μ₂(T₁,T₂)ω₂(T₁,T₂)+μ₃(T₁,T₂)ω₃(T₁,T₂) 式(17)

非边缘点、潜在边缘点和真边缘点这三类像素的类间方差σ²(T₁,T₂)为：

σ²(T₁,T₂)＝ω₁(T₁,T₂)(μ₁(T₁,T₂)-μ_S)²+ω₂(T₁,T₂)(μ₂(T₁,T₂)-μ_S)²+ω₃(T₁,T₂)(μ₃(T₁,T₂)-μ_S)²

式(18)

使得σ²(T₁,T₂)的输入参数T₁、T₂在总范围[1,G-1]依次取值，当类间方差σ²(T₁,T₂)为最大值时，此时的T₁、T₂的值即为最优低阈值T₁和最优高阈值T₂；

34)：用双阈值算法检测和边缘连接；

根据上一步33)中得到的最优低阈值T₁和最优高阈值T₂对第三灰度图像进行过滤，像素点的梯度幅值如果大于高阈值则认为当前像素点必然是边缘点，即真边缘点，如果当前处理像素点的梯度幅值小于低阈值，则认为当前像素点必然不是边缘点；处于低阈值和高阈值之间的像素点是潜在的边缘点，即假边缘；将第三灰度图像根据高阈值过滤得到的边缘点链接成轮廓，当到达轮廓的端点时，在端点的8邻点寻找潜在的边缘点，再根据所找到的潜在的边缘点,在其8邻点中循环处理收集新的边缘点，如此处理整个第三灰度图像，直至边缘闭合，得到第四灰度图像；

步骤4)、基于Freeman链码的直线提取方法：

单条直线L与Y轴的夹角θ范围为0～π,共有八种情况的直线L：

4)-①θ＝0或θ＝π；

4)-②

4)-③

4)-④

4)-⑤

4)-⑥

4)-⑦

4)-⑧

此外还有两条直线有交叉的情况；

设p_o表示左右均无目标像素的单像素；

p_s表示连续像素个数≥2的像素组；

整型常量R表示目标之间可以重叠的像素数；

整型常量L表示一条直线内不包含重叠像素的最少像素数；

整型变量W(P₀)表示单像素p_o在行或列内的位置，行内顺序为从左至右，列内顺序为从上到下；

整型变量S(P_s)表示p_s在行或列内的起始位置，行内顺序为从左至右，列内顺序为从上到下；

整型变量E(P_s)表示p_s在行或列内的终止位置；

数组A_x用于存储x轴方向上的所有p_o和p_s信息，结构包含：目标类型、单像素位置，其中目标类型为p_o或p_s；

数组A_y用于存储y轴方向上的所有p_o和p_s信息，结构包含：目标类型、单像素位置、像素组在列内的起始点位置，其中目标类型为p_o或p_s；

整型变量R(P_o)、R(P_s)分别表示单像素p_o、像素组p_s在数组A_x内的行号；

整型变量C(P_o)、C(P_s)分别表示p_o、p_s在A_y内的列号；

数组AL表示存储一条直线内的所有p_o和p_s信息；

数组GL用于以(θ,ρ)的形式存储检索到的所有直线信息；

第四灰度图像中的直线具备以下特点：

4)-(一)：直线都是由p_o和p_s组成；

4)-(二)：θ角越接近或其线条中的p_o和p_s越多；θ角等于或时，直线仅由p_o组成；θ角越接近0、或π，直线中p_o和p_s越少，但p_s内的像素个数越多；θ角等于0、或π时，直线仅由p_s组成；

4)-(三)：交叉的直线，共用其中一个p_o或p_s；

4)-(四)：直线都是由p_o和p_s同时沿x轴、y轴两个方向连接形成；x轴方向上直线的长度逐渐增加，y轴方向上逐渐产生一定的斜率，或者，y轴方向上直线的长度逐渐增加，x轴方向上逐渐产生一定的斜率；

4)-(五)：当θ角为0、或π时，直线仅包含一个p_s；

令R＝1，以检索的直线，第四灰度图像中的直线提取通过下述步骤实现：

4)-⑴从第四灰度图像的左上角像素开始，按照逐行、逐列的顺序，即从左到右、从上到下的顺序，依次扫描第四灰度图像中的每一个像素点，来检索第四灰度图像中的目标点，目标点即边缘点；提取到的结果保存在A_x和A_y中；

4)-⑵取出A_x中的一个目标A_t，将A_t插入数组AL，其中，第一次抽取的目标A_t是随机抽取，以后抽取依照下步即第4)-⑶步所示，将抽取的目标A_t按照抽取的先后顺序插入数组AL；

4)-⑶沿x轴负方向，查找A_t的下一个目标A_t+1，A_t+1必须同时符合以下条件：

4)-⑶-①R(A_t+1)＝R(A_t)+1；

4)-⑶-②A_t为p_o时：

4)-⑶-②-Ⅰ：如果A_t+1为单像素，W(A_t+1)＝W(A_t)-R；

4)-⑶-②-Ⅱ：如果A_t+1为像素组，W(A_t)-R≤E(A_t+1)≤W(A_t)；

4)-⑶-③A_t为像素组时：

4)-⑶-③-Ⅰ：如果A_t+1为单像素，W(A_t+1)＝S(A_t)-R；

4)-⑶-③-Ⅱ：如果A_t+1为像素组，S(A_t)-R≤E(A_t+1)≤S(A_t)；

4)-⑶-④A_t+1未访问过；

如果找到A_t+1，将A_t+1标记为已访问且将其插入数组AL；令A_t＝A_t+1，递归执行该步骤4)-⑶处理过程；如果找不到A_t+1，结束查找；

4)-⑷令数组AL中的第一个元素为AL₀，最后一个元素为AL_n，如果数组AL中只有一个元素时：若该元素长度大于等于L，该元素长度即像素数，此条直线(此条直线为AL中该元素(θ,ρ)所表示的一条直线)起点为S(AL₀)，终点为E(AL₀)；若长度小于L,此次检索失败；如果数组AL中不少于2个元素，此条直线起点为E(AL₀)，终点为S(AL_n)；

4)-⑸沿x轴正方向，查找A_t的下一个目标A_t+1，A_t+1必须同时符合以下条件：

4)-⑸-①R(A_t+1)＝R(A_t)+1；

4)-⑸-②A_t为单像素时：

4)-⑸-②-Ⅰ：如果A_t+1为单像素，W(A_t+1)＝W(A_t)+R；

4)-⑸-②-Ⅱ：如果A_t+1为像素组，W(A_t)≤S(A_t+1)≤W(A_t)+R；

4)-⑸-③A_t为像素组时：

4)-⑸-③-Ⅰ：如果A_t+1为单像素，W(A_t+1)＝E(A_t)+R；

4)-⑸-③-Ⅱ：如果A_t+1为像素组，E(A_t)≤S(A_t+1)≤E(A_t)+R；

4)-⑸-④A_t+1未访问过；

如果找到A_t+1，将A_t+1标记为已访问且将其插入数组AL；令A_t＝A_t+1，递归执行该步骤4)-⑸的处理过程；如果找不到A_t+1，结束查找，执行下一步4)-⑹；

4)-⑹令数组AL中的第一个元素为AL₀，最后一个元素为AL_n；最终提取的直线分别以下述三种情况存储：

4)-⑹-①如果数组AL中只有一个元素且其长度大于等于L，该元素长度即为像素数，此条直线起点为S(AL₀)，终点为E(AL₀)；

4)-⑹-②如果数组AL中的元素数量≥2且E(AL₀)-S(AL₀)≥L，则此条直线起点为S(AL₀)，终点为E(AL_n)；

4)-⑹-③如果数组AL中的元素数量≥2且E(AL₀)-S(AL₀)＜L，放弃此条直线；

4)-⑺将检索到的直线插入到数组GL中,清空数组AL；

4)-⑻抽取数组A_x中一个未访问过的目标，令其等于A_t,递归执行上述对第四灰度图像中的直线提取步骤中的第4)-⑶步，直至数组A_x中所有目标都被标记为已访问；

4)-⑼扫描数组A_y中的目标，检索和内的直线；取出数组A_y中的一个目标A_t，将A_t插入数组AL；

4)-⑽沿y轴正方向，查找A_t的下一个目标A_t+1，A_t+1必须同时符合以下条件：

4)-⑽-①C(A_t+1)＝C(A_t)+1；

4)-⑽-②A_t为p_o时：

4)-⑽-②-Ⅰ：如果A_t+1为单像素，W(A_t+1)＝W(A_t)-R；

4)-⑽-②-Ⅱ：如果A_t+1为像素组，W(A_t)+R≤E(A_t+1)≤W(A_t)；

4)-⑽-③A_t为像素组时：

4)-⑽-③-Ⅰ：如果A_t+1为单像素，W(A_t+1)＝S(A_t)-R；

4)-⑽-③-Ⅱ：如果A_t+1为像素组，S(A_t)+R≤E(A_t+1)≤S(A_t)；

4)-⑽-④A_t+1未访问过；

如果找到A_t+1，将A_t+1标记为已访问且将其插入数组AL；令A_t＝A_t+1，递归执行该步4)-⑽；如果找不到A_t+1，结束查找，执行下一步4)-⑾；

4)-⑾令数组AL中的第一个元素为AL₀，最后一个元素为AL_n，如果数组AL中只有一个元素且该元素长度小于等于L，该元素长度即为像素数，则此条直线起点为S(AL₀)，终点为E(AL₀)；否则，如果数组AL中有至少2个元素，此条直线起点为E(AL₀)，终点为S(AL_n)；

4)-⑿沿y轴负方向，查找A_t的下一个目标A_t+1，A_t+1必须同时符合以下条件：

4)-⑿-①C(A_t+1)＝C(A_t)+1；

4)-⑿-②A_t为p_o时：

4)-⑿-②-Ⅰ：如果A_t+1为单像素，W(A_t+1)＝W(A_t)+R；

4)-⑿-②-Ⅱ：如果A_t+1为像素组，W(A_t)-R≤S(A_t+1)≤W(A_t)；

4)-⑿-③A_t为像素组时：

4)-⑿-③-Ⅰ：如果A_t+1为单像素，W(A_t+1)＝S(A_t)+R；

4)-⑿-③-Ⅱ：如果A_t+1为像素组，E(A_t)-R≤S(A_t+1)≤E(A_t)；

4)-⑿-④A_t+1未访问过；

如果找到A_t+1，将A_t+1标记为已访问且将其插入AL；令A_t＝A_t+1，递归执行该步4)-⑿；如果找不到A_t+1，结束查找，执行下一步4)-⒀；

4)-⒀令数组AL中的第一个元素为AL₀，最后一个元素为AL_n；最终提取的直线以下述三种情况存储：

4)-⒀-①如果数组AL中只有一个元素且其长度大于等于L，该元素长度即为像素数，此条直线起点为S(AL₀)，终点为E(AL₀)；

4)-⒀-②如果数组AL中的元素数量≥2且E(AL₀)-S(AL₀)≥L，则此条直线起点为S(AL₀)，终点为E(AL_n)；

4)-⒀-③如果AL中的元素数量≥2且E(AL₀)-S(AL₀)＜L，放弃此条直线；

4)-⒁将AL中的内容以Hough变换方法处理求出θ角和ρ值，插入到数组GL中,清空数组AL；

4)-⒂抽取A_y中一个未访问过的目标，令其等于A_t,递归执行上述第4)-⑽步，直至A_y中所有目标都被标记为已访问；

至此，直线的提取全部完成，数组GL中存储了第四灰度图像中0≤θ≤π的所有直线；

步骤5)、直线的归并

设阈值L_R、θ_R、ρ_R分别表示两条直线的端点、θ角、ρ值的容差，令直线L₁的起点、终点、θ角和ρ值分别为S₁、E₁、θ₁和ρ₁，直线L₂的起点、终点、θ角、ρ值分别为S₂、E₂、θ₂和ρ₂；从第四灰度图像中提取到的直线会有部分具备以下特点之一：

5)-⑴θ₁＝θ₂，ρ₁＝ρ₂，MIN{|S₁-E₂|,|S₂-E₁|}≤L_R；

5)-⑵|θ₁-θ₂|≤θ_R，|ρ₁-ρ₂|≤ρ_R；

5)-⑶ρ₁＝ρ₂，|θ₁-θ₂|≤θ_R；

5)-⑷θ₁＝θ₂，|ρ₁-ρ₂|≤ρ_R；

这四类直线需要进行归并，以减少重复的直线；

设合并后的直线为L₃，其起点、终点、θ角和ρ值为S₃、E₃、θ₃和ρ₃，针对以上四种情况，直线归并如下处理：

5)-①当MIN{|S₁-E₂|,|S₂-E₁|}＝|S₁-E₂|时，取E₃＝E₁，S₃＝S₂；否则，当MIN{|S₁-E₂|,|S₂-E₁|}≠|S₁-E₂|时，E₃＝E₂，S₃＝S₁；

5)-②当ρ₁＝ρ₂，|θ₁-θ₂|≤θ_R时，取ρ₃＝ρ₁＝ρ₂,

5)-③当L₁‖L₂时，取θ₃＝θ₁＝θ₂,符号“‖”表示直线间的平行；

四类直线归并前，两条直线分别表示为L₁(S₁,E₁)和L₂(S₂,E₂)，归并后的直线表示为L₃(S₃,E₃)；

步骤6)、图像中的角点提取

(61)：角度阈值过滤

定义角度阈值T_θ，判断两条直线的交叉角度在之间时，认为是直角，直线的交点即为建筑物的角点；

(62)：间断直线过滤

在背景复杂且噪声较多的遥感影像中，部分边缘直线不会被完整的获取到，出现间断的线段，设置阈值T_L，确定在间断长度为阈值T_L以内时，认为此角点可用，否则，间断的直线不做为直角边使用；

(63)：角分裂；

根据直线交叉的情况，分裂成多个直角；

(64)：计算角点；

已知两条直线L₁:A₁x+B₁y+C₁＝0；L₂:A₂x+B₂y+C₂＝0，两直线交点是坐标为(x,y)的像素点；

通过式(19)求出两直线的交点即角点的坐标x值、y值；

为了便于后续的角点归并和排序，在记录角点时，还需要考虑到角的开口方向；在得到角点的坐标值后，令：角的两条边为e₁和e₂，以角点为起点、远离角点的端点为终点的两个向量分别为和角方向如下表示：

的方向即为角的方向，所述角以θ_r表示，角的范围为：0≤θ_r≤2π；

(65)：角点过滤和归并；

由于遥感影像中的复杂背景，提取目标周围会出现一些伪边缘，需要进行角点过滤和归并：

如果两角点的开口方向相同：

令角点P1的两边长和夹角θ_r分别为E₁、E₂、θ_r1，角点P2的两边长和夹角分别为E₃、E₄、θ_r2，设置距离差阈值T_D、角度容差T_θ；当角点P1和角点P2的距离D₁≤T_D且|θ_r1-θ_r2|≤T_θ时，对两角点进行归并；如果E₁E₂≥E₃E₄，则删除角点P2；如果E₁E₂＜E₃E₄，删除角点P1；

如果两角点的开口方向相对，且两角点之间出现了第三个角点P3：

令角点P3的两边长、夹角分别为E₅、E₆、θ_r3，当角点P1和角点P3的距离D₂＞T_D，且|θ_r1-θ_r3|≤T_θ时，以角点P2的θ_r2方向为检索方向，在θ_r2方向上寻找θ_r2-T_θ≤θ_r≤θ_r2+T_θ的角点，找到角点P1和角点P3后，如果E₁E₂≥E₅E₆，则删除角点P3；否则，当E₁E₂小于E₅E₆时，删除角点P1；

步骤7)、归属同一建筑物的角点匹配

建筑物角点的检索可以使用如下步骤完成：

(71)从角点集A中取一个角点A_s,如果A_s已访问过，则重新取；角点集A为经过角点过滤和归并后得到的角点集合；

(72)：沿角点A_s的一条边的方向寻找下一个与角点A_s匹配的角点A_s+1，匹配的条件是：

角点A_s+1的一条边与的方向相反，且误差在设定的角度容差θ_E以内；

(73)：将当前角点A_s标记为已访问，并将其保存在数组B内；

(74)：令A_s＝A_s+1，e₁＝e₄，递归执行步骤(72),直到匹配到最初的A_s角点为止；

(75)：遍历角点集A中所有角点，直到所有角点均为已访问；匹配结束；

对于未完成循环匹配的一组角点，由如下步骤处理推断出其余角点，构成由至少4个角点组成的矩形建筑：

⑴单独角点：由角点A_s在角点集A中匹配不到A_s+1时，以角点A_s的两条边以及角方向的三个端点，构建一个四边形的建筑物轮廓；

⑵两个角点：角点A_s在角点集A中匹配到A_s+1，而A_s+1匹配不到其它角点时，以A_s和A_s+1的两条边的端点为另两个角点，构建一个四边形的建筑物轮廓；

⑶多于两个角点：匹配到两个以上角点，但没有形成环路时，以最初的角点A_s的边与最后一个角点的边的交点做为最后一个角点；

步骤8)、扫描结束，输出处理后的图像，结束程序。