数学建模网络挑战赛

数学建模网络挑战赛

承 诺 书

我们仔细阅读了第九届“认证杯”数学中国数学建模网络挑战赛的竞赛规则。 我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们接受相应处理结果。

我们允许数学中国网站(www.madio.net)公布论文，以供网友之间学习交流，数学中国网站以非商业目的的论文交流不需要提前取得我们的同意。

我们的参赛队号为：1559 参赛队员 (签名) ：

队员2：詹菊英

队员3：林增

参赛队教练员 (签名)： 陈亮

参赛队伍组别（例如本科组）：本科组

队员1：苏阿芬

数学建模网络挑战赛 编 号 专 用 页

参赛队伍的参赛队号：（请各个参赛队提前填写好）： 1559

竞赛统一编号（由竞赛组委会送至评委团前编号）：

竞赛评阅编号（由竞赛评委团评阅前进行编号）：

2016年第九届“认证杯”数学中国 数学建模网络挑战赛第一阶段论文

题 目 基于Canny算子和光流法的图像运动方向识别

关 键 词 灰度变换 Canny算子 Gaussian平滑 光流法

摘 要：

数码摄像技术被广泛应用于多种场合中。有时由于客观条件的限制，拍摄设备只能在较低的分辨率下成像。本文通过建立数学模型，运用多种数学方法和软件对低分辨率下的图像序列进行处理分析，来判断视野的运动方向。

首先采用RGB颜色模型和加权法Y0.299R0.587G0.114B，用Matlab软件将彩色照片经过灰度变换成计算机能处理的灰度图像。

然后基于Canny算子的边缘提取，采用Gaussian平滑函数进行图像去噪，将Gaussian平滑函数与原始图像进行卷积运算，即可得到抑制噪声后的图像，对去噪后的图像采用22大小的卷积模板作为x方向和y方向偏微分的一阶近似进行卷积运算，得到各像素的近似偏导数，然后利用一阶偏导的有限差分来计算像素的灰度梯度的大小G(x,y)和方向；对得到的每个像素灰度梯度的图像通过双阈值1和2来完成算法的提取和边缘连接，即可完成图像的边缘提取。

最后建立直角坐标系，应用光流法来判断运动方向，光流包含了目标的重要信息，可以用来检测图像序列中的运动目标。将图像灰度的变化和二维速度场联系起来，用光流场反应像素点运动的方向和速度，并根据光流场的分布特征，提取出运动目标的区域。假设像素点的灰度值不变，对图像中的每一个像素点都建立一个光流方程f(x,y,t)，进行泰勒展开，因为每两帧图像之间相隔的时间极短，所以我们令t0，得到

fffuv0xyt。 为了解出每个像素点的水平速度矢量u和垂直速度矢量v，引入Horn-Schunck算

法，即引入全局平滑性约束；为了减少光照变化与噪声影响，建立当前像素的一个3×3的邻域，用邻域中像素的加权值和代替当前像素的两个方向矢量并通过迭代运算，

fx[fxnfynft]n1n

得到的稳定值就是二维运动矢量，即得到运动的水平分量u

2fx2fy2

和垂直分量v

n1



n

fx[fxnfynft]

fxfy

222

。

本文还通过获取一序列分辨率为3264的缓慢运动图像，应用本文建立的模型，判断出运动的方向为向右上角移动，验证了模型的可行性与有效性。

参赛队号: 1559

所选题目: B 题

Abstract

Digital imaging technology is widely used in many occasions. Sometimes, because of objective conditions, the imaging device can only shoot at a lower resolution. In this paper, a mathematical model, using a variety of mathematical methods and software for low-resolution image sequences under processing and analysis, to determine the direction of motion vision.

Firstly RGB color model and weightingY0.299R0.587G0.114B, Matlab software with a color picture into a grayscale image after gradation transformation computer can handle.

Then, based on Canny operator edge extraction, using Gaussian smoothing function denoising, the Gaussian smoothing function convolution with the original image, the image can be obtained after the noise suppression, denoised image using22 convolution template size andx

direction as they direction of the first-order approximation of partial differential convolution, each pixel is obtained approximate partial derivatives, then use the first-order partial derivatives of finite difference to calculate the size of the pixel gray gradientand direction; Each image pixel

G(x,y)

gray gradient obtained by the dual-threshold and to complete the extraction

1

and

2

edge

connection algorithms, to complete the image edge extraction. Finally, the establishment of a Cartesian coordinate system, the application of optical flow method to determine the direction of motion, optical flow contains important information about the target can be used to detect moving objects in image sequences. The gray-scale image and two-dimensional velocity field changes linked with light direction and speed of the flow field response pixel movement, and according to the distribution characteristics of the optical flow field region extracted moving object. Suppose pixel gray value remains unchanged,

The image of each pixel have established an optical flow equation

f(x,y,t)

taylor expansion, because between every two images separated by a very short time, so we maket0, get

. fff

x

u

y

v

t

0

In order to understand the vertical and horizontal speed vector velocity vector for each pixel,

the introduction of Horn-Schunck algorithm, namely the introduction of a global smoothness constraints; To reduce the influence of noise and illumination changes, to establish a

neighborhood of the current pixel of the 3 × 3, and replaced by a pixel weighting values in the neighborhood of the current pixel in two directions by the vector and iteration, The resulting stable value is two-dimensional motion vectors, to obtain the horizontal

un1n

and vertical components of movement

fx[fxnfynft]

2fx2fy2

.

v

n1



n

fx[fxnfynft]

2fx2fy2

The article also obtain a sequence of slow motion picture resolution, the application of the model established in this paper, it is determined that the direction of movement of the upper right corner to move to verify the feasibility and effectiveness of the model.

目录

一．问题重述..........................................................................................................................2 二．问题分析..........................................................................................................................2 三．符号说明..........................................................................................................................2 四．模型的假设......................................................................................................................4 五．模型的建立......................................................................................................................4

5.1.灰度变换.....................................................................................................................4 5.2.基于Canny算子的边缘提取.....................................................................................5 5.2.1.Canny算子的边缘检测基理...................................................................................5 5.2.2..Canny算法实现过程..............................................................................................7 5.3.运动方向.....................................................................................................................9 5.4..坐标选取....................................................................................................................9 5.5.光流法.........................................................................................................................10 5.5.1.光流法算法实现......................................................................................................11 5.5.2.模型结果及分析......................................................................................................13 六.模型的求解.........................................................................................................................13

6.1.灰度变换.....................................................................................................................14 6.2.边缘提取.....................................................................................................................15 6.3光流法.........................................................................................................................16 七.模型的优缺点.....................................................................................................................17 八..参考文献............................................................................................................................18 九.附件.....................................................................................................................................19

一.问题重述

在科技发展快速的时代中，很多科学技术应用广泛，数码摄像技术已被普遍使用。由于镜头，光源环境，图像处理技术等客观原因的影响，有些拍摄设备拍摄图像像素低，图像模糊，无法快速辨析物体。

根据低分辨率成像的单色灰度图像序列判断运动方向：运动方向的基本含义是物体的运动你好趋势。具体来说，我们将拍摄的彩色图像序列更改其分辨率为32×64。紧接着建立合理地数学模型和算法，进行图像处理。先将图像灰度处理，然后边缘提取，建立坐标系，最后通过短时间的位移矢量来判断其运动方向。

二.问题分析

2.1.总的问题分析

判断低分辨率的单色图像序列的运动方向涉及分辨率处理，灰度处理，边缘提取，

坐标选取及位移矢量计算等多个问题结合。如果可以将自然拍摄的图像序列经过图像处理得到分辨率为32×64灰度图像，并建立合理模型及算法解决边缘提取，坐标选取及位移矢量的计算，那么，就可以通过位移矢量判断物体运动方向。 2.2对具体的问题分析

1、将拍摄的系列图像实时传输到计算机中，应用Photoshop使图像序列的分辨率皆为32×64，即将整个矩形视野划分为32×64个相同大小的矩形格子，其图像中每个像素的取值为对应格子的亮度平均值

2、将通过photoshop处理过的图像应用matlab软件灰度变换，让拍摄的彩色图像变换成单色图像。图像经过灰度处理后的效果图再运用canny算法进行图像边缘提取，得到效果图。其二值矩阵作为坐标选取的依据，简言之，为坐标选取提供方式方法。

三.符号说明

四.模型的假设

1、假设灰度变换后的图像与原图像保持一致。

2、假设canny算法中恰取最优阈值。 3、假设灰度值恒定，图像运动缓慢。

五.模型的建立

数码摄像技术被广泛使用于多种场合中。有时由于客观条件的限制，拍摄设备只能在较低的分辨率下成像。现在，有一序列在低分辨率下的成像，分辨率为3264。通过图像的处理，来判断视野的运动方向。图像的预处理是图像分析的前期工作，因为图像输入的过程会使图像的质量下降，无论是在视觉效果还是识别方面都存在很多问题。 5.1.灰度变换

随着科学技术的发展和显著提高，目前的摄像头所捕捉到的图像大都是彩色的，但计算机处理需要灰度图像，所以需要将原始图像转换成灰度图像，一般采用RGB，HIS两种颜色模型[1]，本文采用RGB模型，将彩色图像用R/G/B三个通道来表示[2]。一幅彩色图像的每个像素均有R, G, B三个分量，三个分量的各个分量都有255种取值，这样一来每一个素点就有1600多万种颜色变化，图像的存诸空间变大对后续的处理和计算影响很大。彩色图像的灰度化后每个像素只有256级灰度值，在使用算法和处理的过程中能很好的减少存储空间，还能使运算量减少，提高图像处理的速度最终提高运动目标检测与跟踪的速度[3]。彩色图像可采用加权法变换为灰度图像，如下式1-1所示：

Y0.299R0.587G0.114B （1-1）

[4]

式中:Y是根据颜色分量R，G，B以及YUV中亮度信号Y之间的关系算出的亮度，R，G，B分别表示红色、绿色、蓝色分量，其中彩色分量信号对转换的结果影响了红、绿、蓝三种颜色的加权系数。

用matlab软件，可实现图像的灰度变换，编程后的代码如附件一所示。 5.2.基于Canny算子的边缘提取 边缘信息是重要的图像特征信息[5]。 5.2.1.Canny算子的边缘检测基理：

Canny边缘提取是一种利用局部极值提取边缘的算子，Canny把边缘检测问题转换为检测单位函数极大值的问题。在高斯噪声中，一个典型的边缘代表一个阶跃的强度变化，根据这个模型，一个好的边缘检测算子应具有三个指标：一是低失误概率；二是高定位精度，；三是对每个边缘有唯一的响应，得到的边缘为单像素宽[6]。基于Canny边缘提取，推出了下面的三个定义准则:

1）信噪比准则。真正的边缘不发生没有被提取出来的现象，不是边缘点不会被误提取到，即检测出的边缘信息的漏检率最小，误检率最小。 经过f(x)滤波后，边缘点处的图像信号的响应是：



HGG(x)f(x)dx （2-1）



而噪声的响应的平方根是：



Hnf2(x)dx （2-2）





12

是单位长度上噪声振幅的均方。 信噪比的数学表达式为：



SNR

HGHn



G(x)f(x)dx



（2-3）

f2(x)dx







其中，f(x)是边界为,的滤波器的脉冲响应；G(x)代表边缘函数；是高斯

噪声的均方差，信噪比越大，提取的边缘质量越高。

2）定位准则。提取的边缘点的位置距实际的边缘点的位置最接近。即最终获取的边缘

点位置与实际位置之间的距离差异越小，边缘定位精度越高，得到的结果越好。 设检测出的边缘位置在x0（实际的边缘在x=0），则有：

HG(x)Hn(x)在x0处取得最大值，所以

''

(x0)Hn(x0)0 （2-4） HG '

HG(x)在x=0取得最大值，所以HG(0)0

于是就有：

''''2'' HG（2-5） (x0)HG(0)HG(0)(x0)HG(0)x0 '''即 HG(0)x0Hn(x0)

从而



E(x2)0

E(H(x0)) （2-6） 22(H(0))

G'(x)f'(x)dx



''

n''G

2



2f'2(x)dx

这里的E(x)表示的是x的期望[6]。

因为x0越小定位越准确，所以定位精度的数学表达式为：



''G(x)f(x)dx

L



（2-7）

f'(x)dx

2







式中，G'(x)和f'(x)分别表示G(x)及f(x)的一阶导数，Localization值越大，表

明定位精度越高。

3）单边缘响应准则。任何真实具有的边缘点和提取的边缘点是一一对照的关系。即边缘线的宽度应该是一个像素，而不是多个像素。要保证单边缘只有一个像素响应，检测算子的脉冲响应导数的零交叉点平均距离D(f')应满足fn''(x)为f(x)的二阶导数。伪边界平均距离D(f')描述了随机噪声与检测函数卷积后伪边界出现的平均距离，伪边界平均距离越长，测量结果中出现伪边界的个数越少[6]。

'2f(x)dx D(f')f''2(x)dx（2-8）  12

以上面的指标和准则为基础，利用泛函数求导的方法导出一个由边缘定位精度和信噪比乘积组成的表达式，这个表达式近似于Gaussian函数的一阶导数，此即为该最佳函数的最好近似。

设n为任意方向，Gn为Gaussian函数在这个方向上的一阶导数，即

GnG/nnG （2-9）

 边缘强度(G*I)G*IG*Inn12

边缘的方向为： 22 （2-10） 12

(G*I)G*InG*In （2-11） n 12(G*I)nn12

符号*表示卷积，I为输入图像。边缘点定义为Gn作用于图像I后的局部极大值点，在这个点上有(/n)Gn*I0，把式（2-4）代入，得到

(2/2n)Gn*I0 （2-12）

5.2.2.Canny算法实现过程

Canny边缘检测算法实现过程主要包括以下几个过程，即首先采用高斯滤波进行图像去噪，然后利用一阶偏导的有限差分来计算梯度的幅值和方向；最终通过双阈值来完成算法初次提取和连接边缘，得到最终的边缘提取结果[7]。

1）Gaussian平滑

Canny方法首先通过高斯平滑函数来完成去噪过程。高斯平滑函数如下：

x2y2

exp2) （2-13） H(x,y) 2221

式中，：为平滑参数，较小，抑制噪声性能较弱；较大，噪声抑制效果较好，

高斯平滑时模板会变大，计算量增大。(x,y)为像点坐标。

将高斯平滑函数与原始图像进行卷积运算，即可得到抑制噪声后的图像。

2）一阶差分运算

在去除噪声后的图像中，采用以下水平和竖直模板对图像进行卷积运算，得到各像素的近似偏导数。采用22大小的卷积模板作为对x方向和y方向偏微分的一阶近似：

1 H1 （2-14） 1

 H2 （2-15）  假设像素的灰度值为f(x,y)，与模板进行卷积得到水平方向的梯度Gx(x,y)和竖直方向的梯度Gy(x,y)：

Gx(x,y)f(x,y)H1 （2-16） Gy(x,y)f(x,y)H2 （2-17） 由以上两式可得该像素的灰度梯度的大小G(x,y)和方向：

G(x,y)x(x,y)2Gy(x,y)2 （2-18）

H(x,y) (i,j)arct2 （2-19） H1(x,y)

3）边缘点

用Caussian函数卷积图像后，计算式(2-12),从计算结果中找到过零点，即边缘点。边缘点的强度由Gn*I(G*I)来估计。

4）双阈值边缘提取

对非极大值抑制幅值进行阈值化后的结果是一个图像的边缘阵列。阈值化后得到的边缘阵列仍然有假边缘存在，如果用单阈值来处理，选择合适的阈值是很困难的，常需要经过反复试验。有效的方法是选用两个阈值[10]。设非极大值抑制图像为N(i,j)，对得

2，到的每个像素灰度梯度的图像采用两个阈值1、且2*1作用于N(i,j)进行过滤。

得到两个阈值边缘图像，把梯度值小于1的像素的灰度设为0，得到边缘图像T1，同理，可得到图像T2。由于图像T2是使用较大阈值检测得到的，因此，它含有很少的假边缘，去除了大部分噪声，但是也损失了有用的边缘信息。而图像T1是通过较小阈值检测得到保留着较多边缘信息的影像，因此可以以图像T2为基础，以图像T1为补充来连接被损失

的图像边缘。如果图像信号的响应大于高阈值，那么它一定是边缘；如果低于低阈值，那么它一定不是边缘;如果再低阂值和高阂值之间，就看它的8个邻域像素有没有大于高阈值的边缘。[6]

首先在图像T1中扫描，当遇到一个非零灰度的像素P时，跟踪以P为开始点的轮廓线，直到该线的终点Q。接着在图像T1中比较与图像T2中Q点位置对应的Q'点的8一邻近区域。如果在Q'点的8一邻近区域有非零像素R'存在，则将其包括到图像T2中，作为点R，同理，重复在图像T2中继续寻找跟踪以R为开始点的轮廓线，这样循环下去直到在图像T1和图像T2中都无法继续为止。包含P的轮廓线的连接已经完成，可标记为已访问过。然后依次可以重复寻找图像中的每一条轮廓线，直到在图像T2中再也找不到新的轮廓线为止[12]，这样就得到最终的边缘提取结果。

Matlab软件，自带Canny算法，所以编写代码，直接调用Canny算法，即可提取边缘。

5.3.运动方向

运动方向检测主要是利用相邻的两帧或多帧的图像来计算运动物体在二维投影上坐标的变化情况来计算运动的一些参数。主要是计算目标物相对于参照物的位移矢量，即物体的运动方向和速度。

5.4.坐标选取

为了更加方便描述图像运动方向，现以序列图像中第一张图像作为方向基准，建立以像素32×64分为横轴与纵轴的直角坐标系，如图1所示。由灰度二值矩阵中图像边缘最底端作为图像的特征点。为了减少选取误差，取特征点的3×3的邻域。如图2所示，截取特征点邻域内的变化情况，确定并计算其运动方向。

图5-1 图像边缘彩色图

图5-2 二值图像数值图

5.5光流法

自然界中的任何物体，只要存在运动就必然存在运动方向的问题。运动可以用运动场来描述，反映真实世界的三维运动。运动场由图像中每个点的运动矢量构成，从序列图像中近似计算出不能直接得到运动场[8]，所以在这里引入光流场。

所谓光流是指图像中灰度模式的表面运动，它是物点的三维速度矢量在成像平面上的投影，它表示了物点在图像中位置的瞬间变化[8]。

光流反映的是空间运动物体在成像平面上像素的运动，光流法是一种基于帧间象素分析的一种方法，是利用图像像素在时间域上的变化相邻帧之间的运动的连续性找到相邻两帧间的某种联系，从而计算出相邻帧之间运动信息的一种方法。一般情况下，相机运动、场景中目标运动或两者的共同运动这几种情况都会产生光流。

光流包含了目标的重要信息。光流信息可以用来检测图像序列中的运动目标；恢复目标的三维结构信息以及目标与摄像机之间的相对运动[9]。光流法运动目标的检测主要是基于光流场的检测，指的是现实中物体像素点的三维速度矢量反映在视频中的二维瞬

时速度场，即图像上每一点的灰度变化趋势。光流场直观地实现了现实物质运动的信息，也是当前计算机视觉相关领域研究的一个重要内容[10]。

光流场指图像灰度上求解的表面运动，而运动场指三维物体的实际运动在图像平面上的投影。光流场是基于在如下的假设下推导出光流方程的：运动物体做刚体运动(物体在运动过程中没有形变或只有微小形变);物体的反射光变化是平滑的；光照变化是连续的。满足上述条件的情况下可以使用光流方程得到图像的光流场。理想状况下，光流场和运动场互相吻合，但实际上并非如此[11]。

光流法运动目标检测的基本原理是对图像中的每一个像素点都建立一个光流方程，通过计算出来的光流场来模拟运动场，图像上的每一个像素点都是对应现实三维空间物体上的点，用以对运动物体进行动态分析。当图像中不存在运动物体时，那么图像上的光流矢量则是连续变化的，而当图像中存在着运动目标时，运动目标的光流场与相邻背景的光流矢量值是存在差异的，不同光流值相交的一系列点即是运动目标的边缘，此方法还可以计算出运动物体在图像中的相对位置[10]。

光流法的核心就是求解出运动目标的光流，即速度[9]。将图像灰度的变化和二维速度场联系起来，用光流场反应像素点运动的方向和速度，并根据光流场的分布特征，提取出运动目标的区域。根据视觉感知原理，客观物体在空间上一般是相对连续运动，在运动过程中，投射到传感器平面上的图像实际上也是连续变化的，即灰度不变性假设。

5.5.1光流法算法实现

光流计算法的数学算法分析过程如下：假设f(x,y,t) 是像素点 (x,y)在 t 时刻的灰度值，当tt时刻时，该像素点的位置则是xx,yy,假设灰度值不变，那么光流方程可表示为：f(xx,yy,tt)

进行泰勒展开如下：

f(x,y,t)f(x,y,t)fffxyt(2)xyt （3-1）

化简公式，并高次项得:

fffxyt0yt x （3-2）

当t0时，有

fffuv0xyt （3-3）

其中，u是每个像素点的水平速度矢量，v是每个像素点的垂直速度矢量,它们的矢量和

则表示该像素点的光流方向和大小。

udydxvdt （3-5） dt，

为了解出定解u和v，附加其他约束条件，从不同的角度引入不同约束条件，从而产生不同光流分析方法。下面引入Horn-Schunck算法，即引入全局平滑性约束，假设光流在整个图像上光滑变化，即相邻的图像像素点之间具有相似的光流。图像中光流的逐像素变化可以通过速度矢量分量的空间梯度的模的平方和来定量表示[11]。简言之，其速度的变化率为零。

全局平滑约束条件的表达式方程为：

2u(u2u2 ()0 （3-6）xy

v2v2 )(0 （3-7）xy 2v(

式中u和v分别是速度矢量两个分量的空间梯度，此处假设空间和时间坐标是连续变量。

光滑约束方程的误差如下：

bfxufyvft （3-8） 偏离光滑性的误差为：

c2(u2u2v2v2)()(()xyxy （3-9） 设总误差为，则表达式为：’

2(2c2b2)dxdy （3-10） 其中，为图像的平滑程度，由实验经验得到。

将上述公式转换为拉格朗日公式，其偏微分方程如下：

22uf2xufxfyvfxft （3-11） 22vfxfyuf2yfyft （3-12） 为了减少光照变化与噪声影响，建立当前像素的一个3×3的邻域，用邻域中像素的加权值和代替当前像素的两个方向矢量[10]。

设uu、vv，则：

(2fx2fy2)(u)fx[fxfyft] （3-13） (2fxfy)(v)fy[fxfyft] （3-14） 求解u和v的值如下：

un1n2222fx[fxnfynft]

2fx2fy2

fx[fxnfynft] （3-15） vn1n2fx2fy2 （3-16）

式中的n是迭代次数。在点(x,y,t)上计算所有的偏微分。速度矢量的两个分量的初始值分别为u0(x,y,t),v0(x,y,t)通常取值为零。通过迭代运算，得到稳定值就是在点(x,y,t)处的二维运动矢量。

5.5.2模型结果及分析

图像经过灰度处理和边缘处理后，通过光流法分析与计算且运用迭代法，通过循迭代得到位移矢量的水平与竖直分量u和v。u和v的表达式如下：

un1nfx[fxnfynft]

2fx2fy2

vn1nfx[fxnfynft]

fxfy222

当u，v都为正，则图像的位移方向为右上；当u，v都为负，则图像的位移方向为左下； 当u为正，v为负，则图像的位移方向为右下；当u为负，v为正，则图像的位移方向为左上。

短时间缓慢运动图像的位移方向近似代替图像的运动方向。

六.模型的求解

数码摄像技术被广泛使用于多种场合中，有时由于客观条件的限制，设备只能在较低的分辨率下成像，现在，整个视野区域向某个方向缓慢移动，每隔一定时间拍摄一帧

图像，运动的画面体现为图像的序列，拍摄到的系列图像实时地传输到计算机中，通过上文所建立的数学模型，来判断运动的方向。

6.1.灰度变换

用matlab软件，将低分辨率下的系列图像用加权法进行灰度变换。如下图所示。

图6.1.1原图

图6.1.3原图

图6.1.5原图

图6.1.2灰度变换后的图像

图6.1.4灰度变换后的图像

图6.1.6灰度变换后的图像

图6.1.7原图

图6.1.8灰度变换后的图像

图6.1.9原图

图6.1.10灰度变换后的图像

6.2.边缘提取

用matlab软件，调用canny算子，编写程序，如附件二所示，提取灰度变换图像中的边缘。如下图所示：

图6.2.1灰度变换后的图像

图6.2.2边缘提取后的图像

图6.2.3灰度变换后的图像

图6.2.5灰度变换后的图像

图6.2.7灰度变换后的图像

图6.2.9灰度变换后的图像

6.3.光流法

图6.2.4边缘提取后的图像

图6.2.6边缘提取后的图像

图6.2.8边缘提取后的图像

图6.2.10边缘提取后的图像

经过MATLAB代码调试后，（程序如附录四所示）得到位移矢量的水平及竖直分量u和v的值及相应位移变化图。如下图所示：

图6.3.1位移变化图 图6.3.2位移变化图

图6.3.3位移变化图 图6.3.4位移变化图

由程序运行结果可知，图像运动方向为左下，则视野方向与图像运动方向相反，即为右上，与实际相符合。证明所建立的数学模型具有实时性和有效性。

七.模型的优缺点

7．1优点：

1．所采用的Canny算子具有较好的抗噪性能和较高的边缘定位精度，有利于图像的边缘提取。

2．所采用的光流法携带大量的运动信息，支持各类行摄像机，且当不知道场景先验知识时，可以检测出目标区域。

3．Canny算子提取边缘不丢失重要的边缘，没有虚假的边缘；实际边缘与检测到的边缘位置之间的偏差最小

缺点：

1．Canny准则是一个连续准则，也就是说在假设图像和滤波器都是一个连续函数的情形下给出的，但实际上数字图像是离散的。即无法确定在连续域所谓的最优的滤波器在离散的数字图像上还是不是最优的，具有不确定性。

2．应用光流法算法若没有相应的硬件支持时，它的计算耗时、检测实时性差。

3．Canny检子依赖物体的大小且对噪声敏感。

八.参考文献

[1]卞琼.基于worldview2的高分辨率遥感[J].卫星影像河流提取研究，2013.05：30

[2]胡进刚.一种面向对象的高分辨影像道路提取方法[J].遥感技术与应用，2006,21（3）：184-188

[3]贾艳丽.基于视频图像序列的运动目标检测与跟踪[D].硕士学位论文，哈尔滨工程大学2012.03：07

[4]魏玉强，运动目标检测与跟踪算法研究[D].燕山大学，2010:30-31

[5]林卉，赵长胜，舒宁，基于Canny算子的边缘检测及评价[J].黑龙江工程学院学报，第17卷第2期，2003.06：04

[6]连静.基于多尺度融合的图像边缘检测[D].硕士学位论文，吉林大学，2005.05，33-39

[7]徐辛超，付晨，徐爱功.一种改进的Canny边缘提取方法[J].遥感信息.2016.10

[8]Beauchemin S S ,Barron J L.The computation of optical flow,ACM Computing Surveys[J],1995,27(3):433-467

[9]裴巧娜.基于光流法的运动目标检测与跟踪技术[J].计算机应用技术.2009.05:07-08

[10]金国伟，基于混合高斯建模的目标检测和光流法的方向判断[D]，硕士学位论文，杭州电子科技大学，2014.12,21-23

[11]贾海涛，运动目标检测与识别算法的研究[D]，硕士学位论文，电子科技大学，2006.01,7-8,22-23

九.附件

附件一：图像灰度变换

f=imread('1-1.jpg');%读取图像,图像要保存在matlab所在的路径中,根据所要灰度处理的照片，读取对应的照片。

figure(1);

imshow(f); %显示图像

[n m a]=size(f);%判断图像的大小

for x=1:n %通过双循环对图像进行灰度化处理

for y=1:m

p(x,y)=0.3*f(x,y,1)+0.59*f(x,y,2)+0.11*f(x,y,3);

end

end

figure(2);

imshow(p);%显示处理后的图像

imwrite(p,'abc.jpg'); %保存处理后的图像

附件二：边缘提取

f=imread('1-2.jpg');% 根据所要提取的图片，读取对应的图像

imshow(f)

figure

bw=edge(f,'canny');

imshow(bw)

附件三：3*3矩阵

附件四：光流法程序

程序一：

I1=rgb2gray(imread('44.jpg'));

I2=rgb2gray(imread('55.jpg'));

imshow(I1)

figure

imshow(I2)

figure

im1=single(I1);

im2=single(I2);

numLevels=3;

window=9;

iterations=1;

alpha = 0.001;

hw = floor(window/2);

t0=clock;

pyramid1 = im1;

pyramid2 = im2;

for i=2:numLevels

im1 = impyramid(im1, 'reduce');

im2 = impyramid(im2, 'reduce');

pyramid1(1:size(im1,1), 1:size(im1,2), i) = im1;

pyramid2(1:size(im2,1), 1:size(im2,2), i) = im2;

end;

for p = 1:numLevels

im1 = pyramid1(1:(size(pyramid1,1)/(2^(numLevels - p))),

1:(size(pyramid1,2)/(2^(numLevels - p))), (numLevels - p)+1);

im2 = pyramid2(1:(size(pyramid2,1)/(2^(numLevels - p))),

1:(size(pyramid2,2)/(2^(numLevels - p))), (numLevels - p)+1);

if p==1

u=zeros(size(im1));

v=zeros(size(im1));

else

u = 2 * imresize(u,size(u)*2,'bilinear');

v = 2 * imresize(v,size(v)*2,'bilinear');

end

for r = 1:iterations

u=round(u);

v=round(v);

for i = 1+hw:size(im1,1)-hw

for j = 1+hw:size(im2,2)-hw

patch1 = im1(i-hw:i+hw, j-hw:j+hw);

lr = i-hw+v(i,j);

hr = i+hw+v(i,j);

lc = j-hw+u(i,j);

hc = j+hw+u(i,j);

if (lr size(im1,1))||(lc size(im1,2)) else

patch2 = im2(lr:hr, lc:hc);

fx = conv2(patch1, 0.25* [-1 1; -1 1]) + conv2(patch2, 0.25*[-1 1; -1 1]);

fy = conv2(patch1, 0.25* [-1 -1; 1 1]) + conv2(patch2, 0.25*[-1 -1; 1 1]);

ft = conv2(patch1, 0.25*ones(2)) + conv2(patch2, -0.25*ones(2));

Fx = fx(2:window-1,2:window-1)';

Fy = fy(2:window-1,2:window-1)';

Ft = ft(2:window-1,2:window-1)';

A = [Fx(:) Fy(:)];

G=A'*A;

G(1,1)=G(1,1)+alpha; G(2,2)=G(2,2)+alpha;

U=1/(G(1,1)*G(2,2)-G(1,2)*G(2,1))*[G(2,2) -G(1,2);-G(2,1)

G(1,1)]*A'*-Ft(:);

u(i,j)=u(i,j)+U(1); v(i,j)=v(i,j)+U(2);

end

end

end

end

etime(clock,t0)

end

u=u(window:size(u,1)-window+1,window:size(u,2)-window+1)

v=v(window:size(v,1)-window+1,window:size(v,2)-window+1)

RGB1=showmap3(u,v,5);

imshow(RGB1);

程序二：

function RGB = showmap3(u,v,vmax)

anglestep = pi/3;

colorstep=1/anglestep;

RGBpart=[ 0 1 255 255 -1 0;

-1 0 0 1 255 255;

255 255 -1 0 0 1];

RGB=zeros(size(u,1),size(u,2),3);

for i=1:size(u,1)

for j=1:size(u,2)

angle=atan2(u(i,j),v(i,j));

if (angle

angle = angle+2*pi;

end

part=angle/anglestep;

if ( part >= 6.0 )

part = 0.0;

end

length=sqrt( u(i,j)*u(i,j) + v(i,j)*v(i,j) );

if(length>vmax)

RGB(i,j,1)=1;

RGB(i,j,2)=1;

RGB(i,j,3)=1;

else

for k=1:3

switch RGBpart(k,fix(part+1))

case 0

RGB(i,j,k)=0;

case 1

RGB(i,j,k)=colorstep*mod(part,1)*length/vmax;

case -1

RGB(i,j,k)=(1-colorstep*mod(part,1))*length/vmax;

case 255

RGB(i,j,k)=length/vmax;

end

end

end

end

end

end

数学建模网络挑战赛

承 诺 书

我们仔细阅读了第九届“认证杯”数学中国数学建模网络挑战赛的竞赛规则。 我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们接受相应处理结果。

我们允许数学中国网站(www.madio.net)公布论文，以供网友之间学习交流，数学中国网站以非商业目的的论文交流不需要提前取得我们的同意。

我们的参赛队号为：1559 参赛队员 (签名) ：

队员2：詹菊英

队员3：林增

参赛队教练员 (签名)： 陈亮

参赛队伍组别（例如本科组）：本科组

队员1：苏阿芬

数学建模网络挑战赛 编 号 专 用 页

参赛队伍的参赛队号：（请各个参赛队提前填写好）： 1559

竞赛统一编号（由竞赛组委会送至评委团前编号）：

竞赛评阅编号（由竞赛评委团评阅前进行编号）：

2016年第九届“认证杯”数学中国 数学建模网络挑战赛第一阶段论文

题 目 基于Canny算子和光流法的图像运动方向识别

关 键 词 灰度变换 Canny算子 Gaussian平滑 光流法

摘 要：

数码摄像技术被广泛应用于多种场合中。有时由于客观条件的限制，拍摄设备只能在较低的分辨率下成像。本文通过建立数学模型，运用多种数学方法和软件对低分辨率下的图像序列进行处理分析，来判断视野的运动方向。

首先采用RGB颜色模型和加权法Y0.299R0.587G0.114B，用Matlab软件将彩色照片经过灰度变换成计算机能处理的灰度图像。

然后基于Canny算子的边缘提取，采用Gaussian平滑函数进行图像去噪，将Gaussian平滑函数与原始图像进行卷积运算，即可得到抑制噪声后的图像，对去噪后的图像采用22大小的卷积模板作为x方向和y方向偏微分的一阶近似进行卷积运算，得到各像素的近似偏导数，然后利用一阶偏导的有限差分来计算像素的灰度梯度的大小G(x,y)和方向；对得到的每个像素灰度梯度的图像通过双阈值1和2来完成算法的提取和边缘连接，即可完成图像的边缘提取。

最后建立直角坐标系，应用光流法来判断运动方向，光流包含了目标的重要信息，可以用来检测图像序列中的运动目标。将图像灰度的变化和二维速度场联系起来，用光流场反应像素点运动的方向和速度，并根据光流场的分布特征，提取出运动目标的区域。假设像素点的灰度值不变，对图像中的每一个像素点都建立一个光流方程f(x,y,t)，进行泰勒展开，因为每两帧图像之间相隔的时间极短，所以我们令t0，得到

fffuv0xyt。 为了解出每个像素点的水平速度矢量u和垂直速度矢量v，引入Horn-Schunck算

法，即引入全局平滑性约束；为了减少光照变化与噪声影响，建立当前像素的一个3×3的邻域，用邻域中像素的加权值和代替当前像素的两个方向矢量并通过迭代运算，

fx[fxnfynft]n1n

得到的稳定值就是二维运动矢量，即得到运动的水平分量u

2fx2fy2

和垂直分量v

n1



n

fx[fxnfynft]

fxfy

222

。

本文还通过获取一序列分辨率为3264的缓慢运动图像，应用本文建立的模型，判断出运动的方向为向右上角移动，验证了模型的可行性与有效性。

参赛队号: 1559

所选题目: B 题

Abstract

Digital imaging technology is widely used in many occasions. Sometimes, because of objective conditions, the imaging device can only shoot at a lower resolution. In this paper, a mathematical model, using a variety of mathematical methods and software for low-resolution image sequences under processing and analysis, to determine the direction of motion vision.

Firstly RGB color model and weightingY0.299R0.587G0.114B, Matlab software with a color picture into a grayscale image after gradation transformation computer can handle.

Then, based on Canny operator edge extraction, using Gaussian smoothing function denoising, the Gaussian smoothing function convolution with the original image, the image can be obtained after the noise suppression, denoised image using22 convolution template size andx

direction as they direction of the first-order approximation of partial differential convolution, each pixel is obtained approximate partial derivatives, then use the first-order partial derivatives of finite difference to calculate the size of the pixel gray gradientand direction; Each image pixel

G(x,y)

gray gradient obtained by the dual-threshold and to complete the extraction

1

and

2

edge

connection algorithms, to complete the image edge extraction. Finally, the establishment of a Cartesian coordinate system, the application of optical flow method to determine the direction of motion, optical flow contains important information about the target can be used to detect moving objects in image sequences. The gray-scale image and two-dimensional velocity field changes linked with light direction and speed of the flow field response pixel movement, and according to the distribution characteristics of the optical flow field region extracted moving object. Suppose pixel gray value remains unchanged,

The image of each pixel have established an optical flow equation

f(x,y,t)

taylor expansion, because between every two images separated by a very short time, so we maket0, get

. fff

x

u

y

v

t

0

In order to understand the vertical and horizontal speed vector velocity vector for each pixel,

the introduction of Horn-Schunck algorithm, namely the introduction of a global smoothness constraints; To reduce the influence of noise and illumination changes, to establish a

neighborhood of the current pixel of the 3 × 3, and replaced by a pixel weighting values in the neighborhood of the current pixel in two directions by the vector and iteration, The resulting stable value is two-dimensional motion vectors, to obtain the horizontal

un1n

and vertical components of movement

fx[fxnfynft]

2fx2fy2

.

v

n1



n

fx[fxnfynft]

2fx2fy2

The article also obtain a sequence of slow motion picture resolution, the application of the model established in this paper, it is determined that the direction of movement of the upper right corner to move to verify the feasibility and effectiveness of the model.

目录

一．问题重述..........................................................................................................................2 二．问题分析..........................................................................................................................2 三．符号说明..........................................................................................................................2 四．模型的假设......................................................................................................................4 五．模型的建立......................................................................................................................4

5.1.灰度变换.....................................................................................................................4 5.2.基于Canny算子的边缘提取.....................................................................................5 5.2.1.Canny算子的边缘检测基理...................................................................................5 5.2.2..Canny算法实现过程..............................................................................................7 5.3.运动方向.....................................................................................................................9 5.4..坐标选取....................................................................................................................9 5.5.光流法.........................................................................................................................10 5.5.1.光流法算法实现......................................................................................................11 5.5.2.模型结果及分析......................................................................................................13 六.模型的求解.........................................................................................................................13

6.1.灰度变换.....................................................................................................................14 6.2.边缘提取.....................................................................................................................15 6.3光流法.........................................................................................................................16 七.模型的优缺点.....................................................................................................................17 八..参考文献............................................................................................................................18 九.附件.....................................................................................................................................19

一.问题重述

在科技发展快速的时代中，很多科学技术应用广泛，数码摄像技术已被普遍使用。由于镜头，光源环境，图像处理技术等客观原因的影响，有些拍摄设备拍摄图像像素低，图像模糊，无法快速辨析物体。

根据低分辨率成像的单色灰度图像序列判断运动方向：运动方向的基本含义是物体的运动你好趋势。具体来说，我们将拍摄的彩色图像序列更改其分辨率为32×64。紧接着建立合理地数学模型和算法，进行图像处理。先将图像灰度处理，然后边缘提取，建立坐标系，最后通过短时间的位移矢量来判断其运动方向。

二.问题分析

2.1.总的问题分析

判断低分辨率的单色图像序列的运动方向涉及分辨率处理，灰度处理，边缘提取，

坐标选取及位移矢量计算等多个问题结合。如果可以将自然拍摄的图像序列经过图像处理得到分辨率为32×64灰度图像，并建立合理模型及算法解决边缘提取，坐标选取及位移矢量的计算，那么，就可以通过位移矢量判断物体运动方向。 2.2对具体的问题分析

1、将拍摄的系列图像实时传输到计算机中，应用Photoshop使图像序列的分辨率皆为32×64，即将整个矩形视野划分为32×64个相同大小的矩形格子，其图像中每个像素的取值为对应格子的亮度平均值

2、将通过photoshop处理过的图像应用matlab软件灰度变换，让拍摄的彩色图像变换成单色图像。图像经过灰度处理后的效果图再运用canny算法进行图像边缘提取，得到效果图。其二值矩阵作为坐标选取的依据，简言之，为坐标选取提供方式方法。

三.符号说明

四.模型的假设

1、假设灰度变换后的图像与原图像保持一致。

2、假设canny算法中恰取最优阈值。 3、假设灰度值恒定，图像运动缓慢。

五.模型的建立

数码摄像技术被广泛使用于多种场合中。有时由于客观条件的限制，拍摄设备只能在较低的分辨率下成像。现在，有一序列在低分辨率下的成像，分辨率为3264。通过图像的处理，来判断视野的运动方向。图像的预处理是图像分析的前期工作，因为图像输入的过程会使图像的质量下降，无论是在视觉效果还是识别方面都存在很多问题。 5.1.灰度变换

随着科学技术的发展和显著提高，目前的摄像头所捕捉到的图像大都是彩色的，但计算机处理需要灰度图像，所以需要将原始图像转换成灰度图像，一般采用RGB，HIS两种颜色模型[1]，本文采用RGB模型，将彩色图像用R/G/B三个通道来表示[2]。一幅彩色图像的每个像素均有R, G, B三个分量，三个分量的各个分量都有255种取值，这样一来每一个素点就有1600多万种颜色变化，图像的存诸空间变大对后续的处理和计算影响很大。彩色图像的灰度化后每个像素只有256级灰度值，在使用算法和处理的过程中能很好的减少存储空间，还能使运算量减少，提高图像处理的速度最终提高运动目标检测与跟踪的速度[3]。彩色图像可采用加权法变换为灰度图像，如下式1-1所示：

Y0.299R0.587G0.114B （1-1）

[4]

式中:Y是根据颜色分量R，G，B以及YUV中亮度信号Y之间的关系算出的亮度，R，G，B分别表示红色、绿色、蓝色分量，其中彩色分量信号对转换的结果影响了红、绿、蓝三种颜色的加权系数。

用matlab软件，可实现图像的灰度变换，编程后的代码如附件一所示。 5.2.基于Canny算子的边缘提取 边缘信息是重要的图像特征信息[5]。 5.2.1.Canny算子的边缘检测基理：

Canny边缘提取是一种利用局部极值提取边缘的算子，Canny把边缘检测问题转换为检测单位函数极大值的问题。在高斯噪声中，一个典型的边缘代表一个阶跃的强度变化，根据这个模型，一个好的边缘检测算子应具有三个指标：一是低失误概率；二是高定位精度，；三是对每个边缘有唯一的响应，得到的边缘为单像素宽[6]。基于Canny边缘提取，推出了下面的三个定义准则:

1）信噪比准则。真正的边缘不发生没有被提取出来的现象，不是边缘点不会被误提取到，即检测出的边缘信息的漏检率最小，误检率最小。 经过f(x)滤波后，边缘点处的图像信号的响应是：



HGG(x)f(x)dx （2-1）



而噪声的响应的平方根是：



Hnf2(x)dx （2-2）





12

是单位长度上噪声振幅的均方。 信噪比的数学表达式为：



SNR

HGHn



G(x)f(x)dx



（2-3）

f2(x)dx







其中，f(x)是边界为,的滤波器的脉冲响应；G(x)代表边缘函数；是高斯

噪声的均方差，信噪比越大，提取的边缘质量越高。

2）定位准则。提取的边缘点的位置距实际的边缘点的位置最接近。即最终获取的边缘

点位置与实际位置之间的距离差异越小，边缘定位精度越高，得到的结果越好。 设检测出的边缘位置在x0（实际的边缘在x=0），则有：

HG(x)Hn(x)在x0处取得最大值，所以

''

(x0)Hn(x0)0 （2-4） HG '

HG(x)在x=0取得最大值，所以HG(0)0

于是就有：

''''2'' HG（2-5） (x0)HG(0)HG(0)(x0)HG(0)x0 '''即 HG(0)x0Hn(x0)

从而



E(x2)0

E(H(x0)) （2-6） 22(H(0))

G'(x)f'(x)dx



''

n''G

2



2f'2(x)dx

这里的E(x)表示的是x的期望[6]。

因为x0越小定位越准确，所以定位精度的数学表达式为：



''G(x)f(x)dx

L



（2-7）

f'(x)dx

2







式中，G'(x)和f'(x)分别表示G(x)及f(x)的一阶导数，Localization值越大，表

明定位精度越高。

3）单边缘响应准则。任何真实具有的边缘点和提取的边缘点是一一对照的关系。即边缘线的宽度应该是一个像素，而不是多个像素。要保证单边缘只有一个像素响应，检测算子的脉冲响应导数的零交叉点平均距离D(f')应满足fn''(x)为f(x)的二阶导数。伪边界平均距离D(f')描述了随机噪声与检测函数卷积后伪边界出现的平均距离，伪边界平均距离越长，测量结果中出现伪边界的个数越少[6]。

'2f(x)dx D(f')f''2(x)dx（2-8）  12

以上面的指标和准则为基础，利用泛函数求导的方法导出一个由边缘定位精度和信噪比乘积组成的表达式，这个表达式近似于Gaussian函数的一阶导数，此即为该最佳函数的最好近似。

设n为任意方向，Gn为Gaussian函数在这个方向上的一阶导数，即

GnG/nnG （2-9）

 边缘强度(G*I)G*IG*Inn12

边缘的方向为： 22 （2-10） 12

(G*I)G*InG*In （2-11） n 12(G*I)nn12

符号*表示卷积，I为输入图像。边缘点定义为Gn作用于图像I后的局部极大值点，在这个点上有(/n)Gn*I0，把式（2-4）代入，得到

(2/2n)Gn*I0 （2-12）

5.2.2.Canny算法实现过程

Canny边缘检测算法实现过程主要包括以下几个过程，即首先采用高斯滤波进行图像去噪，然后利用一阶偏导的有限差分来计算梯度的幅值和方向；最终通过双阈值来完成算法初次提取和连接边缘，得到最终的边缘提取结果[7]。

1）Gaussian平滑

Canny方法首先通过高斯平滑函数来完成去噪过程。高斯平滑函数如下：

x2y2

exp2) （2-13） H(x,y) 2221

式中，：为平滑参数，较小，抑制噪声性能较弱；较大，噪声抑制效果较好，

高斯平滑时模板会变大，计算量增大。(x,y)为像点坐标。

将高斯平滑函数与原始图像进行卷积运算，即可得到抑制噪声后的图像。

2）一阶差分运算

在去除噪声后的图像中，采用以下水平和竖直模板对图像进行卷积运算，得到各像素的近似偏导数。采用22大小的卷积模板作为对x方向和y方向偏微分的一阶近似：

1 H1 （2-14） 1

 H2 （2-15）  假设像素的灰度值为f(x,y)，与模板进行卷积得到水平方向的梯度Gx(x,y)和竖直方向的梯度Gy(x,y)：

Gx(x,y)f(x,y)H1 （2-16） Gy(x,y)f(x,y)H2 （2-17） 由以上两式可得该像素的灰度梯度的大小G(x,y)和方向：

G(x,y)x(x,y)2Gy(x,y)2 （2-18）

H(x,y) (i,j)arct2 （2-19） H1(x,y)

3）边缘点

用Caussian函数卷积图像后，计算式(2-12),从计算结果中找到过零点，即边缘点。边缘点的强度由Gn*I(G*I)来估计。

4）双阈值边缘提取

对非极大值抑制幅值进行阈值化后的结果是一个图像的边缘阵列。阈值化后得到的边缘阵列仍然有假边缘存在，如果用单阈值来处理，选择合适的阈值是很困难的，常需要经过反复试验。有效的方法是选用两个阈值[10]。设非极大值抑制图像为N(i,j)，对得

2，到的每个像素灰度梯度的图像采用两个阈值1、且2*1作用于N(i,j)进行过滤。

得到两个阈值边缘图像，把梯度值小于1的像素的灰度设为0，得到边缘图像T1，同理，可得到图像T2。由于图像T2是使用较大阈值检测得到的，因此，它含有很少的假边缘，去除了大部分噪声，但是也损失了有用的边缘信息。而图像T1是通过较小阈值检测得到保留着较多边缘信息的影像，因此可以以图像T2为基础，以图像T1为补充来连接被损失

的图像边缘。如果图像信号的响应大于高阈值，那么它一定是边缘；如果低于低阈值，那么它一定不是边缘;如果再低阂值和高阂值之间，就看它的8个邻域像素有没有大于高阈值的边缘。[6]

首先在图像T1中扫描，当遇到一个非零灰度的像素P时，跟踪以P为开始点的轮廓线，直到该线的终点Q。接着在图像T1中比较与图像T2中Q点位置对应的Q'点的8一邻近区域。如果在Q'点的8一邻近区域有非零像素R'存在，则将其包括到图像T2中，作为点R，同理，重复在图像T2中继续寻找跟踪以R为开始点的轮廓线，这样循环下去直到在图像T1和图像T2中都无法继续为止。包含P的轮廓线的连接已经完成，可标记为已访问过。然后依次可以重复寻找图像中的每一条轮廓线，直到在图像T2中再也找不到新的轮廓线为止[12]，这样就得到最终的边缘提取结果。

Matlab软件，自带Canny算法，所以编写代码，直接调用Canny算法，即可提取边缘。

5.3.运动方向

运动方向检测主要是利用相邻的两帧或多帧的图像来计算运动物体在二维投影上坐标的变化情况来计算运动的一些参数。主要是计算目标物相对于参照物的位移矢量，即物体的运动方向和速度。

5.4.坐标选取

为了更加方便描述图像运动方向，现以序列图像中第一张图像作为方向基准，建立以像素32×64分为横轴与纵轴的直角坐标系，如图1所示。由灰度二值矩阵中图像边缘最底端作为图像的特征点。为了减少选取误差，取特征点的3×3的邻域。如图2所示，截取特征点邻域内的变化情况，确定并计算其运动方向。

图5-1 图像边缘彩色图

图5-2 二值图像数值图

5.5光流法

自然界中的任何物体，只要存在运动就必然存在运动方向的问题。运动可以用运动场来描述，反映真实世界的三维运动。运动场由图像中每个点的运动矢量构成，从序列图像中近似计算出不能直接得到运动场[8]，所以在这里引入光流场。

所谓光流是指图像中灰度模式的表面运动，它是物点的三维速度矢量在成像平面上的投影，它表示了物点在图像中位置的瞬间变化[8]。

光流反映的是空间运动物体在成像平面上像素的运动，光流法是一种基于帧间象素分析的一种方法，是利用图像像素在时间域上的变化相邻帧之间的运动的连续性找到相邻两帧间的某种联系，从而计算出相邻帧之间运动信息的一种方法。一般情况下，相机运动、场景中目标运动或两者的共同运动这几种情况都会产生光流。

光流包含了目标的重要信息。光流信息可以用来检测图像序列中的运动目标；恢复目标的三维结构信息以及目标与摄像机之间的相对运动[9]。光流法运动目标的检测主要是基于光流场的检测，指的是现实中物体像素点的三维速度矢量反映在视频中的二维瞬

时速度场，即图像上每一点的灰度变化趋势。光流场直观地实现了现实物质运动的信息，也是当前计算机视觉相关领域研究的一个重要内容[10]。

光流场指图像灰度上求解的表面运动，而运动场指三维物体的实际运动在图像平面上的投影。光流场是基于在如下的假设下推导出光流方程的：运动物体做刚体运动(物体在运动过程中没有形变或只有微小形变);物体的反射光变化是平滑的；光照变化是连续的。满足上述条件的情况下可以使用光流方程得到图像的光流场。理想状况下，光流场和运动场互相吻合，但实际上并非如此[11]。

光流法运动目标检测的基本原理是对图像中的每一个像素点都建立一个光流方程，通过计算出来的光流场来模拟运动场，图像上的每一个像素点都是对应现实三维空间物体上的点，用以对运动物体进行动态分析。当图像中不存在运动物体时，那么图像上的光流矢量则是连续变化的，而当图像中存在着运动目标时，运动目标的光流场与相邻背景的光流矢量值是存在差异的，不同光流值相交的一系列点即是运动目标的边缘，此方法还可以计算出运动物体在图像中的相对位置[10]。

光流法的核心就是求解出运动目标的光流，即速度[9]。将图像灰度的变化和二维速度场联系起来，用光流场反应像素点运动的方向和速度，并根据光流场的分布特征，提取出运动目标的区域。根据视觉感知原理，客观物体在空间上一般是相对连续运动，在运动过程中，投射到传感器平面上的图像实际上也是连续变化的，即灰度不变性假设。

5.5.1光流法算法实现

光流计算法的数学算法分析过程如下：假设f(x,y,t) 是像素点 (x,y)在 t 时刻的灰度值，当tt时刻时，该像素点的位置则是xx,yy,假设灰度值不变，那么光流方程可表示为：f(xx,yy,tt)

进行泰勒展开如下：

f(x,y,t)f(x,y,t)fffxyt(2)xyt （3-1）

化简公式，并高次项得:

fffxyt0yt x （3-2）

当t0时，有

fffuv0xyt （3-3）

其中，u是每个像素点的水平速度矢量，v是每个像素点的垂直速度矢量,它们的矢量和

则表示该像素点的光流方向和大小。

udydxvdt （3-5） dt，

为了解出定解u和v，附加其他约束条件，从不同的角度引入不同约束条件，从而产生不同光流分析方法。下面引入Horn-Schunck算法，即引入全局平滑性约束，假设光流在整个图像上光滑变化，即相邻的图像像素点之间具有相似的光流。图像中光流的逐像素变化可以通过速度矢量分量的空间梯度的模的平方和来定量表示[11]。简言之，其速度的变化率为零。

全局平滑约束条件的表达式方程为：

2u(u2u2 ()0 （3-6）xy

v2v2 )(0 （3-7）xy 2v(

式中u和v分别是速度矢量两个分量的空间梯度，此处假设空间和时间坐标是连续变量。

光滑约束方程的误差如下：

bfxufyvft （3-8） 偏离光滑性的误差为：

c2(u2u2v2v2)()(()xyxy （3-9） 设总误差为，则表达式为：’

2(2c2b2)dxdy （3-10） 其中，为图像的平滑程度，由实验经验得到。

将上述公式转换为拉格朗日公式，其偏微分方程如下：

22uf2xufxfyvfxft （3-11） 22vfxfyuf2yfyft （3-12） 为了减少光照变化与噪声影响，建立当前像素的一个3×3的邻域，用邻域中像素的加权值和代替当前像素的两个方向矢量[10]。

设uu、vv，则：

(2fx2fy2)(u)fx[fxfyft] （3-13） (2fxfy)(v)fy[fxfyft] （3-14） 求解u和v的值如下：

un1n2222fx[fxnfynft]

2fx2fy2

fx[fxnfynft] （3-15） vn1n2fx2fy2 （3-16）

式中的n是迭代次数。在点(x,y,t)上计算所有的偏微分。速度矢量的两个分量的初始值分别为u0(x,y,t),v0(x,y,t)通常取值为零。通过迭代运算，得到稳定值就是在点(x,y,t)处的二维运动矢量。

5.5.2模型结果及分析

图像经过灰度处理和边缘处理后，通过光流法分析与计算且运用迭代法，通过循迭代得到位移矢量的水平与竖直分量u和v。u和v的表达式如下：

un1nfx[fxnfynft]

2fx2fy2

vn1nfx[fxnfynft]

fxfy222

当u，v都为正，则图像的位移方向为右上；当u，v都为负，则图像的位移方向为左下； 当u为正，v为负，则图像的位移方向为右下；当u为负，v为正，则图像的位移方向为左上。

短时间缓慢运动图像的位移方向近似代替图像的运动方向。

六.模型的求解

数码摄像技术被广泛使用于多种场合中，有时由于客观条件的限制，设备只能在较低的分辨率下成像，现在，整个视野区域向某个方向缓慢移动，每隔一定时间拍摄一帧

图像，运动的画面体现为图像的序列，拍摄到的系列图像实时地传输到计算机中，通过上文所建立的数学模型，来判断运动的方向。

6.1.灰度变换

用matlab软件，将低分辨率下的系列图像用加权法进行灰度变换。如下图所示。

图6.1.1原图

图6.1.3原图

图6.1.5原图

图6.1.2灰度变换后的图像

图6.1.4灰度变换后的图像

图6.1.6灰度变换后的图像

图6.1.7原图

图6.1.8灰度变换后的图像

图6.1.9原图

图6.1.10灰度变换后的图像

6.2.边缘提取

用matlab软件，调用canny算子，编写程序，如附件二所示，提取灰度变换图像中的边缘。如下图所示：

图6.2.1灰度变换后的图像

图6.2.2边缘提取后的图像

图6.2.3灰度变换后的图像

图6.2.5灰度变换后的图像

图6.2.7灰度变换后的图像

图6.2.9灰度变换后的图像

6.3.光流法

图6.2.4边缘提取后的图像

图6.2.6边缘提取后的图像

图6.2.8边缘提取后的图像

图6.2.10边缘提取后的图像

经过MATLAB代码调试后，（程序如附录四所示）得到位移矢量的水平及竖直分量u和v的值及相应位移变化图。如下图所示：

图6.3.1位移变化图 图6.3.2位移变化图

图6.3.3位移变化图 图6.3.4位移变化图

由程序运行结果可知，图像运动方向为左下，则视野方向与图像运动方向相反，即为右上，与实际相符合。证明所建立的数学模型具有实时性和有效性。

七.模型的优缺点

7．1优点：

1．所采用的Canny算子具有较好的抗噪性能和较高的边缘定位精度，有利于图像的边缘提取。

2．所采用的光流法携带大量的运动信息，支持各类行摄像机，且当不知道场景先验知识时，可以检测出目标区域。

3．Canny算子提取边缘不丢失重要的边缘，没有虚假的边缘；实际边缘与检测到的边缘位置之间的偏差最小

缺点：

1．Canny准则是一个连续准则，也就是说在假设图像和滤波器都是一个连续函数的情形下给出的，但实际上数字图像是离散的。即无法确定在连续域所谓的最优的滤波器在离散的数字图像上还是不是最优的，具有不确定性。

2．应用光流法算法若没有相应的硬件支持时，它的计算耗时、检测实时性差。

3．Canny检子依赖物体的大小且对噪声敏感。

八.参考文献

[1]卞琼.基于worldview2的高分辨率遥感[J].卫星影像河流提取研究，2013.05：30

[2]胡进刚.一种面向对象的高分辨影像道路提取方法[J].遥感技术与应用，2006,21（3）：184-188

[3]贾艳丽.基于视频图像序列的运动目标检测与跟踪[D].硕士学位论文，哈尔滨工程大学2012.03：07

[4]魏玉强，运动目标检测与跟踪算法研究[D].燕山大学，2010:30-31

[5]林卉，赵长胜，舒宁，基于Canny算子的边缘检测及评价[J].黑龙江工程学院学报，第17卷第2期，2003.06：04

[6]连静.基于多尺度融合的图像边缘检测[D].硕士学位论文，吉林大学，2005.05，33-39

[7]徐辛超，付晨，徐爱功.一种改进的Canny边缘提取方法[J].遥感信息.2016.10

[8]Beauchemin S S ,Barron J L.The computation of optical flow,ACM Computing Surveys[J],1995,27(3):433-467

[9]裴巧娜.基于光流法的运动目标检测与跟踪技术[J].计算机应用技术.2009.05:07-08

[10]金国伟，基于混合高斯建模的目标检测和光流法的方向判断[D]，硕士学位论文，杭州电子科技大学，2014.12,21-23

[11]贾海涛，运动目标检测与识别算法的研究[D]，硕士学位论文，电子科技大学，2006.01,7-8,22-23

九.附件

附件一：图像灰度变换

f=imread('1-1.jpg');%读取图像,图像要保存在matlab所在的路径中,根据所要灰度处理的照片，读取对应的照片。

figure(1);

imshow(f); %显示图像

[n m a]=size(f);%判断图像的大小

for x=1:n %通过双循环对图像进行灰度化处理

for y=1:m

p(x,y)=0.3*f(x,y,1)+0.59*f(x,y,2)+0.11*f(x,y,3);

end

end

figure(2);

imshow(p);%显示处理后的图像

imwrite(p,'abc.jpg'); %保存处理后的图像

附件二：边缘提取

f=imread('1-2.jpg');% 根据所要提取的图片，读取对应的图像

imshow(f)

figure

bw=edge(f,'canny');

imshow(bw)

附件三：3*3矩阵

附件四：光流法程序

程序一：

I1=rgb2gray(imread('44.jpg'));

I2=rgb2gray(imread('55.jpg'));

imshow(I1)

figure

imshow(I2)

figure

im1=single(I1);

im2=single(I2);

numLevels=3;

window=9;

iterations=1;

alpha = 0.001;

hw = floor(window/2);

t0=clock;

pyramid1 = im1;

pyramid2 = im2;

for i=2:numLevels

im1 = impyramid(im1, 'reduce');

im2 = impyramid(im2, 'reduce');

pyramid1(1:size(im1,1), 1:size(im1,2), i) = im1;

pyramid2(1:size(im2,1), 1:size(im2,2), i) = im2;

end;

for p = 1:numLevels

im1 = pyramid1(1:(size(pyramid1,1)/(2^(numLevels - p))),

1:(size(pyramid1,2)/(2^(numLevels - p))), (numLevels - p)+1);

im2 = pyramid2(1:(size(pyramid2,1)/(2^(numLevels - p))),

1:(size(pyramid2,2)/(2^(numLevels - p))), (numLevels - p)+1);

if p==1

u=zeros(size(im1));

v=zeros(size(im1));

else

u = 2 * imresize(u,size(u)*2,'bilinear');

v = 2 * imresize(v,size(v)*2,'bilinear');

end

for r = 1:iterations

u=round(u);

v=round(v);

for i = 1+hw:size(im1,1)-hw

for j = 1+hw:size(im2,2)-hw

patch1 = im1(i-hw:i+hw, j-hw:j+hw);

lr = i-hw+v(i,j);

hr = i+hw+v(i,j);

lc = j-hw+u(i,j);

hc = j+hw+u(i,j);

if (lr size(im1,1))||(lc size(im1,2)) else

patch2 = im2(lr:hr, lc:hc);

fx = conv2(patch1, 0.25* [-1 1; -1 1]) + conv2(patch2, 0.25*[-1 1; -1 1]);

fy = conv2(patch1, 0.25* [-1 -1; 1 1]) + conv2(patch2, 0.25*[-1 -1; 1 1]);

ft = conv2(patch1, 0.25*ones(2)) + conv2(patch2, -0.25*ones(2));

Fx = fx(2:window-1,2:window-1)';

Fy = fy(2:window-1,2:window-1)';

Ft = ft(2:window-1,2:window-1)';

A = [Fx(:) Fy(:)];

G=A'*A;

G(1,1)=G(1,1)+alpha; G(2,2)=G(2,2)+alpha;

U=1/(G(1,1)*G(2,2)-G(1,2)*G(2,1))*[G(2,2) -G(1,2);-G(2,1)

G(1,1)]*A'*-Ft(:);

u(i,j)=u(i,j)+U(1); v(i,j)=v(i,j)+U(2);

end

end

end

end

etime(clock,t0)

end

u=u(window:size(u,1)-window+1,window:size(u,2)-window+1)

v=v(window:size(v,1)-window+1,window:size(v,2)-window+1)

RGB1=showmap3(u,v,5);

imshow(RGB1);

程序二：

function RGB = showmap3(u,v,vmax)

anglestep = pi/3;

colorstep=1/anglestep;

RGBpart=[ 0 1 255 255 -1 0;

-1 0 0 1 255 255;

255 255 -1 0 0 1];

RGB=zeros(size(u,1),size(u,2),3);

for i=1:size(u,1)

for j=1:size(u,2)

angle=atan2(u(i,j),v(i,j));

if (angle

angle = angle+2*pi;

end

part=angle/anglestep;

if ( part >= 6.0 )

part = 0.0;

end

length=sqrt( u(i,j)*u(i,j) + v(i,j)*v(i,j) );

if(length>vmax)

RGB(i,j,1)=1;

RGB(i,j,2)=1;

RGB(i,j,3)=1;

else

for k=1:3

switch RGBpart(k,fix(part+1))

case 0

RGB(i,j,k)=0;

case 1

RGB(i,j,k)=colorstep*mod(part,1)*length/vmax;

case -1

RGB(i,j,k)=(1-colorstep*mod(part,1))*length/vmax;

case 255

RGB(i,j,k)=length/vmax;

end

end

end

end

end

end