随着航空摄影技术和空间荷载技术的飞速发展，遥感作为一门科学技术也展现出日新月异的面貌。随着遥感和测绘技术的发展，摄影与遥感测绘正在越来越多的场合互补甚至替代了传统的数字化测图，但是，摄影与遥感测绘仍然存在着局部精度较低(大比例尺地形图)，影像数据获取成本较高，时效性无法保证(部分小比例尺地图)等缺点，因此，为了弥补上述缺点，进一步提高测绘工作效率，如果能在基于空中三角测量加密的遥感影像的测绘系统中实现部分地物要素甚至全要素的计算机软件自动绘图，将具有重大意义^[1]。

地表地物都会发射电磁波，在同一时空条件下，特定的地物反射、发射、折射、吸收电磁波的特性与波长形成相对应的函数关系。而无论影像是多光谱、全色、远红外、近红外，当我们将相应的符合影像规律的这种函数关系用数学曲线形式表现出来时，就形成了直观的地物波谱，而不同的地物具有不同的波谱曲线和形态。特定地物波谱的反射率和它的波长之间的函数规律称为特定地物反射的波谱特性。而正因为不同的地物具有不同的波谱特性，使得基于计算机算法的自动绘图系统成为可能^[2]。

而想要实现上述的自动绘图，关键在于影像边缘检测算法的选择，特定影像局部区域的亮度和色彩等解译要素明显变化的部分就是影像的边缘，影像边缘的灰度剖面可以看作是一个阶跃，即其从一个灰度值在一个很小的缓冲区域内急剧变化到另一个灰度，而且两者间相差较大的灰度值。而边缘检测算法主要是通过计算机程序语言来实现影像灰度变化的度量、检测和定位，从而得出正确的影像边缘。自从1959年提出边缘检测概念以来，经过半个世纪的发展，现今已有许多种不同的边缘检测算法，如Roberts边缘检测算法、Sobel边缘检测算法、Prewitt边缘检测算法、Laplacian边缘检测算法、Marr (LoG)边缘检测算法以及本文选择的坎尼边缘检测算法(以下简称坎尼算法)等^[3-7]。

1 坎尼边缘检测算法分析

检测影像阶跃边缘的基本思路就是在影像中通过算法找出拥有局部最大梯度幅值的像素点。影像边缘检测过程有两个要求：首先，能有效地抑制干扰波谱；其次，必须尽量精确确定真实影像边缘的位置。而处理流程的困难之处在于选定算法在提高了对相关地物正确边缘的敏感性的同时，也提高了对干扰波谱边缘的敏感(如房顶上的反光色带、水田中的积水区域边缘，甚至遥感过程中的云雾边缘)。

1.1 坎尼边缘检测算法的基本原理

成功的影像边缘检测基于合理的算法和数学基础，本文的方法需要用到离散化梯度逼近函数，使其根据二维灰度矩阵的梯度向量来寻找待处理影像灰度矩阵的灰度跃变的位置，以跃变位置像素为基础，将影像上的角点、纹理、边缘等基元图的点连接起来，就可以构成最初的影像边缘。

然而实际作业过程中完全符合理论设想的灰度阶跃影像及线条状影像边缘是非常罕见的，而且大多数的遥感传感器又具有低频滤波的特性，这样会使得阶跃边缘变为斜坡性边缘，从而使灰度变化看起来跨越了一定的缓冲区域而不是瞬间的。因此，影像滤波工作在影像边缘检测中具有十分重要的意义。

本文之所以选择坎尼算法来采集影像边缘，是因为坎尼算法在边缘检测滤波方法的选择上相较于其他算法更适合遥感影像的边缘特性。其滤波方法能在去除干扰波谱与保存影像正确边缘之间找到较理想的平衡点，算法在原始影像的基础上通过边缘检测滤波器得到初步平滑的影像；然后通过全局梯度值和局部梯度值的计算进一步得到平滑的影像；再通过抑制影像局部梯度非极大值来得出阈值处理之前的过程数据；最后，通过大小区域阈值算法及循环算法得到最终的影像边缘。本文的边缘检测流程图如图 1所示。

1.2 坎尼算法处理顺序

1) 使用高斯滤波器平滑影像，其数学描述如下。二维高斯函数为：

$G\left( x,y \right)=\frac{1}{2\text{ }\!\!\pi\!\!\text{ }{{\delta }^{\text{2}}}}\exp \left( -\frac{{{x}^{2}}+{{y}^{2}}}{2{{\delta }^{2}}} \right)$

(1)

而高斯函数G_n(x, y)在某个方向向量上的一阶方向导数的数学表达式如下(设方向矢量为n，梯度矢量为∇G)：

${{G}_{n}}=\frac{\partial \mathit{\boldsymbol{G}}}{\partial \mathit{\boldsymbol{n}}}=\mathit{\boldsymbol{n}}\nabla \mathit{\boldsymbol{G}}$

(2)

将影像f(x, y)与G_n作卷积(将两者相乘后求和)的数学运算，通过计算G_n×f(x, y)的最大值得到与检测边缘正交垂直的方向，也就是本步骤需要的n的方向。

再利用一阶偏导的有限差分来计算梯度的幅值和方向，数学描述为：

$ \begin{align} & {{E}_{x}}=\frac{\partial \mathit{\boldsymbol{G}}}{\partial x}\times f\left( x,y \right),{{E}_{y}}=\frac{\partial \mathit{\boldsymbol{G}}}{\partial y}\times f\left( x,y \right) \\ & A\left( x,y \right)=\sqrt{E_{x}^{2}+E_{y}^{2}},\mathit{\boldsymbol{ }}\!\!\theta\!\!\rm{ }=\arctan \left( \frac{{{E}_{x}}}{{{E}_{y}}} \right) \\ \end{align} $

式中，A(x, y)反映了影像(x, y)点处的边缘强度；θ是影像(x, y)点处的法向矢量。

2) 确定影像边缘。需要同时考虑全局梯度与局部梯度之间的关系, 影像边缘必须包含局部梯度的最大值点，同时应该有效地舍弃(最理想状态是完全不包含)局部梯度非极大值(non-maxima suppression, NMS)的点，因此，需要引进梯度方向的概念。局部梯度非极大值抑制的原理如图 2所示。

在图 2中，人为地将一个圆用东-西、西北-东南、北-南、西南-东北4个方向的直径8等分，再取其中一个方向为界，则两侧各有4个45°的扇形区域，分别把4组对角扇形用0、1、2、3数字标记为4组，根据其空间位置和排列组合原则，任选3个相邻的扇形区域作为对象，则有4种有效的组合，以图 2中间的子图为例，即有0-3-2、3-2-1、2-1-0、1-0-3共4种可能(循环重复)。而每个组合都有一个中心像素(即中间的扇形)，设其值为M，如果M值不大于其两侧相邻扇形区域的标记值，则令M=0(依上例，则1-0-3的组合满足要求)。

3) 对目标影像的边缘进行检测处理和连接处理。本文选择大小阈值区域限制的算法，这个算法也是一个有效控制目标虚假边缘的方法，即对G(x, y)事先人为设定一个阈值，然后将低于阈值的所有值都赋零值。而阈值选取方法分为以下两步。

首先，找到和确定影像的边缘位置，在此可以引入区域集合概念，由上述可知，边缘强度为A(x, y), 设高阈值为α，低阈值为β，三者相比较有以下3种可能：①A(x, y)>α>β，依照上述数学原理可判定本特征点可作为边缘点；②α>β>A(x, y)，依照上述数学原理可排除本特征点作为边缘点的可能；③α>A(x, y)>β，依照上述数学原理，需要另外分析其相邻两个扇形区域的M值，分别设为M₁与M₂，如果M₁>α>A(x, y)>β或M₂>α>A(x, y)>β，则判定特征点是边缘点，如果α>M₁且α>M₂，则判定特征点不是边缘点。

其次，对目标影像边缘进行连接融合处理，依据阈值区域限制算法的程序，设高阈值为δ₁，低阈值为δ₂，且人为规定δ₁/2≈δ₂，根据上述数学原理计算出其二维高斯函数分别为G₁(x, y)、G₂(x, y)，因为G₁(x, y)是高阈值的函数，依上述原理，影像虚假边缘几率很小，但是有可能产生影像边缘间断的情况。阈值区域限制算法首先要使用高阈值函数G₁(x, y)来确定目标影像边缘的最外围轮廓(有可能不完全连接)在断裂的端点位置；再利用低阈值的高斯函数G₂(x, y)进行循环计算，依据图 2所示原理，通过G₂(x, y)的邻域位置和聚类分析技术把断裂的影像边缘进行连接，从而得到需要的影像边缘。

综上所述，坎尼算法具有如下优点：①极低的误码率，基本不会把边缘点误认为非边缘点；②相比同类算法来说拥有极高的定位精度，能够精确地把边缘点定位到灰度变化最大的像素上；③坎尼算法可以较高效地抑制影像的虚假边缘。

2 坎尼算法的代码设计及实例分析 2.1 代码设计

依上所述，依靠VC++的完整的坎尼算法代码需要以下步骤：①影像读取和灰度化；②影像的高斯滤波：根据一维高斯核进行两次滤波，生成一维高斯滤波系数，分别进行x向和y向的一维加权滤波，保存数据；根据二维高斯核进行滤波，采用高斯核进行高斯滤波，保存数据；③计算影像梯度及其方向；④非极大值抑制：定义相关的参数，对边界进行初始化，进行局部最大值寻找，再判定最优解；⑤双阈值检测实现：定义相应参数，构造灰度图的统计直方图，获取最大梯度幅值及潜在的边缘点个数，计算两个阈值，进行边缘检测；⑥生成矢量文件。

以下是非极大值抑制过程中相关参数的定义部分代码(其余代码略)：

01.unsigned char* N=new unsigned char[n Width*n Height]; //计算非极大值抑制结果

02.int g1=0, g2=0, g3=0, g4=0;//用于通过插值来得到亚像素点坐标值

03.double d Tmp1=0.0, d Tmp2=0.0;//存储亚像素点插值的灰度数据

04.double d Weight=0.2;//设置插值的权重和精度，其与影像质量成正比

2.2 实例分析

从数字影像中提取影像边缘是一个确定几何实体位置和方位的测绘过程，利用上述代码，以房屋为例作为特定地物，从数字影像中提取房屋需要以下过程：首先，根据影像定向参数确定出居民区所在整体影像的可能位置，以像素为单位，根据成图所要求的比例尺不同，依据精度要求选择不同的像素范围作为缓冲区；然后，基于坎尼算法与目标地物的波谱特征自动提取范围内所有地物的影像边缘，并生成矢量数据文件；接着，对以上影像边缘进行聚类分析，找出所有房屋边缘的近似位置；再利用稳健参数估计的方法对房屋边缘初始位置进行分析，尽量剔除非房屋边缘；最后，利用模板匹配算法对房屋边缘进行精确定位，得到房屋轮廓。

本文根据以上代码及过程提取了一小块影像较好的房屋区域(例图来自2014年获取的湖南省沅江市真彩色影像)，原图如图 3所示，图 4至图 6为具体的实验结果，由此比较可得，在影像较清楚、参数设置合适、算法进行过优化的前提下，房屋自动绘制的精度及效果都是良好的。

3 误差分析 3.1 算法本身的误差分析

1) 从数字影像中提取要素边缘的不确定性及误差分析。对于数字影像来说，可以使用上文列举的各种算法来提取影像边缘(算法间相比较各有优劣)，其通用流程基本上是首先提取影像边缘像素，参照给定的格网来描述局部边缘的位置、方位及强度的相似性，然后，按照测区要求的连续性将以上所有像素综合起来成为边缘。由此可以知道，所提取边缘的不确定性随梯度的增大而减小，也就是说，算法所提取边缘的精度随梯度的增大而更高。

2) 直线边缘部分的误差分析。对所提取的边缘像素进行加权最小二乘拟合后，可以获得子像素级的直线边缘。通过上述数学解释可知，每一边缘像素的权的精度与其灰度梯度的平方成正比。

3) 聚类技术形成的误差分析。采取聚类分析技术可以将所提取的非目标对象边缘删去，留下的少量影像边缘经过对特定目标模板进行平移、旋转和缩放变换，就可以进一步剔除错误的非目标对象边缘，并将影像边缘规划到我们需要的位置上。这一过程的实质是基于稳健估计的方法来完成的，通过估计来确定变换参数。其基本思路是将最小二乘平差后具有较大残差或较大统计量的观测值赋予很小的权，从而减少其对平差结果的影响。

在特征匹配中，影像边缘的两个端点是分别采集的，拟合得到的最终结果必须进行一致性检验。如果提取的影像边缘被接受为我们需要的影像边缘，那么一致性检验中未发现的粗差以及隐含在特定目标模板中的系统误差对最终结果是有影响的，但其影响不会超过算法中事先由人为确定的阈值。

3.2 其他因素的误差分析

遥感影像的绝对定向精度完全依赖于空中三角测量的精度，本文涉及的最终矢量文件同样依赖于空中三角测量的精度。如文中所述，实际情况中理想的灰度阶跃及其线条边缘影像极其罕见，因此，越接近这个目标，结果矢量文件的精度越高。在影像干扰因素较多或遮挡较多时，人工干预成为必要，在算法自动绘图过程中，算法相关参数的设置和干扰因素的判读都最终影响着成果文件的精度。

4 结束语

当前遥感影像自动绘图还处在探索阶段，需要有人工的干预和检查，此外，需要有一种可靠而快捷的影像边缘处理算法来作为计算机程序自动绘图的基础。基于坎尼算法，计算机系统可以在遥感影像清晰、遮挡较少时进行质量较高的自动绘图工作，本文提出的自动绘图方法的精度从整体来说，与空中三角测量加密成果的精度、影像质量成正比。从本文所提方法本身来说，与算法参数设置的合理性及人机交互的正确性，甚至操作人员对影像的补充判读分析都有一定的关系。