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

摘要

本发明提供了一种高分辨率遥感影像中的建筑物快速提取方法，先对图像进行灰度提取，然后通过高斯滤波进行去噪，接着使用基于Otsu方法的Canny自适应处理方法提取到最符合实际情况的边缘二值图像，再使用基于Freeman链码的方法快速提取二值图像中的直线，最后取到角点，连接形成建筑物轮廓，最终完成高分辨率遥感影像中建筑物的快速自动提取；本发明提供的高分辨率遥感影像中的建筑物快速提取方法，可实现快速提取遥感影像中的建筑物，处理过程中不需人为干预，提取速度快，精度高。

主权项

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个元素的模板值分别为m1、m2、m3、m4、m5、m6、m7、m8、m9；第一灰度图像中，除左、上、右、下四条宽度为一像素的边缘外任一像素点P点的8邻像素，即紧邻像素点P的左上、上、右上、左、右、左下、下、右下8个方向上的像素；像素点P与8邻像素的灰度值按照第一灰度图像中从左至右、从上到下的顺序依次为g1、g2、g3、g4、g5、g6、g7、g8、g9；用该高斯模板U(x，y)按照从左到右、从上到下的顺序依次扫描第一灰度图像中除左、上、右、下四条宽度为一像素的边缘外的每一个像素P，各像素点P由高斯模板U(x，y)做高斯滤波的方法见式(3)；所述的U0表示以高斯模板U(m，n)做高斯滤波后所到的像素点P的灰度值；S为中间变量，表示像素点P和其8邻像素的灰度值的总和；按照式(3)依从左到右、从上到下的顺序对第一灰度图像中像素点P进行依次扫描时，对当前所扫描的像素点P由高斯模板U(x，y)所确定的邻域内所有像素的加权平均灰度值替代当前像素点P的灰度值，扫描结束后得到第二灰度图像；步骤3)、使用基于Otsu方法的自适应Canny算法提取图像中的边缘；31)：用一阶偏导的有限差分来计算梯度的幅值和方向；采用Canny算法,其中所采用的卷积模板为式(4)；其中M1、M2分别为第二灰度图像x、y轴方向上的一阶偏导数矩阵,依照所述一阶偏导数矩阵M1、M2按照从左到右、从上到下的顺序依次扫描第二灰度图像中的每一个像素点，得到每个像素点x、y方向上的梯度幅值Gx、Gy：Gx＝[f(x+1,y)‑f(x,y)+f(x+1,y+1)‑f(x,y+1)]/2    式(5)Gy＝[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)；其中，G5为当前处理像素的梯度幅值，Gk(1≤k≤4,6≤k≤9)表示8邻像素的梯度幅值，Gp表示梯度方向上前一个像素的梯度幅值，Gn表示梯度方向上后一个像素的梯度幅值；当梯度方向α的值在上述相应的区间内时，当且仅当目前所处理像素的梯度幅值G5大于前述①～④中各插值方法计算出的Gp和Gn时，G5为极大值，也就表明当前处理像素是边缘点；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₂为最优高阈值；非边缘点总个数N1(T1,T2)为：其中，Ni是梯度幅值等于i的像素的个数；非边缘点所占比例值ω1(T1,T2)为：潜在边缘点总个数N2(T1,T2)为：潜在边缘点所占比例值ω2(T1,T2)为：真边缘点总个数N3(T1,T2)为：真边缘点所占比例值ω3(T1,T2)为：非边缘点的平均梯度μ1(T1,T2)、潜在边缘点的平均梯度μ2(T1,T2)和真边缘点的平均梯度μ3(T1,T2)如下：第三灰度图像总的梯度均值μs为：μs＝μ1(T1,T2)ω1(T1,T2)+μ2(T1,T2)ω2(T1,T2)+μ3(T1,T2)ω3(T1,T2)    式(17)非边缘点、潜在边缘点和真边缘点这三类像素的类间方差σ2(T1,T2)为：σ2(T1,T2)＝ω1(T1,T2)(μ1(T1,T2)‑μS)2+ω2(T1,T2)(μ2(T1,T2)‑μS)2+ω3(T1,T2)(μ3(T1,T2)‑μS)2                                              式(18)使得σ2(T1,T2)的输入参数T1、T2在总范围[1,G‑1]依次取值，当类间方差σ2(T1,T2)为最大值时，此时的T1、T2的值即为最优低阈值T1和最优高阈值T2；34)：用双阈值算法检测和边缘连接；根据上一步33)中得到的最优低阈值T1和最优高阈值T2对第三灰度图像进行过滤，像素点的梯度幅值如果大于高阈值则认为当前像素点必然是边缘点，即真边缘点，如果当前处理像素点的梯度幅值小于低阈值，则认为当前像素点必然不是边缘点；处于低阈值和高阈值之间的像素点是潜在的边缘点，即假边缘；将第三灰度图像根据高阈值过滤得到的边缘点链接成轮廓，当到达轮廓的端点时，在端点的8邻点寻找潜在的边缘点，再根据所找到的潜在的边缘点,在其8邻点中循环处理收集新的边缘点，如此处理整个第三灰度图像，直至边缘闭合，得到第四灰度图像；步骤4)、基于Freeman链码的直线提取方法：单条直线L与Y轴的夹角θ范围为0～π,共有八种情况的直线L：4)‑①θ＝0或θ＝π；4)‑②4)‑③4)‑④4)‑⑤4)‑⑥4)‑⑦4)‑⑧此外还有两条直线有交叉的情况；设po表示左右均无目标像素的单像素；ps表示连续像素个数≥2的像素组；整型常量R表示目标之间可以重叠的像素数；整型常量L表示一条直线内不包含重叠像素的最少像素数；整型变量W(P0)表示单像素po在行或列内的位置，行内顺序为从左至右，列内顺序为从上到下；整型变量S(Ps)表示ps在行或列内的起始位置，行内顺序为从左至右，列内顺序为从上到下；整型变量E(Ps)表示ps在行或列内的终止位置；数组Ax用于存储x轴方向上的所有po和ps信息，结构包含：目标类型、单像素位置，其中目标类型为po或ps；数组Ay用于存储y轴方向上的所有po和ps信息，结构包含：目标类型、单像素位置、像素组在列内的起始点位置，其中目标类型为po或ps；整型变量R(Po)、R(Ps)分别表示单像素po、像素组ps在数组Ax内的行号；整型变量C(Po)、C(Ps)分别表示po、ps在Ay内的列号；数组AL表示存储一条直线内的所有po和ps信息；数组GL用于以(θ,ρ)的形式存储检索到的所有直线信息；第四灰度图像中的直线具备以下特点：4)‑(一)：直线都是由po和ps组成；4)‑(二)：θ角越接近或其线条中的p_o和p_s越多；θ角等于或时，直线仅由p_o组成；θ角越接近0、或π，直线中p_o和p_s越少，但p_s内的像素个数越多；θ角等于0、或π时，直线仅由p_s组成；4)‑(三)：交叉的直线，共用其中一个po或ps；4)‑(四)：直线都是由po和ps同时沿x轴、y轴两个方向连接形成；x轴方向上直线的长度逐渐增加，y轴方向上逐渐产生一定的斜率，或者，y轴方向上直线的长度逐渐增加，x轴方向上逐渐产生一定的斜率；4)‑(五)：当θ角为0、或π时，直线仅包含一个p_s；令R＝1，以检索的直线，第四灰度图像中的直线提取通过下述步骤实现：4)‑⑴从第四灰度图像的左上角像素开始，按照逐行、逐列的顺序，即从左到右、从上到下的顺序，依次扫描第四灰度图像中的每一个像素点，来检索第四灰度图像中的目标点，目标点即边缘点；提取到的结果保存在Ax和Ay中；4)‑⑵取出Ax中的一个目标At，将At插入数组AL，其中，第一次抽取的目标At是随机抽取，以后抽取依照下步即第4)‑⑶步所示，将抽取的目标At按照抽取的先后顺序插入数组AL；4)‑⑶沿x轴负方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：4)‑⑶‑①R(At+1)＝R(At)+1；4)‑⑶‑②At为po时：4)‑⑶‑②‑Ⅰ：如果At+1为单像素，W(At+1)＝W(At)‑R；4)‑⑶‑②‑Ⅱ：如果At+1为像素组，W(At)‑R≤E(At+1)≤W(At)；4)‑⑶‑③At为像素组时：4)‑⑶‑③‑Ⅰ：如果At+1为单像素，W(At+1)＝S(At)‑R；4)‑⑶‑③‑Ⅱ：如果At+1为像素组，S(At)‑R≤E(At+1)≤S(At)；4)‑⑶‑④At+1未访问过；如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步骤4)‑⑶处理过程；如果找不到At+1，结束查找；4)‑⑷令数组AL中的第一个元素为AL0，最后一个元素为ALn，如果数组AL中只有一个元素时：若该元素长度大于等于L，该元素长度即像素数，此条直线(此条直线为AL中该元素(θ,ρ)所表示的一条直线)起点为S(AL0)，终点为E(AL0)；若长度小于L,此次检索失败；如果数组AL中不少于2个元素，此条直线起点为E(AL0)，终点为S(ALn)；4)‑⑸沿x轴正方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：4)‑⑸‑①R(At+1)＝R(At)+1；4)‑⑸‑②At为单像素时：4)‑⑸‑②‑Ⅰ：如果At+1为单像素，W(At+1)＝W(At)+R；4)‑⑸‑②‑Ⅱ：如果At+1为像素组，W(At)≤S(At+1)≤W(At)+R；4)‑⑸‑③At为像素组时：4)‑⑸‑③‑Ⅰ：如果At+1为单像素，W(At+1)＝E(At)+R；4)‑⑸‑③‑Ⅱ：如果At+1为像素组，E(At)≤S(At+1)≤E(At)+R；4)‑⑸‑④At+1未访问过；如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步骤4)‑⑸的处理过程；如果找不到At+1，结束查找，执行下一步4)‑⑹；4)‑⑹令数组AL中的第一个元素为AL0，最后一个元素为ALn；最终提取的直线分别以下述三种情况存储：4)‑⑹‑①如果数组AL中只有一个元素且其长度大于等于L，该元素长度即为像素数，此条直线起点为S(AL0)，终点为E(AL0)；4)‑⑹‑②如果数组AL中的元素数量≥2且E(AL0)‑S(AL0)≥L，则此条直线起点为S(AL0)，终点为E(ALn)；4)‑⑹‑③如果数组AL中的元素数量≥2且E(AL0)‑S(AL0)＜L，放弃此条直线；4)‑⑺将检索到的直线插入到数组GL中,清空数组AL；4)‑⑻抽取数组Ax中一个未访问过的目标，令其等于At,递归执行上述对第四灰度图像中的直线提取步骤中的第4)‑⑶步，直至数组Ax中所有目标都被标记为已访问；4)‑⑼扫描数组A_y中的目标，检索和内的直线；取出数组A_y中的一个目标A_t，将A_t插入数组AL；4)‑⑽沿y轴正方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：4)‑⑽‑①C(At+1)＝C(At)+1；4)‑⑽‑②At为po时：4)‑⑽‑②‑Ⅰ：如果At+1为单像素，W(At+1)＝W(At)‑R；4)‑⑽‑②‑Ⅱ：如果At+1为像素组，W(At)+R≤E(At+1)≤W(At)；4)‑⑽‑③At为像素组时：4)‑⑽‑③‑Ⅰ：如果At+1为单像素，W(At+1)＝S(At)‑R；4)‑⑽‑③‑Ⅱ：如果At+1为像素组，S(At)+R≤E(At+1)≤S(At)；4)‑⑽‑④At+1未访问过；如果找到At+1，将At+1标记为已访问且将其插入数组AL；令At＝At+1，递归执行该步4)‑⑽；如果找不到At+1，结束查找，执行下一步4)‑⑾；4)‑⑾令数组AL中的第一个元素为AL0，最后一个元素为ALn，如果数组AL中只有一个元素且该元素长度小于等于L，该元素长度即为像素数，则此条直线起点为S(AL0)，终点为E(AL0)；否则，如果数组AL中有至少2个元素，此条直线起点为E(AL0)，终点为S(ALn)；4)‑⑿沿y轴负方向，查找At的下一个目标At+1，At+1必须同时符合以下条件：4)‑⑿‑①C(At+1)＝C(At)+1；4)‑⑿‑②At为po时：4)‑⑿‑②‑Ⅰ：如果At+1为单像素，W(At+1)＝W(At)+R；4)‑⑿‑②‑Ⅱ：如果At+1为像素组，W(At)‑R≤S(At+1)≤W(At)；4)‑⑿‑③At为像素组时：4)‑⑿‑③‑Ⅰ：如果At+1为单像素，W(At+1)＝S(At)+R；4)‑⑿‑③‑Ⅱ：如果At+1为像素组，E(At)‑R≤S(At+1)≤E(At)；4)‑⑿‑④At+1未访问过；如果找到At+1，将At+1标记为已访问且将其插入AL；令At＝At+1，递归执行该步4)‑⑿；如果找不到At+1，结束查找，执行下一步4)‑⒀；4)‑⒀令数组AL中的第一个元素为AL0，最后一个元素为ALn；最终提取的直线以下述三种情况存储：4)‑⒀‑①如果数组AL中只有一个元素且其长度大于等于L，该元素长度即为像素数，此条直线起点为S(AL0)，终点为E(AL0)；4)‑⒀‑②如果数组AL中的元素数量≥2且E(AL0)‑S(AL0)≥L，则此条直线起点为S(AL0)，终点为E(ALn)；4)‑⒀‑③如果AL中的元素数量≥2且E(AL0)‑S(AL0)＜L，放弃此条直线；4)‑⒁将AL中的内容以Hough变换方法处理求出θ角和ρ值，插入到数组GL中,清空数组AL；4)‑⒂抽取Ay中一个未访问过的目标，令其等于At,递归执行上述第4)‑⑽步，直至Ay中所有目标都被标记为已访问；至此，直线的提取全部完成，数组GL中存储了第四灰度图像中0≤θ≤π的所有直线；步骤5)、直线的归并设阈值LR、θR、ρR分别表示两条直线的端点、θ角、ρ值的容差，令直线L1的起点、终点、θ角和ρ值分别为S1、E1、θ1和ρ1，直线L2的起点、终点、θ角、ρ值分别为S2、E2、θ2和ρ2；从第四灰度图像中提取到的直线会有部分具备以下特点之一：5)‑⑴θ1＝θ2，ρ1＝ρ2，MIN{|S1‑E2|,|S2‑E1|}≤LR；5)‑⑵|θ1‑θ2|≤θR，|ρ1‑ρ2|≤ρR；5)‑⑶ρ1＝ρ2，|θ1‑θ2|≤θR；5)‑⑷θ1＝θ2，|ρ1‑ρ2|≤ρR；这四类直线需要进行归并，以减少重复的直线；设合并后的直线为L3，其起点、终点、θ角和ρ值为S3、E3、θ3和ρ3，针对以上四种情况，直线归并如下处理：5)‑①当MIN{|S1‑E2|,|S2‑E1|}＝|S1‑E2|时，取E3＝E1，S3＝S2；否则，当MIN{|S1‑E2|,|S2‑E1|}≠|S1‑E2|时，E3＝E2，S3＝S1；5)‑②当ρ₁＝ρ₂，|θ₁‑θ₂|≤θ_R时，取ρ₃＝ρ₁＝ρ₂,5)‑③当L₁‖L₂时，取θ₃＝θ₁＝θ₂,符号“‖”表示直线间的平行；四类直线归并前，两条直线分别表示为L1(S1,E1)和L2(S2,E2)，归并后的直线表示为L3(S3,E3)；步骤6)、图像中的角点提取(61)：角度阈值过滤定义角度阈值T_θ，判断两条直线的交叉角度在之间时，认为是直角，直线的交点即为建筑物的角点；(62)：间断直线过滤在背景复杂且噪声较多的遥感影像中，部分边缘直线不会被完整的获取到，出现间断的线段，设置阈值TL，确定在间断长度为阈值TL以内时，认为此角点可用，否则，间断的直线不做为直角边使用；(63)：角分裂；根据直线交叉的情况，分裂成多个直角；(64)：计算角点；已知两条直线L1:A1x+B1y+C1＝0；L2:A2x+B2y+C2＝0，两直线交点是坐标为(x,y)的像素点；通过式(19)求出两直线的交点即角点的坐标x值、y值；为了便于后续的角点归并和排序，在记录角点时，还需要考虑到角的开口方向；在得到角点的坐标值后，令：角的两条边为e₁和e₂，以角点为起点、远离角点的端点为终点的两个向量分别为和角方向如下表示：的方向即为角的方向，所述角以θ_r表示，角的范围为：0≤θ_r≤2π；(65)：角点过滤和归并；由于遥感影像中的复杂背景，提取目标周围会出现一些伪边缘，需要进行角点过滤和归并：如果两角点的开口方向相同：令角点P1的两边长和夹角θr分别为E1、E2、θr1，角点P2的两边长和夹角分别为E3、E4、θr2，设置距离差阈值TD、角度容差Tθ；当角点P1和角点P2的距离D1≤TD且|θr1‑θr2|≤Tθ时，对两角点进行归并；如果E1E2≥E3E4，则删除角点P2；如果E1E2＜E3E4，删除角点P1；如果两角点的开口方向相对，且两角点之间出现了第三个角点P3：令角点P3的两边长、夹角分别为E5、E6、θr3，当角点P1和角点P3的距离D2＞TD，且|θr1‑θr3|≤Tθ时，以角点P2的θr2方向为检索方向，在θr2方向上寻找θr2‑Tθ≤θr≤θr2+Tθ的角点，找到角点P1和角点P3后，如果E1E2≥E5E6，则删除角点P3；否则，当E1E2小于E5E6时，删除角点P1；步骤7)、归属同一建筑物的角点匹配建筑物角点的检索可以使用如下步骤完成：(71)从角点集A中取一个角点As,如果As已访问过，则重新取；角点集A为经过角点过滤和归并后得到的角点集合；(72)：沿角点A_s的一条边的方向寻找下一个与角点A_s匹配的角点A_s+1，匹配的条件是：角点A_s+1的一条边与的方向相反，且误差在设定的角度容差θ_E以内；(73)：将当前角点As标记为已访问，并将其保存在数组B内；(74)：令As＝As+1，e1＝e4，递归执行步骤(72),直到匹配到最初的As角点为止；(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)、扫描结束，输出处理后的图像，结束程序。