1 引 言

遥感影像水体提取在水资源调查和洪涝灾害监测中发挥着重要作用，目前有水体指数法^[1]、谱间关系法^[2]和面向对象分类^[3]等多种方法，但是要获取准确而连续的水体边界依然需要进一步的处理。文献^[4]提出的主动轮廓搜索模型(Snake)能直接获取目标连续性轮廓，在道路提取^[5]和GIS矢量边界更新^[6]中得到了成熟应用。文献^[7]利用Snake模型提取了双线型河流，文献^[8]利用掩膜Snake模型获取了密云水库面积。但Snake模型对凹边界提取效果不佳，对复杂拓扑结构的处理能力较弱。文献^[9]提出的GVF Snake模型改进了Snake模型的外部力，对浅凹边界具有更强的收敛性能。文献^{[10, 11]}分别利用GVF Snake模型提取了多光谱影像和SAR影像中的海陆边界，文献^[12]研究了GVF Snake的初始轮廓对边界提取结果的影响，并提取了遥感影像中不规则水域边界。但该模型同样需手工获取靠近边缘的初始轮廓，并且也存在对噪声敏感、无拓扑变换能力等不足之处，难以提取有深凹边界和复杂拓扑结构(如包含河中岛)的复杂河流边界。

文献^[13]提出的拓扑自适应动态轮廓(topology adaptive Snake，T-Snake)模型是对Snake模型的另一种改进，用三角网格划分图像，通过判断和处理三角网格与Snake曲线的空间关系来实现拓扑变形，成功应用在医学影像边缘提取中。针对该模型需计算节点移动方向和位置的缺点，文献^[14]用正交网格代替三角网格，省略对节点移动方向的判断，在大脑皮层重建和细胞提取中取得了良好效果。文献^{[15, 16, 17]}进一步将节点约束在网格顶点，成功提取了血管影像边缘。文献^[18]将T-Snake模型与同态滤波技术相结合，文献^[19]改进了T-Snake模型的外部能量，这些改进均降低了图像斑点噪声对T-Snake边缘提取精度的影响。文献^[20]在T-Snake模型中增加了自适应性强的膨胀力和区域能量，使之有效提取图像的弱边界和不连续边界。但这些文献均未考虑Snake曲线在遇到目标内部岛状空洞时的拓扑冲突问题，而这一点却是遥感影像河中岛边界提取时需要考虑的重点。另外，这些研究均采取手动方式进行T-Snake模型初始化，自动化程度不高。

正交T-Snake模型在目标边缘处能够自适应地改变轮廓曲线的拓扑结构，对于有复杂拓扑结构边界的提取能力较强，目前尚未有学者将其应用在遥感影像地物提取中。本文以复杂水体边界自动提取为目标，引入正交T-Snake模型并加以改进，针对包含河中岛的河流边界提取会产生的拓扑冲突，提出了一种检测和处理机制，并依据遥感影像分形特征，提出了水体初始轮廓自动生成方法。

正交T-Snake模型采用一定大小的正交网格划分影像，并用顺序相连的网格顶点所构成的闭合轮廓来表示某一时刻的Snake曲线，如图 1(a)所示，实心点为曲线的内部点，空心点为曲线的节点和外部点，箭头表示节点方向。正交T-Snake模型中对节点有两个约束条件：一是节点必须沿网格线移动，二是节点必须位于网格顶点处。因此，曲线节点总是沿网格线从一个网格顶点移动到下一个网格顶点。

图 1 Snake曲线的运动与拓扑变形 Fig. 1 Movement and topology deformation of Snake curve

图选项

在曲线每次变形前先进行节点拆分。如图 1(b)-(d)所示，对每一个节点，除了其来源方向之外，在剩余的3个方向上各生成一个新节点，新节点在影像上的位置坐标均与原节点相同，不同的是各自的移动方向；若新节点的方向指向曲线内部则弃之，并将剩余节点按顺序插入原曲线节点序列中去，此时这部分曲线的拓扑结构就发生了改变。

能量函数的内部能量由弹力和张力组成，弹力可避免曲线出现曲率过大的部分，张力可使曲线相邻节点间的距离趋于相同。由于正交T-Snake模型的节点约束在网格顶点，节点间距固定为k或k(k为网格边长)，因此本文忽略张力。为保证曲线向外膨胀变形，添加一个膨胀力。外部能量主要指图像力，可保证Snake曲线停留在梯度极大值处(目标边缘)。因此，节点V_i的局部能量由弹力、膨胀力和图像力组成，函数表达式为

点击浏览公式

点击浏览公式

以这样一个由弹力、膨胀力和图像力组成的能量函数为原则，引导Snake曲线检测目标边界的过程是：①当节点在河流内部并远离边缘时，|▽I(V_i)|较小而e^-ε(σ+▽)值较大，此时膨胀力占主导作用。随着曲线向外扩张，节点与曲线几何中心的距离|V_i-C|逐渐变大导致膨胀力逐渐变小，因而能量函数值逐渐变小，使曲线继续变形以追求能量函数最小值；②当节点在河流边缘时，e^-ε(σ+▽)的值迅速衰减至零而|▽I(V_i)|突然变大，此时图像力占主导作用，能量函数值将不再继续减小，促使曲线停止运动，从而停留在该处。

当要提取的目标边界由两个轮廓曲线组成，且其中一个轮廓曲线在另一个的内部时，将内部轮廓所构成的区域称为岛。在河流影像上常常出现这种情况，例如河中岛或沙洲等。这种情况下Snake曲线一定会在某个地方发生自相交，如图 2(a)所示，节点V_i与V_j将在曲线下一次变形中发生重合，从而发生自相交拓扑冲突。解决办法是将V_i的前一节点与V_j的后一节点相连，V_i的后一节点与V_j的前一节点相连，然后将V_i、V_j抛弃，如图 2(b)所示。这样原曲线分裂为两条新曲线，内部曲线将检测内部岛屿边界，外部曲线将检测目标边界。

图 2 岛状空洞拓扑冲突的处理 Fig. 2 Handling for topological conflict caused by island shaped hollow

图选项

假设当前Snake曲线是S={V₁,V₂,…,V_n}，节点个数是n。对于当前要移动的节点V_i，如果它的目标位置是曲线上另一节点V_j，而j-i=m(即V_i、V_j之间的节点个数为m-1)，那么进行以下操作：从节点V_i、V_j处将曲线S分裂成两条新闭合曲线S₁、S₂,其中，S₁是内部曲线，由V_i和V_j之间的节点序列构成，节点个数是m-1；S₂是外部曲线，由V₁到V_i-1、和V_j+1到V_n两个节点序列合并构成，节点个数是n-m-1。曲线S₁构造为

点击浏览公式

曲线S₂构造为

点击浏览公式

外部曲线继续向外扩张变形，最终停留位置就是河流自身边界，内部曲线向内运动变形，最终停留位置就是河中岛的边界。

正交T-Snake模型只要求初始轮廓是目标内部一条闭合曲线，无形状和大小限制，一般手工选取目标内任一点，将距该点最近的网格顶点依次相连作为初始轮廓。文献^[21]分析了遥感影像中自然地物纹理的分形特征，用分形维数(fractal dimension，FD)度量，认为水体的分形维数是最低的。据此，本文提出水体轮廓自动初始化方法，如图 3所示。

图 3 轮廓初始化流程 Fig. 3 Flowchart of contour initialization

图选项

首先，将影像划分成一系列的影像块，计算每块影像的分形维数；其次，选择分形维数最小的影像块作为水体内部区域，若多个区域同时拥有最小值，则随机选择其中一个；然后，对该区域的每一个边界点，搜索与它距离最近的T-Snake模型网格顶点，根据网格顶点与该点之间的位置关系，为网格顶点赋予方向属性：若网格顶点在相应边界点的上方，则网格顶点的方向是“向上”，同理，可判断方向为“向下”“向左”或“向右”的网格顶点；最后，将这些网格顶点按顺时针顺序依次相连，所构成的闭合曲线作为水体初始轮廓。

每个影像块的分形维数采用分形布朗运动模型^[22]计算，该模型将图像的灰度值当作是随机游走的结果：对于任一点x，其邻域(x+△x)内各像素灰度对应的空间曲面可以用一个分形布朗曲面f(x)来拟合。f(x)具有一个基本性质，即

点击浏览公式

根据H可以计算分形维数FD=3-H。求H值的方法是对式(5)取对数，有

点击浏览公式

水体边界提取是初始轮廓随着能量函数最小化而逐渐优化的过程，算法如下：

(1) 模型初始化。初始化网格大小r；用3.2节方法初始化水体轮廓，记录各节点为(x,y,O)，其中，(x,y)是节点坐标，O={‘up’,‘down’,‘left’,‘right’}是节点方向；计算初始轮廓的能量E_before。

(2) 节点拆分，并记录新的节点序列中各点的坐标和方向。

(3) 遍历节点序列，对节点V_i(x_i,y_i,O_i)进行以下操作：①判断V_i的目标位置V′_i，若O_i=‘up’，则 V′_i(x_i,y_i-r,‘up’);若O_i=‘down’，则V′_i(x_i,y_i+r,‘down’);若O_i=‘left’,则V′_i(x_i-r,y_i,‘left’)；若O_i=‘right’,则V′_i(x_i+r,y_i,‘right’)。②检测岛状空洞拓扑冲突，若存在一个曲线节点V_j，坐标上满足V_j=V′_i，则说明发生拓扑冲突，执行③；若不存在，则执行④。③按照3.1节方法将原曲线分裂成为两条新曲线S₁和S₂，然后对S₁和S₂分别执行(2)。④计算V_i、V′_i的局部能量E_i、E′_i，并比较大小：若E′_i＜E_i，则令V_i=V′_i，E_i=E′_i。

(4)计算新曲线的能量。比较曲线变形前后的能量大小：若E_after＜E_before，则令E_before=E_after，继续执行(2)；否则，输出Snake曲线的节点坐标序列。最终输出的一条或多条Snake曲线就是提取到的水体边界。

试验数据为国产GF-1卫星PMS相机近红外波段的一幅湖泊影像(1100像素×1100像素)，如图 4(a)所示，和一幅包含河中岛的长江某段影像(2244像素×1339像素)，如图 6(a)所示，空间分辨率为8m。采取改进的正交T-Snake模型进行水体边界提取：正交网格大小r=2；弹力的权重系数η=0.1；膨胀力的权重系数γ=-0.2，膨胀力的衰减速度ε=2，在膨胀力中计算节点邻域灰度标准差和极差的邻域大小k=2；图像力的权重系数λ=0.4。

图 4 3种模型提取的湖泊边界 Fig. 4 The lake boundary extracted by three models

图选项

图 6 经典Snake模型和GVF Snake模型的河流边界提取结果 Fig. 6 The results of the river boundary extraction based on classic Snake model and GVF Snake model

图选项

采用经典Snake模型^[7]和GVF Snake模型^[12]的提取结果进行对比，考虑到它们对初始轮廓的敏感性，试验中保证初始轮廓尽量贴近水体边缘，且几何形状尽量与水体边界形状一致。试验在Intel酷睿i5、CPU主频2.27GHz、内存4GB的计算机上，eclipse+pydev编译环境下用Python编程实现。

第1幅影像的湖泊边缘有多处深凹边界，如图 4(a)所示。图 4(b)、(c)分别为传统Snake和GVF Snake的提取结果，黄色曲线是初始轮廓，红色曲线是检测的湖泊边界。图 4(d)为本文方法的提取结果，黄色区域是自动获取的分形维数最小区域，红色曲线是检测的湖泊边界。

图 5 不同的曲线变形方式 Fig. 5 Different curves deformation mode

图选项

从图 4可以看出，传统Snake和GVF Snake的边界提取效果均不如改进正交T-Snake模型，表现为对凹状边界的收敛效果不佳。因为它们是通过搜索邻域能量最小点作为节点新位置，如图 5(a)、(b)所示，阴影区域为点P的5×5邻域的局部最优点，在狭长凹区域会导致节点易在局部能量极小值的作用下停留在凹区域内某一侧方边界上，这将使得Snake曲线不能收敛到凹区域真正边界。GVF Snake由于采取GVF场作为外部力，扩大了节点搜索范围，故对浅凹边界的检测效果优于传统Snake，但对于狭长深凹边界依然会出现跟传统Snake一样的缺点。改进的正交T-Snake模型由于通过节点拆分和搜索最优的网格顶点来实现曲线变形，故不会陷入局部极小值，如图 5(c)、(d)所示，所以能够很好地收敛到各深凹区域的边缘。

第2幅影像中河流内有两个河中岛，对于传统Snake和GVF Snake，如图 6所示，考虑到它们对于初始轮廓的敏感性，传统Snake的初始轮廓完全包裹影像中整个河流，GVF Snake模型的初始轮廓尽量与河流形状相似，但即便如此，这两种模型均不能获取完整而精确的河流边界：传统Snake模型不能检测到河流内部河中岛的边界，GVF Snake模型无法提取河流的所有支流边界。特别是GVF Snake模型，受其曲线变形方式的局限性，容易出现局部极小值，因检测不到岛区域背面边界，曲线在继续变形过程中会不断演变成U型，因此无法提取完整的岛区域边界。

改进正交T-Snake模型提取河流边界的过程如图 7所示，图 7(a)是用23×23的窗口计算分形维数，并均衡化拉伸到0-255区间后得到的分形维数分布图。可以看到，纹理光滑平坦的河流内部分形维数偏低，表现为连续均匀的暗区域，其中黄色区域为计算得到的一个分形维数最小的区域，可以看到该区域位于河流的内部某一处。从图 7(b)-(e)可以看出，Snake曲线在能量函数的作用下不断地膨胀和向河流边缘运动，在此过程中，其形状也在不断地发生改变；在遇到河中岛时，由于设计了合理的拓扑冲突检测机制和处理办法，Snake曲线的两端在河中岛边界某处融合，并且整体分裂为两条曲线，其中内部曲线停留在该处从而获取了河中岛的边界，而外部曲线则继续运动检测到河流自身边界。这说明本文方法在提取含河中岛河流时是有效的。

图 7 基于改进正交T-Snake模型的河流边界提取结果 Fig. 7 The results of the river boundary extraction based on improved orthogonal T-Snake model

图选项

采用Pratt品质因数(Pratt’s figure of merit，PFM)^[23]和形状误差^[6]相结合的方式对3种模型的水体边界提取精度进行客观评价。

Pratt品质因数可评价边界的定位精度，计算公式为

点击浏览公式

点击浏览公式

在湖泊影像和河流影像上手动提取湖泊边界和河流边界，将它们作为各自的理想边界，然后利用式(7)或式(8)分别计算经典Snake模型、GVF Snake模型和本文方法的Pratt品质因数(PFM)和形状误差(SE)。另外统计了3种模型在两个案例中的时间开销。精度和效率的统计结果如表 1所示。

从表 1可以看出，两个案例中正交T-Snake模型的PFM最高、SE最低，说明正交T-Snake模型所提取的水体边界在定位精度和形状准确度上都是最佳的。GVF Snake的湖泊边界提取精度优于传统Snake，是因为它能够检测到部分浅凹边界，但它们均遗漏掉很大一部分河流边界，因此在河流案例中的精度都非常低。该结果亦验证了4.2.1 节的定性分析。从表 1中还可以看出，改进正交T-Snake模型的耗时最少，计算效率最高，这是因为：①本文方法实现了轮廓自动初始化因而在构造初始轮廓上耗时最少；②另两种模型需对节点邻域的所有像素都计算一次能量函数，若邻域大小M×M，节点数n，则一次迭代的时间复杂度是O(n×M²)，而本文方法的每个节点最多计算3个方向的能量函数，一次迭代的时间复杂度小于O(3n)；③曲线两次变形之间，另两种模型需进行的节点内插操作也增加了时间开销，而本文方法无需此项操作。

图 8统计了当网格宽度分别为1、2、3、4、5个像素时的两幅影像的水体提取效率和精度。总体上，随着网格大小的增加，水体提取的时间效率有所提高，而水体提取精度则会下降(Pratt品质因数降低，形状误差变大)。其中，当网格宽度大于2时，随着网格宽度变大，Pratt品质因数明显降低，形状误差大幅增加，而耗时的变化比较缓慢。因此，本试验的网格大小设置为2×2，此时基于改进正交T-Snake模型的水体提取精度较好，计算效率也较高。

图 8 改进正交T-Snake模型的网格大小对精度和耗时的影响 Fig. 8 The grid size of the T-snake model impact on accuracy and time cost

图选项

(1) 本文提出了一种可自动提取复杂水体边界的改进正交T-Snake模型，在国产GF-1卫星数据的水体边界提取中取得了不错的效果。试验证明，该模型所提取的水体边界的Pratt品质因数超过93%、形状误差在5%以下，定位精度和形状准确度很高，对于深凹边界的提取效果优于传统Snake和GVF Snake模型。

(2) 改进正交T-Snake模型特别针对河中岛现象设计了曲线变形冲突检测与处理算法，使遥感影像中包含沙洲或岛屿的水域边界的直接提取变得简单化，且可以直接得到连续的水体边界线。

(3) 改进正交T-Snake模型无需手动设置初始轮廓，自动化程度高；曲线变形过程简单使得模型计算效率较高。

(4) 在改进正交T-Snake模型中，划分影像的网格尺寸越大，边界提取的效率越高而精度越低，因为一些小于网格大小的影像特征易被忽略。实际应用中需灵活设置网格大小，在边界提取精度和模型计算效率两者之间取得平衡。

综上所述，本文改进的正交T-Snake模型具有自动化程度高、拓扑变形能力强、深凹边界收敛性能好、可检测含岛目标和计算效率高等特点，能够对包含河中岛的拓扑结构极其复杂的河流边界进行精确提取，亦可在其他地物边界提取中推广应用，如海上岛屿提取和道路网提取等。