本发明公开了一种多层次知识驱动的全色遥感影像的道路变化信息提取方法，首先在道路提取层中充分应用先验知识，自动生成多尺度模板，实现道路候选段的自动提取，根据人对道路模型的感知知识，采用感知编组连接候选道路段，形成初始道路网；其次在变化道路检测层中加入知识判断原则，对提取道路网和旧道路网进行缓冲区分析，得到道路变化检测结果；再次在新增道路检测层中，结合道路网模型进行推理和假设，得到符合道路模型的知识约束规则，从视觉分割得到分割块反映影像属性的角度出发，对候选新增道路进行优化和处理；最后在半自动提取层中针对道路复杂性，在人的知识的驱动下完成新增道路的半自动跟踪，产生完整的新增道路网。
多层次知识驱动的全色遥感影像的道路变化信息提取方法
技术领域
[0001] 本发明涉及一种多层次知识驱动的全色遥感影像的道路变化信息提取方法，属于影像处理领域，特别是属于遥感影像处理与目标提取技术领域。
背景技术
[0002] 道路作为最重要的人工地理实体，在地理基础信息中占有极大的比重。从遥感影像中自动发现道路变化不仅是摄影测量与遥感领域的难题，也是计算机视觉与图像理解研究的重点之一。
[0003] 在过去的二十多年中，人们提出许多从遥感图像中自动或半自动提取道路变化信息的方法，大致可以分为三类：第一类方法是基于道路本身特性，建立合适的数学模型来描述理想道路，然后运用数学方法求解，典型算法包括模板匹配法(McKeown，Denlinger.1988，Vosselman，Knecht.1995)、平行边跟踪连接法(Nevatia，Babu.1980)、动态规划法(Gruen，Li.1995)、多尺度分析法(Heipke etal.1996，Mayer，Steger.1998)等等；第二类方法是基于上下文信息，分析道路与其它目标之间的相互关系检测道路变化，典型算法包括使用知识库自动调整参数的方法(Trinder.1998)、利用上下文信息(车辆、树木等)支持或已有地图引导对道路做出判断和有效补充的方法(Stilla 1995，Baumartner et al.1999)等等；第三类方法是前两者的混合，是道路变化检测研究的趋势，典型算法包括模拟人的感知，以段的几何与灰度属性为主，顾及局部上下文线索的感知编组算法(胡翔云，2001)、基于线特征检测道路段，通过道路网的全局约束条件进行道路网连接实现变化道路检测(钟家强，2007)等等。
[0004] 目前还没有一种针对所有道路类型和各种分辨率图像都通用的道路自动变化检测策略和算法。
发明内容
[0005] 综合利用现有技术并克服其局限性，本发明提供了一种多层次知识驱动的全色遥感影像的道路变化信息提取方法，可以快速有效地检测出道路的变化，将检测结果提供给用户，提高地图修测的自动化程度。
[0006] 实现本发明目的采用的技术方案是：一种多层次知识驱动的全色遥感影像的道路变化信息提取方法，首先在道路提取层中充分应用先验知识，自动生成多尺度模板，实现道路候选段的自动提取，根据人对道路模型的感知知识，采用感知编组连接候选道路段，形成初始道路网；其次在变化道路检测层中加入知识判断原则，对提取道路网和旧道路网进行缓冲区分析，得到道路变化检测结果；再次在新增道路检测层中，结合道路网模型进行推理和假设，得到符合道路模型的知识约束规则，从视觉分割得到分割块反映影像属性的角度出发，对候选新增道路进行优化和处理；最后在半自动提取层中针对道路复杂性，在人的知识的驱动下完成新增道路的半自动跟踪，产生完整的新增道路网。其包括以下具体步骤：
[0007] (1)对近期全色遥感影像进行预处理，预处理操作包括：(a)对全色遥感影像的正射纠正，(b)对正射影像的拼接，(c)根据同区域旧矢量数据范围，对拼接的正射影像进行裁减，(d)将裁减结果作为变化检测的输入影像；
[0008] (2)对输入影像进行道路增强和噪声去除，操作包括：(a)对输入影像进行Wallis滤波加大图像反差，突出道路信息，(b)高斯平滑去除噪声；
[0009] (3)在先验知识的引导下，生成道路候选段，操作包括：(a)由已有旧道路矢量地图获得检测区域道路网的先验知识，包括道路的最大路宽、最小路宽，自动生成多尺度模板，在影像的行、列方向分别计算模板相关系数，以相关系数峰值来检测影像断面信号的矩形脉冲，得到候选的道路中心点，生成道路中心点的二值图像；(b)对二值图像进行线条细化得到中心骨架；(c)对中心骨架进行单纯细化，去除交叉点和分支点；(d)对去除交叉点和分支点的结果图像进行数学形态学处理，填补段之间的缝隙，得到连通的道路线；(e)对连通的道路线进行矢量跟踪，删除短线；
[0010] (4)根据人对道路模型的感知知识，采用感知编组连接候选道路段，形成道路网，包括以下步骤：
[0011] (4-a)根据人对道路模型的感知知识，即道路几何属性、灰度属性、道路与其局部上下文关系描述，以道路几何与灰度属性为主，顾及局部上下文线索来进行共线链的编组，将各段连接成无分支的共线长链；
[0012] (4-b)提取交叉和分支形成道路网络；
[0013] (5)在缓冲区检测中加入知识判断原则，结合已有旧道路矢量数据，检测道路网中变化和消失路段，具体检测步骤如下：
[0014] (5-a)旧道路矢量地图与提取道路网配准；
[0015] (5-b)确定检测缓冲区的缓冲距离。忽略数据转换、计算等的误差，旧矢量地图与提取道路网的匹配误差由图形和图像配准精度、图形误差和图像误差决定，由协方差传播率可获得检测缓冲距离值：
[0016]
[0017] 其中，BufferDis为缓冲距离值，σmap为地图成图中误差，σreg为配准精度，σroad为道路平均宽度；
[0018] (5-c)在已有旧道路矢量两侧，用缓冲距离构造缓冲区，将缓冲区与提取道路网叠加，计算道路变化率：
[0019]
[0020] 其中，δ为道路变化率，lroad为缓冲区内道路网长度，ltotal为道路矢量长度，根据判断阈值T，若δ大于T则认为该段道路未变化，置不变标志，否则设置可能变化标志；
[0021] (5-d)加入知识判断原则，包括：
[0022] (5-d-1)道路是平滑的带状曲线，曲率变化平缓；
[0023] (5-d-2)道路附近地物复杂，导致某段因遮盖找不到道路点；
[0024] 针对以上两个判断原则，对设置可能变化标志的道路段进一步判断，判断步骤如下：
[0025] (5-d-a)以该路段走向为中心，根据角度缓冲值，设定正向角度搜索区间和反向角度搜索区间：
[0026] BufferAngle1＝[α-δ，α+δ]
[0027] BufferAngle2＝[π+α-δ，π+α+δ]
[0028] 其中，BufferAngle1为正向角度搜索区间，BufferAngle2为反向角度搜索区间，α为道路走向与水平方向的夹角，δ为角度缓冲阈值；(5-d-b)在旧矢量地图上，延路段走向正向，在正向角度搜索区间内搜索具有不变标志的道路，将距离最近的路段作为当前路段的前路段R1，若未找到前路段，则R1赋无结果标志；延路段走向反向，在反向角度搜索区间内搜索具有不变标志的道路，将距离最近的路段作为当前路段的后路段R2，若未找到后路段，则R2赋无结果标志；(5-d-c)若该路段为首路段，后道路R2被赋予无结果标志，则为该路段赋消失标志；若该路段为尾路段，前道路R1被赋予无结果标志，则为该路段赋消失标志；若该路段为中间路段，前道路R1或后道路R2被赋予无结果标志，则为该路段赋消失标志，进入(5-e)，否则，计算参数：
[0029]
[0030] 其中，P为存在测度，|γR1-γR2|为R1和R2的曲率差，dR1R2为R1和R2临近端点之间的距离，Tγ为曲率差阈值，Td为距离阈值，ωγ和ωd分别为曲率差和距离的权。若P小于1，则该路段没有发生变化，置不变标志；否则，该路段变化，置变化标志；
[0031] (5-e)根据无变化标志，删除道路网中对应路段；根据变化标志，将变化的旧道路矢量段输出；
[0032] (6)结合道路网模型进行推理和假设，得到符合道路模型的知识约束规则，从视觉分割得到分割块反映影像属性的角度出发，对候选新增道路进行优化和处理，处理步骤如下：
[0033] (6-a)对步骤(2)得到的图像进行视觉分割；
[0034] (6-b)从道路网模型的角度对影像分割块属性进行合理推理，得到道路模型的知识约束规则，提出四个假设：
[0035] (6-b-1)分割块的属性唯一，即分割得到的区域道路属性值为true或false：
[0036] R(segi)＝true or false
[0037] segi为第i个分割区域，R(segi)为区域道路属性；
[0038] (6-b-2)非道路区域内部属性相同，即非道路区域不包含道路段：
[0039]
[0040] p为影像像素；
[0041] (6-b-3)一分割块只包含一条道路，若道路区域的分割块内存在多条道路，则道路区域上提取出无法连接的道路中线：
[0042]
[0043]
[0044] road为该区域内唯一道路；
[0045] (6-b-4)道路分割块内道路线反映块的骨架，即道路区域的走向应与道路中心线的走向一致，道路中心线必须能够代表道路区域与道路功能有关的几何特征：
[0046]
[0047] “≈”表示道路必须块的骨架形态；
[0048] (6-c)从视觉分割结果和四个假设，对候选新增路段进行优化和处理，步骤如下：
[0049] (6-c-a)利用分割块宽度和由旧矢量地图得到的道路宽度范围比较，判断分割块道路属性；
[0050] (6-c-b)判断移除非道路区域内短线。由假设(6-b-2)，判断道路段R上的点p所在分割块seg，若R(seg)＝false，则该路段位于非道路区域，从道路网内移除该点。
[0051] (6-c-c)结合假设(6-b-3)和(6-b-4)，处理道路区域内多条道路的情况，具体步骤如下：
[0052] (6-c-c-1)道路网为roadnet，输入满足p∈roadnet的点p，判断所在分割块seg，并判断所在roadi∈{road1，...，roadn}，记录seg包含的roadi序号；
[0053] (6-c-c-2)判断R(seg)，如果R(seg)＝false，则从道路网上移除该点，并移除seg包含的roadi序号记录；
[0054] (6-c-c-3)对于道路网上下一点p，转到(6-c-c-2)，直到处理所有道路点；
[0055] (6-c-c-4)对记录了包含道路序号的分割块seg，判断块内道路数量；
[0056] (6-c-c-5)如果seg包含道路数量大于1，则从道路网上移除seg内包含的较短道路部分；
[0057] (6-c-d)利用道路区域的延伸性与道路分割块的形态，处理道路区域内道路段的延伸，具体步骤如下；
[0058] (6-c-d-1)输入分割块骨架线Tseg及块seg内道路段road，计算Tseg的端点到道路线段road的距离，设定较小值W，记录距离大于W的端点p，p∈Tseg；
[0059] (6-c-d-2)从点p出发，计算Tseg中下一点到道路线段road的距离，判断距离是否较大，即判断距离是否大于步骤(6-c-d-1)选择的较小值W；
[0060] (6-c-d-3)如果距离较大，则转到步骤(6-c-d-2)，依次处理Tseg下一点。如果距离不大，则中止步骤(6-c-d-2)，转到步骤(6-c-d-4)；
[0061] (6-c-d-4)将Tseg中作过判断的点作为道路点加入到道路网中。采用最近点的连接方法连接首个距离不大的点与道路线段road，即连接该点与该点到road最短距离点；
[0062] (7)在人的知识的驱动下，完成新增道路的半自动提取，步骤如下：
[0063] (7-a)提供点选删除和拉框删除接口，删除错提道路线；
[0064] (7-b)人工给出少量种子点作为知识，驱动半自动跟踪过程，利用基于最小二乘B样条的主动轮廓模型对新增道路进行半自动提取，具体步骤如下：
[0065] (7-b-1)利用多尺度模板匹配的方法，自动设置基于最小二乘B样条的主动轮廓模型参数；
[0066] (7-b-2)自适应模板匹配自动产生新的道路点；
[0067] (7-b-3)采用人工选取和(7-b-2)自适应模板生成的道路点作为种子点，利用基于最小二乘B样条的主动轮廓模型对道路进行提取和拟合；
[0068] (7-c)自动修改拓扑关系，生成新的道路网；
[0069] (7-d)输出道路变化检测结果。
[0070] 本发明充分利用待修测地区旧道路矢量数据对全色遥感影像进行上述操作，可以快速有效地检测出道路的变化，将检测结果提供给用户，并提供了有效的人机交互接口，便于操作人员对检测结果进行实时监测和修改，为道路变化检测提供了有效的方法。
附图说明
[0071] 下面结合附图和具体实施方式对本发明作进一步的说明。
[0072] 图1为本发明流程示意图。
[0073] 图2是某地区IKNOS全色影像经过预处理得到的输入影像的局部影像。
[0074] 图3是图2与矢量数据叠加结果，黑色线为该地区道路矢量。
[0075] 图4是图2所示地区的初始道路网提取结果。
[0076] 图5是图4所示初始道路网提取结果与旧矢量数据进行缓冲区分析和知识判断结果。
[0077] 图6是图2的视觉分割结果。
[0078] 图7是图2所示地区新增道路采用视觉分割优化结果。
[0079] 图8是另一地区全色遥感影像的局部影像。
[0080] 图9是图8影像道路种子点选取图。
[0081] 图10是在图9影像上基于最小二乘B样条主动轮廓模型自动跟踪的道路与图8叠加结果。
具体实施方式
[0082] 本实施例首先通过处理旧矢量数据和经过预处理的遥感影像，在先验知识的指导下提取候选道路段，根据人对道路模型的感知知识，采用感知编组连接候选道路段，形成初始道路网，其次在对提取道路网和旧道路网进行缓冲区分析的过程中，加入知识判断原则，获得消失和变化道路，再次结合道路网模型进行推理和假设，得到符合道路模型的知识约束规则，从视觉分割得到分割块反映影像属性的角度出发，对候选新增道路进行优化和处理，最后由人提供知识，驱动半自动跟踪过程，利用改进的基于最小二乘B样条主动轮廓模型跟踪道路，对新增道路网进行补充，完成变化检测。其流程如图1所示。
[0083] 本实施例选取某地区IKNOS全色影像和另一地区遥感影像进行实验，选取局部典型区域具体说明本发明的实施步骤如下：
[0084] (1)对近期全色遥感影像进行预处理，预处理操作包括：
[0085] (1-a)对全色航空影像的正射纠正；
[0086] (1-b)对正射影像的拼接；
[0087] (1-c)根据同区域旧矢量数据范围，对拼接的正射影像进行裁减；
[0088] (1-d)将裁减结果作为变化检测的输入影像，图2为某地区IKNOS全色影像经过预处理所得输入影像的局部典型区域影像，图8为另一地区遥感影像经过预处理后所得输入影像的局部典型区域影像；
[0089] (2)对输入影像进行道路增强和噪声去除，操作包括：
[0090] (2-a)对输入影像进行Wallis滤波加大图像反差，突出道路信息，[0091] (2-b)高斯平滑去除噪声；
[0092] (3)在先验知识的引导下，生成道路候选段，操作包括：
[0093] (3-a)由已有旧道路矢量地图获得检测区域道路网的先验知识，包括道路的最大路宽Wmax、最小路宽Wmin，自动生成多尺度模板，如下：
[0094] Rwidth＝{W|Wminx≤W≤Wmax，W为正整数}
[0095]
[0096] 其中，Rwidth为道路宽度集合，Modeli为第i个多尺度模板，其宽度为Rwi+4，高度为3；
[0097] (3-b)在影像的行、列方向分别计算模板相关系数，以相关系数峰值来检测影像断面信号的矩形脉冲，得到候选的道路中心点，生成道路中心点的二值图像；
[0098] (3-c)对二值图像进行线条细化得到中心骨架；
[0099] (3-d)对中心骨架进行单纯细化，去除交叉点和分支点；
[0100] (3-e)对去除交叉点和分支点的结果图像进行数学形态学处理，填补段之间的缝隙，得到连通的道路线；
[0101] (3-f)对连通的道路线进行矢量跟踪，删除短线；
[0102] (4)根据人对道路模型的感知知识，采用感知编组连接候选道路段，包括以下步骤：
[0103] (4-a)根据人对道路模型的感知知识，即道路几何属性、灰度属性、道路与其局部上下文关系描述，以道路几何与灰度属性为主顾及局部上下文线索来进行共线链的编组，将各段连接成无分支的共线长链；
[0104] (4-b)提取交叉和分支形成道路网络，如图4中，黑色线表示图2感知编组连接候选道路段后自动提取的道路网；
[0105] (5)在缓冲区检测中加入知识判断原则，结合已有旧道路矢量数据，检测道路网中变化和消失路段，具体检测步骤如下：
[0106] (5-a)旧道路矢量地图与提取道路网配准，旧矢量数据与图2叠加结果如图3所示，其中黑色线表示旧矢量图上的道路；
[0107] (5-b)确定检测缓冲区的缓冲距离。忽略数据转换、计算等的误差，旧矢量地图与提取道路网的匹配误差由图形和图像配准精度、图形误差和图像误差决定，由协方差传播率可获得检测缓冲距离值：
[0108]
[0109] 其中，BufferDis为缓冲距离值，σmap为地图成图中误差，σreg为配准精度，σroad为道路平均宽度。图6与矢量地图的配准精度为1.2像素，矢量地图的成图精度为5像素，影像中道路的宽度平均值为13.6像素。根据缓冲区距离的计算公式，得到缓冲区距离为14.5像素；(5-c)在已有旧道路矢量两侧，用缓冲距离构造缓冲区，将缓冲区与提取道路网叠加，计算道路变化率：
[0110]
[0111] 其中，δ为道路变化率，lroad为缓冲区内道路网长度，ltotal为道路矢量长度，根据判断阈值T，若δ大于T则认为该段道路未变化，置不变标志，否则设置可能变化标志；
[0112] (5-d)加入知识判断原则，包括：
[0113] (5-d-1)道路是平滑的带状曲线，曲率变化平缓；
[0114] (5-d-2)道路附近地物复杂，导致某段因遮盖找不到道路点；
[0115] 针对以上两个判断原则，对设置可能变化标志的道路段进一步判断，判断步骤如下：
[0116] (5-d-a)以该路段走向为中心，根据角度缓冲值，设定正向角度搜索区间和反向角度搜索区间：
[0117] BufferAngle1＝[α-δ，α+δ]
[0118] BufferAngle2＝[π+α-δ，π+α+δ]
[0119] 其中，BufferAngle1为正向角度搜索区间，BufferAngle2为反向角度搜索区间，α为道路走向与水平方向的夹角，δ为角度缓冲阈值；
[0120] (5-d-b)在旧矢量地图上，延路段走向正向，在正向角度搜索区间内搜索具有不变标志的道路，将距离最近的路段作为当前路段的前路段R1，若未找到前路段，则R1赋无结果标志；延路段走向反向，在反向角度搜索区间内搜索具有不变标志的道路，将距离最近的路段作为当前路段的后路段R2，若未找到后路段，则R2赋无结果标志；
[0121] (5-d-c)若该路段为首路段，后道路R2被赋予无结果标志，则为该路段赋消失标志；若该路段为尾路段，前道路R1被赋予无结果标志，则为该路段赋消失标志；若该路段为中间路段，前道路R1或后道路R2被赋予无结果标志，则为该路段赋消失标志，进入(5-e)，否则，计算参数：
[0122]
[0123] 其中，P为存在测度，|γR1-γR2|为R1和R2的曲率差，dR1R2为R1和R2临近端点之间的距离，Tγ为曲率差阈值，Td为距离阈值，ωγ和ωd分别为曲率差和距离的权。若P小于1，则该路段没有发生变化，置不变标志；否则，该路段变化，置变化标志；
[0124] (5-e)根据无变化标志，删除道路网中对应路段；根据变化标志，将道路变化检测结果输出。图5为图2区域与旧矢量数据缓冲区分析和知识判断结果，线a表示变化、消失的道路段，其余黑线表示新增的道路段；
[0125] (6)道路网中剩余路段为候选新增路段，结合道路网模型进行推理和假设，得到符合道路模型的知识约束规则，从视觉分割得到分割块反映影像属性的角度出发，对候选新增道路进行优化和处理。由于道路的复杂性，存在三类错误：第一类，非道路区域内检测到道路中心线；第二类，道路区域上提取出无连接的道路中心线；第三类，道路中断于存在道路延伸的区域。从视觉分割得到分割块反映影像属性的角度，结合高层知识，引入符合道路网模型和推理的假设，对候选新增路段进行优化和处理，处理步骤如下：
[0126] (6-a)对步骤(2)得到的图像进行视觉分割，G＝(V，E)是无向图，V是有待分割的元素集，连接相邻点对构成边(vi，vj)∈E。每条边(vi，vj)∈E都有相应的权值w((vi，vj))，这个权值衡量相邻元素vi和vj间的不同，对于灰度图像，图像上的每个像素有相应的点vi∈V对应，E中的每条边端点为处于八连通区域中的两个相邻像素。以边所连接的两个像素间的绝对强度不同为基础，定义的边权值函数如下：
[0127] w((vi，vj))＝|I(pi)-I(pj)|
[0128] 其中I(pi)是像素pi的灰度值。为了评价在分割得到的两个部分间是否有明显的边界，以边界上元素的相异性为基础上对判断的标准进行了定义。定义如下：
[0129] (6-a-1)定义区域C的内部差异为在这个部分的中的最大权值。
[0130]
[0131] (6-a-2)定义两个部分 “之间的差异”为连接这两个部分的边的最小权值。
[0132]
[0133] 如果没有边连接C1和C2，Dif(C1，C2)＝∞。
[0134] (6-a-3)得到两个部分间是否明有明显的边界的判断标准D。
[0135]
[0136] 其中，最小内部差异MInt如下定义：
[0137] MInt(C1，C2)＝min(Int(C1)+r(C1)，Int(C2)+r(C2))
[0138] γ(C)＝k/|C|
[0139] |C|表示C中元素的个数。k为运行时确定的参数，受图像分辨率和图像细节的影响，对分割结果而言相当于设置了一个观察范围，较大的k值会使分割得到的分割块更大。具体步骤如下：
[0140] 输入图形G’＝(V，E’)，它有n个点，如图2，n为图像像素总数；m条边，输出为分割S，即将V划分为部分S＝(C1，C2，...，Cr)，步骤如下：
[0141] (6-a-a)按照边缘的权重由小到大的顺序，将E排序，得到(e1，...，em)；
[0142] (6-a-b)初始化分割为S0，在这个分割里面每一个点vi都作为一个独立分割块；
[0143] (6-a-c)从q＝1，到m，重复(6-a-d)；
[0144] (6-a-d)由Sq-1构造Sq。vi和vj为序列中(e1，...，em)的第q条边的两点，即eq＝q-1 q-1 q-1 q-1(vi，vj)。Ci 为S 中包含vi的一个部分，而Cj 是S 中包含vj的部分。如果并且 那么Sq可以从Sq-1中通过融合Ciq-1、Cjq-1而获得，否则Sq＝
q-1
S ；
m
[0145] (6-a-e)返回S＝S，对图2图像分割，得到结果如图6；
[0146] (6-b)从道路网模型的角度对影像分割块属性进行合理推理，得到道路模型的知识约束规则，提出四个假设：
[0147] (6-b-1)分割块的属性唯一，即分割得到的区域道路属性值为true或false：
[0148] R(segi)＝true or false
[0149] segi为第i个分割区域，R(segi)为区域道路属性；
[0150] (6-b-2)非道路区域内部属性相同，即非道路区域不包含道路段：
[0151]
[0152] p为影像像素；
[0153] (6-b-3)一分割块只包含一条道路，若道路区域的分割块内存在多条道路，则道路区域上提取出无法连接的道路中线：
[0154]
[0155]
[0156] road为该区域内唯一道路；
[0157] (6-b-4)道路分割块内道路线反映块的骨架，即道路区域的走向应与道路中心线的走向一致，道路中心线必须能够代表道路区域与道路功能有关的几何特征：
[0158]
[0159] “≈”表示道路必须块的骨架形态；
[0160] (6-c)从视觉分割结果和四个假设，对候选新增路段进行优化和处理，步骤如下：
[0161] (6-c-a)利用分割块宽度和由旧矢量地图得到的道路宽度范围比较，判断分割块道路属性；
[0162] (6-c-b)判断移除非道路区域内短线。由假设(6-b-2)，判断道路段R上的点p所在分割块seg，若R(seg)＝false，则该路段位于非道路区域，从道路网内移除该点。
[0163] (6-c-c)结合假设(6-b-3)和(6-b-4)，处理道路区域内多条道路的情况，具体步骤如下：
[0164] (6-c-c-1)道路网为roadnet，输入满足p∈roadnet的点p，判断所在分割块seg，并判断所在roadi∈{road1，...，roadn}，记录seg包含的roadi序号；
[0165] (6-c-c-2)判断R(seg)，如果R(seg)＝false，则从道路网上移除该点，并移除seg包含的roadi序号记录；
[0166] (6-c-c-3)对于道路网上下一点p，转到(6-c-c-2)，直到处理所有道路点；
[0167] (6-c-c-4)对记录了包含道路序号的分割块seg，判断块内道路数量；
[0168] (6-c-c-5)如果seg包含道路数量大于1，则从道路网上移除seg内包含的较短道路部分；
[0169] (6-c-d)利用道路区域的延伸性与道路分割块的形态，处理道路区域内道路段的延伸，具体步骤如下：
[0170] (6-c-d-1)输入分割块骨架线Tseg及块seg内道路段road，计算Tseg的端点到道路线段road的距离，设定较小值W，不妨采用距离道路最大宽度作为该较小值，记录距离大于W的端点p，p∈Tseg；
[0171] (6-c-d-2)从点p出发，计算Tseg中下一点到道路线段road的距离，判断距离是否较大，即判断距离是否大于步骤(6-c-d-1)选择的较小值W；
[0172] (6-c-d-3)如果距离较大，则转到步骤2，依次处理Tseg下一点。如果距离不大，则中止步骤(6-c-d-2)，转到步骤(6-c-d-4)；
[0173] (6-c-d-4)将Tseg中作过判断的点作为道路点加入到道路网中。采用最近点的连接方法连接首个距离不大的点与道路线段road，即连接该点与该点到road最短距离点；
[0174] 对图5中新增道路线，结合视觉分割进行优化和处理，得到准确的道路网，即为图7中黑色道路网；
[0175] (7)由于道路的复杂性，在人的知识的驱动下，对道路进行半自动提取，解决新增道路提取中存在的问题。自动提取的新增道路存在两方面问题，第一，所提取道路并非新增道路，而是几何和辐射特性类似的现状地物；第二，新增道路由于遮挡等干扰未提取出，如图8所示局部典型区域，该段道路在上述操作中没有被检测出。本发明提供人机交互的接口，针对第一类错误提供道路删除操作；针对第二类错误采用改进的基于最小二乘B样条的主动轮廓模型，由人提供知识，驱动半自动跟踪过程，完成新增道路的半自动提取，步骤如下：
[0176] (7-a)提供点选删除和拉框删除接口，删除错提道路线；
[0177] (7-b)人工给出少量种子点作为知识，驱动半自动跟踪过程，利用基于最小二乘B样条的主动轮廓模型对新增道路进行半自动提取，图9为图8区域人工选取种子点结果，其中人工给定的种子点用黑色十字丝表示，具体步骤如下：
[0178] (7-b-1)利用多尺度模板匹配的方法，自动设置基于最小二乘B样条的主动轮廓模型参数，具体步骤如下：
[0179] (7-b-1-1)基于影像，将人工输入的种子点连接成线；
[0180] (7-b-1-2)沿着垂直于线的方向，用多尺度模板匹配的方法进行搜索，取得道路不同地方的宽度W1和相关系数Ti，计算路段平均宽度W和道路类型T：
[0181]
[0182]
[0183] 其中，T1为正相关系数之和，T2为负相关系数之和，T取值为1表示该道路为亮道路，T取值为-1表示该道路为暗道路；
[0184] (7-b-2)自适应模板匹配自动产生新的道路点；
[0185] (7-b-3)采用人工选取和(7-b-2)自适应模板生成的道路点作为种子点，利用基于最小二乘B样条的主动轮廓模型对道路进行提取和拟合，结果如图10所示，曲线A表示跟踪得到道路中心线；
[0186] (7-c)自动修改拓扑关系，生成新的道路网；
[0187] (7-d)输出道路变化检测结果。
