基于立体约束的遥感影像道路提取方法

专利名称：基于立体约束的遥感影像道路提取方法
技术领域：
本发明属于遥感技术领域，具体涉及一种基于立体约束的遥感影像道路提取方 法。
背景技术：
道路是一种重要的专题信息，其对地学分析、城市管理有着重要的价值，是一类基 础的专题信息数据，对各种地理信息系统的建设与管理有重要的意义，特别是城市道路对 于城市的日常管理有着重要价值。利用遥感影像进行道路信息提取是一种快速、且现势性 好的方式。因此成为人们研究的热点。目前利用遥感影像提取道路的方法有人工方法、 半自动跟踪方法以及自动模式分类的方法。其中人工方法包括直接在单幅遥感影像上的 提取、以及在立体像对(多为航空影像对)下基于立体观察的提取方式；半自动跟踪主要 是在单幅遥感影像上利用给定的种子点进行跟踪提取道路；这两种方法都需要大量的人 工操作，且效率低。利用遥感影像上道路的光谱及形态特征进行道路的自动提取是一种重 要方法，但是现有的基于光谱信息的分类识别方法并不是针对道路信息提取而研究的；而 基于形态的道路自动识别提取对于一般分辨率(空间分辨率20m，或更低)的农村地区的 遥感影像比较适用，对于高分辨率(分辨率达到5m及以上)且城市地区的遥感影像，其对 道路的识别提取结果不能满足后续成图及应用要求，由于城市地区建筑物众多其形态与道 路有相似之处，对道路提取有很大的干扰，如何消除建筑物对道路提取的影响，在城市遥感 影像的道路提取中是个问题。目前，在高分辨率遥感影像分类识别信息提取方面，德国的 eCognition系统的分类方法是比较好的一种，其基于多尺度分割，分类时需要有分类协议， 每一次分类的分类协议用于其它分类时需经过修改，最终实现分类过程的自动化，这里的 分类协议类似于知识支持。其特点是分类的细化，但对分类协议有较大依赖(有协议库支 持较好)，并且对道路信息提取时效率不高，需要大量的交互，特别是在城市地区。为提高道 路提取效率，避免提取过程中的人工干扰，需要实现一种自动的道路信息提取方法，以解决 城市地区高分辨率遥感影像的道路提取。

发明内容
本发明的目的在于提供一种基于立体约束的遥感影像道路提取方法，针对上述的 建筑物干扰道路提取的问题，产生利用遥感影像立体像对参数，根据道路候选特征线在立 体像对上的同名信息，计算道路候选特征线的高度信息，利用该信息构成一个约束条件，对 于高度信息超出该地区地面高度值较多的特征线，认为是建筑物特征线予以剔除，有效的 保留真实的地面道路信息；在此约束的基础上，进行道路边线的自动选择提取，得到道路专 题信息。本发明的内容为一种基于立体约束的遥感影像道路提取方法，首先对立体像对中 的一幅影像利用边缘提取算子(CANNY算子)进行边缘提取、再利用细化算法(hilditch算 法)构成单像素宽度的边缘；在上面提取结果的基础上，使用相位编组方法对边缘长度和弯曲程度进行约束以去除噪声边缘及曲线干扰，提取待选道路边缘线；为进一步消除建筑 物边线等对道路提取的干扰，在待选道路边缘线上选取若干点先通过仿射变换确定这些点 在立体像对另一影像上的匹配初值，再利用灰度相关匹配得到这些点另一幅影像上的准确 同名点位置；利用立体像对已知位置及姿态关系参数(高分辨率遥感影像立体像对的位置 及姿态关系可以是严格的外方位元素表示、也可以是利用有理多项式系数表示的参数)，根 据匹配得到的每对同名点的坐标及两幅影像的位置、姿态参数，前方交会确定上述边线点 的高度信息；从道路和建筑物一般具有明显的高度差异出发，建立道路边线高度约束准则， 利用所求高度信息对上述道路候选边缘进行的约束，从而进一步消除建筑物边缘的干扰， 实现遥感影像上最终道路边缘。在此过程中，在尽量少的人工干预和先验知识条件下，进行 道路信息提取，提取获得道路信息矢量数据。具体包括基于立体像对左影像的边缘提取和 细化获得单像素边缘的二值图像；利用改进相位编组方法对上面的细化结果进行形态约束 和矢量化，得到图像候选道路边缘线；基于灰度相关进行边缘线节点匹配，确定右影像上边 缘线节点的同名像点；基于像对RPC模型参数(或内、外方位元素)前方交会解算节点高 程；基于节点高程约束的道路边缘线确定得到最终道路提取结果；具体包含如下步骤(1)基于立体像对左影像的边缘提取和细化获得单像素边缘的二值图像；(2)利用改进相位编组方法对上述单像素边缘二值图像进行形态约束和矢量化， 得到图像候选道路边缘线；(3)基于灰度相关进行边缘线节点匹配，确定边缘线节点在右影像上的同名像占.
^ \\\ (4)基于像对RPC模型参数或内\外方位元素前方交会解算节点高程；(5)基于节点高程约束的道路边缘线确定得到道路提取结果。本发明中，步骤(1)所述的基于立体像对左影像的边缘提取和细化获得单像 素边缘的二值图像，基于CANNY算子的左影像边缘检测，以及在此边缘检测结果上基于 Hilditch细化方法的左影像单像素边缘二值图像生成，具体按如下步骤进行
V(G*/)①将图像f 与尺度为0的高斯函数G做卷积G*f，利用式〃 =|v(Gi[i/)|，对图像的 每个像素，计算局部边缘的法向量n ；②使用非最大值抑制方法找到法向量n方向的方向梯度最大值位置(x，y)，
|tg*/ = °实际做法是，沿着点(x，y)处法向量方向n，在(x，y)两侧各列出一定范围点 on，
的方向梯度值，利用这些散点，通过曲线拟合的方法可以确定梯度极大值，(x，y)的值可更 新为此位置；③利用式
lVG*/(x,y)HG*/(x + l,y)-G*/(x-l,y)l + \G*/(x,y + l)-G*/(x,y-l)jif^mm^
(x,y)的梯度幅值
9④用高低阈值hlOT，hhigh判断该点是否属于边缘，若，则确定为边缘 点；若且在8邻域与某些高阈值边缘点连接也认为是边缘点，将边缘点像素 值置为1，非边缘点置为0;
⑤将CANNY算子提取的边缘图像转化为单像素边缘二值图像，利用Hilditch细化 方法，完成轮廓细化，在步骤④获得的二值图像中，像素点的值为1或0，值为1的点，称为前 景点或边缘，即图形像素，值为0的点，称为背景点，在细化算法中，一个前景点满足下面的 条件才能被删除a.该点为边界点；b.该点非端点；c.该点非断点；d.该点的删除不应该 造成过度侵蚀。点P的8邻域以及邻域模板如下图所示(*表示0或1)。 Hilditch细化方法需要对整幅图像作多次遍历，在每一次遍历过程中，通过对像 素点八邻域与上述模板进行匹配，标记同时满足上述条件a d的像素点，在此次遍历结束 后，再将这些标记的像素点删除即由前景点转化为背景点，然后开始下一次遍历，直到在某 一次遍历过程中，没有像素点被标记，整个细化过程结束；经过细化运算后得到单像素边缘 二值图像，边缘像素点的灰度为0，非边缘像素点的灰度为1。本发明中，可以定义步骤(1)中①中的尺度参数o = 3，步骤④中的高低阈值hlOT —75, hhigh — 225。本发明中，步骤(2)所述的利用改进相位编组方法对上述单像素边缘二值图像矢 量化，得到图像边缘线，具体按如下步骤进行设计三个阈值来限定道路边缘，依次为边缘线的长度阈值，首末距离阈值和弯曲 程度阈值，长度阈值指边缘线上像素个数，首末距离阈值指首末两端点的距离，在边缘弯曲 甚至首尾接近时，首末距离阈值可以用来排除非道路边缘，弯曲程度阈值用首末点和线上 一结点的连线夹角表示，若夹角小于此阈值，则说明弯曲过大，应从道路边缘中排除；具体 为①将步骤(1)中提取的单像素边缘图像上，所有边缘像素点的状态设为0，同时确 定道路边缘线提取的长度阈值Th-L，首末距离阈值Th-D和弯曲程度阈值Th-curve ；②对单像素边缘图像上，状态为0的边缘像点逐像素判断，计算边缘像素8邻域中 灰度为0的点数m，如m= 1，则确定其为边缘线的一个端点；③从该端点起沿着该线段逐点记录坐标，记录到点列表pointList中，直到到达 另一个端点，长度L为该单像素边缘中包含的像素点个数，首末距离D为首末端点利用欧 式距离公式计算得到的结果，D =，如果D彡Th-D, L彡Th-L,则保
留该线段；并对该列线段进行平分，找到3分节点，利用式C0S《=<口,办〉计算每个分
a
b
节点(xiJyi)和首末点的连线夹角9i，设向量a= (Xi-X首，y「y首)，b = (Xi-X末，yi-y末)， 对每个分节点和首末点的连线夹角9 i进行判断，若0i<Th-CUrve阈值，则认为该线段 弯曲予以去除，如果不是则该线段予以保留，作为候选道路边缘；对跟踪判断过的边缘像素点，无论是否作为候选道路边缘，都将其标记为“已跟踪判断”即边缘像素点状态改为1 ；④对步骤(1)中所有未跟踪的边缘即状态为0，重复步骤②、③，直到所有边缘 线都完成上述判断，最后得到候选的矢量化的道路边缘线，这些候选边缘线放入线列表 lineList 中。本发明中，可以定义步骤(2)中①中的边缘线长度阈值Th-L为18/P到48/P像 素范围，P为该卫星图像的像素空间分辨率(一像素对应的实际地面边长)；首末距离阈值 Th-D为18/P到36/P像素范围；弯曲阈值为大于等于120度小于180度。本发明中，步骤(3)所述的基于灰度相关进行边缘线节点匹配，确定边缘线节点 在右影像上的同名像点，具体按如下步骤进行对步骤(2)中提取的候选道路边缘线进行三等分，加上端点共计四个分节点，利 用基于相关系数的灰度相关在右片上搜索出相应的同名点，为增强匹配的可靠性和速度， 通过仿射变换模型确定匹配的初始位置，初始位置的确定方法如下①采用一阶仿射变换模型来确定待匹配点在右片影像的同名点初始位置，
如； 选取若干左右影像同名像点3个以上，利用最小二乘方法解算以上6个参数(ai，
bi)；②依据这些系数求解匹配点的初始位置，利用式
以及
左影像上待匹配点的坐标(&，yi)，计算出其在右影像上同名点的粗略位置(乂，；③确定右片影像上粗略位置后，在以此位置为中心的19*19的搜索窗 口中寻找左影像点的匹配点的精确位置，左片匹配窗口大小为15*15，依据式
计算相关系数，并选取搜索窗口内相关系数最大的点作 V (=1/=1
为匹配最佳位置，即可得到精确的同名点位置(^，y》，式中x，y分别为匹配窗口和右片待 相关窗口内的像素灰度值，[和^分别是两窗口内的灰度均值，N为窗口中的像素个数；利用②、③可以确定每条候选道路边缘线上4个分节点在右影像上的同名点坐 标，并记录下来。本发明中，步骤(4)所述的基于像对RPC模型参数或内\外方位元素前方交会解 算节点高程，具体按如下步骤进行有理多项式函数 其中 0 ^ Ii1 ^ order，0 ^ n2 ^ order，0 ^ n3 ^ order 且 Wn3 ^ order，order 为多项式的阶数，一般认为对于线阵推扫式影像，取order = 3，已可满足精度的要求，这样 每个P (X，Y，Z)是一个20项的三阶多项式； 把& = — Cn = ^^代入上述有理多项式函数得到下式 立体像对中每幅影像的有理多项式函数RPC中的Cijk已知或内、外方位元素已知， 平移参数A、C0, X0,\、Z0和比例参数rs、cs, Xs、Ys、Zs也已知时，根据从步骤(3)中得到的 匹配同名点对，利用前方交会原理，每对同名点可以列出4个等式，解算出该点对应地面物 点的三维坐标(X，Y，Z)，其中Z坐标一般对应该物点的高程信息；一条边缘线三等分点可以得到四个分节点的高程坐标，将其高程坐标平均，得到 本发明中，步骤(5)所述的基于节点高程约束的道路边缘线确定得到道路提取结 果，具体按如下步骤进行根据所处理影像区域地面平均高程情况，给定道路高程约束阈值TH-h，如平均高 程大于道路高程阈值Havwage > TH-h，认为该四点所在边缘是非道路边缘如建筑边缘，保存 满足高程约束条件的道路边缘线在IineList列表中，对每条候选道路边缘线逐个按照步 骤(4)中的方法计算候选道路边缘线高程坐标平均值Harera一并判断其是否是道路边缘线， 对满足要求的道路边缘线矢量数据保留，存在列表RoadList中，并记录输出。本发明中，可以定义步骤(5)中的高程约束阈值TH-h为该地区平均高程值加上 7-8m。本发明的有益效果在于本发明从立体影像入手，对边缘提取、细化、经过形态约 束和矢量化后的道路边缘线的高程信息进行提取，建立对候选道路边缘线的分节点进行匹 配、根据同名点信息进行前方交会获得高程信息的模型，根据道路的高程特点，建立高程约 束条件，直接进行高程的约束，从候选道路边缘线中提取满足要求的道路边缘线，从而实现 遥感影像上的道路专题信息的快速自动提取，本发明受噪声和误差影响较少，可以较好去 除城市区域建筑物影像对道路提取的影响，鲁棒性强，在开始给出两立体影像5-6对同名 点位置后，其处理提取道路过程全自动进行，无须人工干预，能处理的数据广泛，耗时少。


图1是本发明的流程示意图2是原始quickbird立体像对的左右片影像图，其中(a)为左片影像，(b)为右 片影像；图3是左片边缘提取和细化操作结果图；其中(a)为边缘提取结果图，(b)为细 化操作结果图；图4是道路线提取结果，其中左片提取的矢量边缘线(暗灰色线)以及经过高程 约束最终保留的道路边缘线(亮色线)。
具体实施例方式下面的实施例是对本发明的进一步说明，而不是限制本发明的范围。实施例1 要对一高分辨率遥感影像立体像对通过立体约束进行道路边缘线的提取，获得道 路专题信息矢量数据。实验图像为Quickbird全色波段立体像对，其图像大小左影像为6961*5048、右 影像为6606*5419，图像像素分辨率为0. 6m ；由于图像太大，在实验中取了图像中的部分 A (左影像上)、B (右影像上)进行道路边缘线提取实验，A和B的图像大小均为2000*2000 像素，为TIFF格式。A图像左上角在原始quickbird立体像对左图像上的像素坐标为 (2865,1859), B图像左上角在原始quickbird立体像对右图像上的像素坐标为(2802， 1687)，两者基本为同一区域。由于数据量较大，例子中的显示数据是A图像局部自图像像 素行列坐标(438，226)至(454，242)范围的16X 16大小图像块的处理结果。(1)基于立体像对左影像的边缘提取和细化获得单像素边缘的二值图像；①CANNY算子提取边缘 (438，226)至(454，242)范围的灰度值为22242018201919211816191918203138
21181817182019191917191919253641
17171817181917171920191820324544
18191817171817181919191924394745
18191817171819181819192034454544
18171718181919191918182539464545
18171917191917171920203143464644
18181918181818181820273946464645
17181817181818191919324345454543
18171718191918181825394646464645
18171818161922243040444447474846
19212527343747535859545046454646
37445155636971737068666458484647
66697169677470707070697070605551
69687171707070716968707170706760
69707170717070717169707272707071
选择尺度参_ ο=3，高低阈值分别为hOW75, hhigh=225，进行 CANNY缘提取后。其结果图像如下2552552552552552552552552552552552552552550255
25525525525525525525525525525525525525500255
2552552552552552552552552552552552552550255255
25525525525525525525525525525525525500255255
25525525525525525525525525525525500255255255
2552552552552552552552552552552550255255255255
25525525525525525525525525525500255255255255
25525525525525525525525525525500255255255255
2552552552552552552552552552550255255255255255
25525525525525525525525525500255255255255255
25525525525525525525525500255255255255255255
255255255000002552550255255255255255
000255255255255255255255255000255255
255255255255255255255255255255255255000255
25525525525525525525525525525525525525525500
255255255255255255255255255255255255255255255255其中0值对应像素为边缘位置。②利用Hilditch细化方法，对①边缘图像细化获取单像素边缘二值图，细化图像
像素值如下CN 101901343 A说明书8/28 页2552552552552552552552552552552552552552550255255255255255255255255255255255255255255255025525525525525525525525525525525525525525502552552552552552552552552552552552552552552550255255255255255'255255255255255255255255255025525525525525525 52552552552552552552552552550255255255255255255255255255255255255255255025525525525525525525525525525525525525525525502552552552552552552552552552552552552552550255255255255255255255255255255255255255255255255255255255255255-S 2552552552552552552552552550255255255255255255138,23 1^52552552550000255255255255255.25525525500002552552552552552552550255255255255255255255255255255255'2552552552552550002552552552552552552552552552552552552552552552550255255255255255255255255255255255255255255255255其中0值对应像素为边缘位置。(2)利用改进相位编组方法对上述单像素边缘二值图像矢量化，得到图像边缘线选择的边缘线长度60，首末距离可取60像素。弯曲度阈值可取120度.按照前述 方法得到符合条件的边缘线。(1)中二值图的边缘线长度不满足要求，故舍去。(3)基于灰度相关进行边缘线节点匹配，确定左影像上边缘线节点的同名像点；①仿射变换参数获取在左右片上选取下列同名点，利用最小二乘方法解算一阶仿射模型的系数。
ai-bi的解算结果为-26. 75930. 89760. 0064-193. 94300. 09911. 0277②确定左影像上边缘线节点的同名像点；实验图像上，以一条满足候选边缘线为例。该边缘线上的坐标如下 所在位置的三等分节点，坐标分别为点1(111，485)，点2(127，511)，点3(142， 537)，点 4(157562)。以点1为例左片像素坐标(111，485)，目标窗口的大小为15*15，数据如下: 仿射变换解算出右片的初位置为(75，315)，以此为中心，确定的搜索窗口大小为 19*19，灰度数据如下 计算所得相关系数最大值为0.89，对应的右片影像上同名点坐标为(82312) (4)基于像对RPC模型参数(或外方位元素)前方交会解算节点高程；实验立体像对的RPC参数的为左影像为
IineOffset = 12881 ； (r0)sampOffset = 13793 ； (C0)IatOffset = 31. 2219 ； (X0)IongOffset = 121. 5191 ； (Y0)heightOffset = 20 ； (Z0)IineScale = 13172 ； (rs)sampScale = 13812 ； (Cs)IatScale = 0. 0826 ； (Xs)IongScale = 0. 0944 ； (Ys)heightScale = 501 ； (Zs)IineNumCoef =(
0129]-1. 057032E-03,
0130]-1. 213815E-02,
0131]-9. 912392E-01,
0132]+2. 160953E-02,
0133]-1. 266359E-04,
0134]+6. 745135E-06,
0135]+2. 179383E-04,
0136]-6. 908676E-04,
0137]+1. 798719E-03,
0138]-6. 150191E-06,
0139]-1. 171661E-06,
0140]-2. 743666E-07,
0141]+2. 715213E-05,
0142]-1. 330819E-07,
0143]-1. 143990E-05,
0144]-6. 096293E-05,
0145]-1. 141245E-05,
0146]+4. 867008E-07,
0147]-1. 759363E-05,
0148]+4. 353202E-06)；(以上括号中为⑶式中P1Gn, Yn, Zn)中的系数值 ijk)
0149]IineDenCoef =(
0150]+1. 000000E+00,
0151]-1. 477262E-04,
0152]+6. 431291E-04,
0153]+2. 386261E-04,
0154]+1. 571120E-05,
0155]-6. 717437E-07,
-5. 375222E-06,+8. 588398E-06,-2. 672493E-05,+1. 527395E-05,+3. 376996E-07,+4. 208274E-08,+1. 401940E-05,-8. 380114E-08,+1. 918142E-07,+3. 232610E-04,-6. 729629E-06,-1. 418655E-08,+1. 477277E-05,-1. 099389E-05)；(以上括号中为⑶式中P2 (Xn，Yn, Zn)中的系数值 Cijk)sampNumCoef =(-2. 952752E-03,+1. 001093E+00,+3. 223476E-03,+2. 267689E-03,-1. 820613E-03,+6. 311864E-04,-3. 251605E-04,+1. 043426E-03,-1. 303782E-03,+9. 541229E-06,+2. 433079E-06,-1. 604161E-05,-3. 218350E-05,-7. 563561E-06,+6. 586746E-05,-5. 668938E-05,-2. 200462E-06,+8. 402269E-07,+1. 564317E-05,-2. 139418E-06)；(以上括号中为⑶式中P3(Xn，Yn, Zn)中的系数值ciJk)sampDenCoef =(+1. 000000E+00,+1. 915708E-03,0194]+1. 836343E-03,
0195]-6. 241859E-04,
0196]-4. 831275E-05,
0197]+2. 811915E-06,
0198]-1. 959083E-06,
0199]-1. 618778E-08,
0200]+1. 738714E-05,
0201]-8. 081709E-06,
0202]-1. 511757E-07,
0203]-7. 945617E-08,
0204]+1. 434163E-06,
0205]+0. 000000E+00,
0206]-4. 982303E-08,
0207]-2. 951413E-07,
0208]-7. 620251E-08,
0209]+3. 317219E-08,
0210]+2. 282115E-07,
0211]-1. 476494E-08)；(以上括号中为⑶式中P4(Xn，Yn, Zn)中的系数值 iJk)右影像为(其参数意义同左影像)
0212]IineOffset = 12970 ；
0213]sampOffset = 13762 ；
0214]IatOffset = 31. 2274 ；
0215]IongOffset = 121. 5146 ；
0216]heightOffset = 21 ；
0217]IineScale = 13127 ；
0218]sampScale = 14070 ；
0219]IatScale = 0. 0884 ；
0220]IongScale = 0. 1058 ；
0221]heightScale = 500 ；
0222]IineNumCoef =(
0223]-6. 908239E-03,
0224]+9. 814168E-02,
0225]-1. 087285E+00,
0226]+1. 143412E-02,
0227]+1. 508245E-03,
0228]-4. 152784E-05,
0229]+1. 434624E-04,
0230]-9. 117079E-04,
0231]+6. 039878E-03,
-2. 775409E-06,+8. 759992E-07,-2. 810502E-06,-8. 566692E-05,-6. 499227E-06,+5. 727487E-05,+2. 214905E-04,+7. 206675E-05,-2. 960914E-07,-1. 123729E-05,-1· 175942E-06) ； IineDenCoef =(+1. 000000E+00,+2. 405598E-03,-9. 496480E-04,+2. 331411E-04,-8. 558844E-05,+3. 254521E-07,-7. 707115E-06,-5. 145213E-05,+2. 075941E-04,-6. 519818E-05,-1. 574049E-06,+5. 853576E-08,+1. 562691E-04,-5. 160127E-07,-1. 218015E-05,-6. 104565E-04,+1. 799919E-06,+3. 781986E-08,+1. 033377E-05,-1. 516558E-06)；sampNumCoef =(-6. 521602E-03,+9. 860974E-01,-8. 363866E-04,-2. 074643E-02,-3. 782817E-03,+4. 173930E-04,-2. 525911E-04,
+7.756150E-03,
-7.525981E-04,
-6.612932E-06,
-3.244307E-06,
+3.736748E-05,
+3.411336E-05,
-4.151620E-06,
+3.774801E-05,
-3.128959E-04,
-1.999274E-07,
+6.750285E-06,
+1.305323E-05,
+9.244454E-08)
sampDenCoef=(
+1.000000E+00,
-1.212741E-03,
+3.719373E-03,
-5.473858E-04,
-4.765092E-05,
-3.145637E-06,
+1.865777E-06,
+6.198560E-06,
-4.246335E-05,
-4.569907E-06,
-1.275657E-07,
-2.463548E-07,
+6.169447E-06,
+2.168261E-08,
-1.195460E-06,
-1.163123E-05,
-3.526989E-08,
-9.533090E-08,
+4.465392E-07,
+0.000000E+00)以(3)中候选边缘线所在位置的三等分节点及对应的同名像点为例如1点左片像素坐标为(111，485)，右片像素坐标为(82312)；根据上面的RPC参 数利用前方交会，计算出该点的WGS84系坐标121. 50731. 23916. 982 ；转化为上海城市坐 标3815. 857347. 6677. 182。(5)基于节点高程约束去除虚假道路边缘，进行道路边缘线确定(得到道路提取结果)；对(3)中得到的边缘线列表，依照(4)中介绍的方法对每条线上的三等分进行匹 配和高程求解。约束道路边缘线上的平均高程不超过阈值10米。下面是保留道路边缘线的坐标和高程约束的细节取值例一对上面的边缘线，其三等分点和端点的高程值为7. 1816522. 4432013. 1067649. 34 8106高程均值为5. 5199，不超限，故保留此线段。线长度为78，沿线坐标值如下表所示 例二该实验影像对内，还有一条边缘线线长度为68，沿线坐标值如下表所示 其三等分点和端点的高程值为7. 06537416. 36848117. 94507623. 847347，高程均 值为16. 3066，超过10，因此根据规则，该线段被删除舍去。对实验的左影像中，通过步骤(3)得到所有的道路边缘线都进行(4)、(5)，得到最 终的道路边缘线矢量结果，如附图4。
权利要求
一种基于立体约束的遥感影像道路提取方法，其特征在于包括如下步骤(1)基于立体像对左影像的边缘提取和细化获得单像素边缘的二值图像；(2)利用改进相位编组方法对上述单像素边缘二值图像进行形态约束和矢量化，得到图像候选道路边缘线；(3)基于灰度相关进行候选道路边缘线节点匹配，确定边缘线节点在右影像上的同名像点；(4)基于像对RPC模型参数或内\外方位元素前方交会解算节点高程；(5)基于节点高程约束的道路边缘线确定得到道路提取结果。
2.根据权利要求1所述的基于立体约束的遥感影像道路提取方法，其特征在于步 骤(1)所述的基于立体像对左影像的边缘提取和细化获得单像素边缘的二值图像，基于 CANNY算子的左影像边缘检测，以及在此边缘检测结果上基于Hilditch细化方法的左影像 单像素边缘二值图像生成，具体按如下步骤进行V(G* f)①将图像f与尺度为σ的高斯函数G做卷积G*f，利用式〃=|V(GV)|，对图像的每个像素，计算局部边缘的法向量η ；②使用非最大值抑制方法找到法向量η方向的方向梯度最大值位置(X，y)，= 0实际做法是，沿着点(x，y)处法向量方向n，在(x，y)两侧各列出一定范围点 δη,的方向梯度值，利用这些散点，通过曲线拟合的方法可以确定梯度极大值，U，y)的值可更 新为此位置；③利用式 I VG=Vk力 I=I G*/(x + l，力-GsV(X-Ij)I+ |G*/(x，少+计算边缘点(x，y)的梯度幅值lv(G*/)L ④用高低阈值hlOT，Iihhigh判断该点是否属于边缘，若|V(G*/)|>，贝_定为边缘点； 若|V(G*/)| 且在8邻域与某些高阈值边缘点连接也认为是边缘点，将边缘点像素值置 为1，非边缘点置为0;⑤将CANNY算子提取的边缘图像转化为单像素边缘二值图像，利用Hilditch细化方 法，完成轮廓细化，在步骤④获得的二值图像中，像素点的值为1或0，值为1的点，称为前景 点或边缘，即图形像素，值为0的点，称为背景点，在细化算法中，一个前景点满足下面的条 件才能被删除a.该点为边界点；b.该点非端点；c.该点非断点；d.该点的删除不应该造 成过度侵蚀；Hilditch细化方法需要对整幅图像作多次遍历，在每一次遍历过程中，通过对像素点 八邻域与模板进行匹配，标记同时满足上述条件a d的像素点，在此次遍历结束后，再将 这些标记的像素点删除即由前景点转化为背景点，然后开始下一次遍历，直到在某一次遍 历过程中，没有像素点被标记，整个细化过程结束；经过细化运算后得到单像素边缘二值图 像，边缘像素点的灰度为0，非边缘像素点的灰度为1。
3.根据权利要求2所述的基于立体约束的遥感影像道路提取方法，其特征在于定义步 骤(1)中①中的尺度参数σ =3，步骤 中的高低阈值111 = 75，、_ = 225。
4.根据权利要求1所述的基于立体约束的遥感影像道路提取方法，其特征在于步骤(2)所述的利用改进相位编组方法对上述单像素边缘二值图像进行形态约束和矢量化，得 到图像候选道路边缘线，具体按如下步骤进行设计三个阈值来限定道路边缘，依次为边缘线的长度阈值，首末距离阈值和弯曲程度 阈值，长度阈值指边缘线上像素个数，首末距离阈值指首末两端点的距离，在边缘弯曲甚至 首尾接近时，首末距离阈值可以用来排除非道路边缘，弯曲程度阈值用首末点和线上一结 点的连线夹角表示，若夹角小于此阈值，则说明弯曲过大，应从道路边缘中排除；具体为①将步骤(1)中提取的单像素边缘图像上，所有边缘像素点的状态设为0，同时确定道 路边缘线提取的长度阈值Th-L，首末距离阈值Th-D和弯曲程度阈值Th-curve ；②对单像素边缘图像上，状态为0的边缘像点逐像素判断，计算边缘像素8邻域中灰度 为0的点数m，如m = 1，则确定其为边缘线的一个端点；③从该端点起沿着该线段逐点记录坐标，记录到点列表PointList中，直到到达另一 个端点，长度L为该单像素边缘中包含的像素点个数，首末距离D为首末端点利用欧式距离公式计算得到的结果，D = -χ^)2+(^m-^)2，如果D彡Th-D, L彡Th-L,则保留该 线段；并对该列线段进行平分，找到3分节点，利用式CGS^ = Hf广计算每个分节点IHIII0II(Xi, Yi)和首末点的连线夹角θ i，设向量a= (Xi_X首，y「y首)，b = (Xi-X末，yi-y末)，对每 个分节点和首末点的连线夹角Θ i进行判断，若ei<Th-curve阈值，则认为该线段弯曲 予以去除，如果不是则该线段予以保留，作为候选道路边缘；对跟踪判断过的边缘像素点， 无论是否作为候选道路边缘，都将其标记为“已跟踪判断”即边缘像素点状态改为1 ；④对步骤(1)中所有未跟踪的边缘即状态为0，重复步骤②、③，直到所有边缘线都完 成上述判断，最后得到候选的矢量化的道路边缘线，这些候选边缘线放入线列表IineList 中。
5.根据权利要求4所述的基于立体约束的遥感影像道路提取方法，其特征在于定义步 骤(2)中①中的边缘线长度阈值Th-L为18/P到48/P像素范围，P为该卫星图像的像素分 辨率；首末距离阈值Th-D为18/P到36/P像素范围；弯曲阈值为大于等于120度小于180度。
6.根据权利要求1所述的基于立体约束的遥感影像道路提取方法，其特征在于步骤(3)所述的基于灰度相关进行边缘线节点匹配，确定边缘线节点在右影像上的同名像点，具 体按如下步骤进行对步骤(2)中提取的候选道路边缘线进行三等分，加上端点共计四个分节点，利用基 于相关系数的灰度相关在右片上搜索出相应的同名点，为增强匹配的可靠性和速度，通过 仿射变换模型确定匹配的初始位置，初始位置的确定方法如下①采用一阶仿射变换模型来确定待匹配点在右片影像的同名点初始位置，— + ci ^pcγ + Cl >2 y^Y<y2 =K+bXxI+bIyi选取若干左右影像同名像点3个以上，利用最小二乘方法解算以上6个参数(a” bi)； ②依据这些系数求解匹配点的初始位置，利用式<1以及左影 像上待匹配点的坐标(Xl，yi)，计算出其在右影像上同名点的粗略位置(4，/)；③确定右片影像上粗略位置后，在以此位置为中心的19*19的搜索窗口 中寻找左影像点的匹配点的精确位置，左片匹配窗口大小为15*15，依据式 _计算相关系数，并选取搜索窗口内相关系数最大的点作 为匹配最佳位置，即可得到精确的同名点位置(^，y》，式中X，y分别为匹配窗口和右片待 相关窗口内的像素灰度值J和^分别是两窗口内的灰度均值，N为窗口中的像素个数；利用②、③确定每条候选道路边缘线上4个分节点在右影像上的同名点坐标，并记录 下来。
7.根据权利要求1所述的基于立体约束的遥感影像道路提取方法，其特征在于步骤 (4)所述的基于像对RPC模型参数或内\外方位元素前方交会解算节点高程，具体按如下步 骤进行有理多项式函数 其中 0 ^ Ii1 ^ order，0 ^ n2 ^ order，0 ^ n3 ^ order 且 Wn3 ^ order, order 为多 项式的阶数，一般对于线阵推扫式影像，取order = 3，已可满足精度的要求，这样每个P (X， Y，Z)是一个20项的三阶多项式； 代入上述有理多项式函数得到下式 立体像对中每幅影像的有理多项式函数RPC中的Cijk已知或内、外方位元素已知，平移 参数IV C0, X0, \、Z0和比例参数rs、Cs, Xs、Ys、Zs也已知时，根据从步骤(3)中得到的匹配 同名点对，利用前方交会原理，每对同名点可以列出4个等式，解算出该点对应地面物点的 三维坐标(X，Y，Z)，其中Z坐标一般对应该物点的高程信息；一条边缘线三等分点可以得到四个分节点的高程坐标，将其高程坐标平均，得到
8.根据权利要求1所述的基于立体约束的遥感影像道路提取方法，其特征在于步骤 (5)所述的基于节点高程约束的道路边缘线确定得到道路提取结果，具体按如下步骤进 行根据所处理影像区域地面平均高程情况，给定道路高程约束阈值TH-h，如平均高程大 于道路高程阈值Haraage > TH-h，认为该四点所在边缘是非道路边缘如建筑边缘，保存满足 高程约束条件的道路边缘线在IineList列表中，对每条候选道路边缘线逐个按照步骤(4) 中的方法计算候选道路边缘线高程坐标平均值Harera一并判断其是否是道路边缘线，对满足 要求的道路边缘线矢量数据保留存在列表RoadList中，并记录输出。
9.根据权利要求8所述的基于立体约束的遥感影像道路提取方法，其特征在于定义步 骤(5)中的高程约束阈值TH-h为该地区平均高程值加上7-8m。
全文摘要
本发明涉及一种基于立体约束的遥感影像道路提取方法，包括基于立体像对左影像的边缘提取和细化获得单像素边缘的二值图像；利用改进相位编组方法对单像素边缘二值图像进行形态约束与矢量化；基于灰度相关进行边缘线节点匹配，确定边缘线节点在右影像上的同名像点；基于像对RPC模型参数或内\外方位元素前方交会解算节点高程；基于节点高程约束的道路边缘线确定得到道路提取结果。本发明从立体影像入手，对边缘提取、细化、经过形态约束和矢量化后的道路边缘线的高程信息进行提取，建立对候选道路边缘线的分节点进行匹配、根据同名点信息进行前方交会获得高程信息的模型，根据道路的高程特点，建立高程约束条件，进行高程约束，从候选道路边缘线中提取满足要求的道路边缘线，从而实现遥感影像上的道路专题信息的快速自动提取。
文档编号G06K9/46GK101901343SQ20101023325
公开日2010年12月1日 申请日期2010年7月20日 优先权日2010年7月20日
发明者叶勤, 张小虎, 王卫安 申请人:同济大学