2008 高教社杯全国大学生数学建模竞赛
数码相机定位
摘要
本文假设数码相机成像原理为小孔成像,在此基础上,通过两种合理的模型
对数码相机定位问题进行了较深入的研究。
针对问题一和二,我们建立了两种不同模型——变换矩阵模型和公切线模
型。在变换矩阵模型中,建立了物、像、相机三个坐标系,分别称为世界坐标系,
像坐标系和光心坐标系。研究世界坐标系向像坐标系的变换矩阵3 4 ( ) ij M a ´ = ,推
导出圆在像坐标系中的像为椭圆。利用灰度检测可以得到像中各椭圆圆周上各点
的坐标,通过多元线性回归拟合出各椭圆方程;对单独一个圆进行研究时,在合
理的近似前提下,以圆心为世界坐标系的原点,可求出该圆心所成像在像坐标系
中的坐标u = a14 , v = a24。最后我们求得5个圆心所成像在光心坐标系中的坐标分
别为(单位:mm):(-50.00,51.32,-417.20)、(-23.54,49.47,-417.20)、(33.86,
45.24,-417.20)、(-60.05,-31.22,-417.20)、(18.52,-31.48,-417.20)。在公
切线模型中,通过简单几何证明,得出在小孔成像时,公切线交点的像就是公切
线像的交点,联系题目中所给标靶的特殊性(所有圆全等),得出像平面中公切
线交点连线的交点就是标靶中对应圆心的像,并设计了一种算法得到5个圆心所
成像在光心坐标系中的坐标分别为(单位:mm):(-49.92,51.36,-417.20)、
(-23.47,49.34,-417.20)、(33.88,45.05,-417.20)、(-60.04,-31.29,-417.20)、
(18.58,-31.56,-417.20)。
在问题三中,我们用计算机模拟的方法,统计和分析了我们模型的在不同的
情况下所得到的结果与理论值之间的误差,并着重研究了相机与标靶的距离和像
平面与圆平面之间的偏角对结果的影响。结果表明在一定的前提下,当相机与标
靶的距离大于200 毫米,-0.5 £a £ 0.5以及-1£ b £ 0.5(单位为弧度)时,我们
的结果与理论值相差不到一个像素,有着较好的稳定性和精度。
问题四中,通过每个相机旋转变换矩阵R和平移向量T,可以得到两相机的
变换关系: 1 1
1 2 1 2 2 1 R = R R -、T = R R - T +T ,即相对位置关系,并理论推导了从两相
机中像在光心坐标系中的参数得到物在世界坐标中的参数,实现双目定位。另外
在相机的光心和像屏中心的连线垂直于象平面基础上,我们还给出另外一种合理
模型,通过矢量的方法求出物相对于光心坐标系的精确位置,从而可以得到两相
机的相对位置。
关键词:相机定位、小孔成像、变换矩阵、公切线、计算机模拟
数学中国论文共享
www.madio.cn
2
一、问题重述
数码相机定位在交通监管(电子警察)等方面有广泛的应用。所谓数码相机
定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的
定位方法是双目定位,即用两部相机来定位。对物体上一个特征点,用两部固定
于不同位置的相机摄得物体的像,分别获得该点在两部相机像平面上的坐标。只
要知道两部相机精确的相对位置,就可用几何的方法得到该特征点在固定一部相
机的坐标系中的坐标,即确定了特征点的位置。于是对双目定位,精确地确定两
部相机的相对位置就是关键,这一过程称为系统标定。
标定的一种做法是:在一块平板上画若干个点,同时用这两部相机照相,分
别得到这些点在它们像平面上的像点,利用这两组像点的几何关系就可以得到这
两部相机的相对位置。然而,无论在物平面或像平面上我们都无法直接得到没有
几何尺寸的“点”。实际的做法是在物平面上画若干个圆(称为靶标),它们的圆
心就是几何的点了。而像平面上的园一般会变形,所以必须从靶标上的这些圆的
像中把圆心的像精确地找到,标定就可实现。
有人设计靶标如下,取1 个边长为100mm 的正方形,分别以四个顶点(对
应为A、C、D、E)为圆心,12mm为半径作圆。以AC边上距离A点30mm处
的B 为圆心,12mm为半径作圆。用一位置固定的数码相机摄得其像。利用所得
图像, 具体解决如下几个问题:
(1) 建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标,
这里坐标系原点取在该相机的焦点,x-y平面平行于像平面;
(2) 对由图2、图3 分别给出的靶标及其像,计算靶标上圆的圆心在像平面上
的像坐标, 该相机的像距(即焦点到像平面的距离)是1577 个像素单位(1
毫米约为3.78个像素单位),相机分辨率为1024×786;
(3) 设计一种方法检验你们的模型,并对方法的精度和稳定性进行讨论;
(4) 建立用此靶标给出两部固定相机相对位置的数学模型和方法。
二、模型假设
a.本题中数码相机成像系统看成是小孔成像;
b.灰度检测前将图像改成黑白图的误差不予考虑。
三、符号说明
W O :世界坐标系
C O :光心坐标系
P O :像坐标系
R:旋转矩阵
T :平移向量
M :空间变换矩阵
数学中国论文共享
www.madio.cn
3
f :焦距(mm)
L:每毫米的像素单位
P:图像矩阵
k B :第k个圆的边缘点集
四、模型的建立与求解
引理:
我们针对小孔模型提出一些内在性质并给予简单证明。
图 1
性质1:直线经小孔后所得到的像仍为直线。
证明:如图,a、g 分别为物平面和像平面,b 为与像平面平行的平面,点O为
小孔,AC为a 上的线段,B为AC上任一点,经小孔O后得g 上直线2 2 A C 。由
AC和O 可以确定一个平面OAC,所以2 2 AC 为平面OAC和平面g的交线;设B
经小孔成像到g上的2 B 点。因为2 B 同时在平面OAC 和平面g上,所以2 B 必在
两平面交线即2 2 AC 上。因为B为AC上任一点,所以直线AC经小孔O成像为g
平面上的直线2 2 AC ,故直线经小孔后所得到的像仍为直线。
性质 2:线段中点经小孔成像所得点不一定还是像线段(经小孔成像后所得到的
线段)的中点。
证明:如图, 1 1 AC 为平面OAC与b平面的交线, 1 B 为B所对应点。因为b 平行
于g,平面OAC分别交b、g与1 1 AC 、2 2 AC ,则1 1 AC 与2 2 AC 相互平行。
当a 与b不重合即q不为0 时,如B 为AC中点,因为1 1 CC与AA 相交,则1 B必然
数学中国论文共享
www.madio.cn
4
不是 1 1 AC 中点,所以2 B 也必然不是2 2 AC 中点,也就是说线段中点经小孔成像所
得点不一定还是像线段的中点,仅当该线段平行于像平面时才仍是中点。
性质 3:两直线交点的像仍是两直线像的交点。
证明:如下图:a 平面上两直线AB与CD交于点E,经小孔O两直线分别得到b 平
面上的像A1B1和C1D1,点E在b 平面上的像为1 E ,则平面1 1 1 1 ABB A与CDDC 交于
直线OE,平面1 1 1 1 ABB A与CDDC 分别交平面b 于直线1 1 1 1 A B与C D ,则1 E
同时在平
面1 1 1 1 ABB A、CDDC与b 上,则1 E 同时在1 1 1 1 A B和C D 上,即b 平面上, 1 1
1 1 A B和C D
交于1 E 点,两直线交点的像仍是两直线像的交点。
我们将此结论推广到曲线的情形,很显然结论也必然成立。只要在曲线相交处对
两曲线取无限小段,那么就可以看成是直线相交的情形。
图二
性质 4:圆经小孔成像为椭圆.
图三
数学中国论文共享
www.madio.cn
5
证明:如图所示,圆O经小孔成像如图。在一个与圆O所在平面平行的平面上所
成像仍然是圆,如图中圆2 O (这点很容易通过以上三个性质得到)。那么对于在
与圆O所在平面不平行的平面上的像,则可以看成是用一个平面去截图中的圆锥
(当然图中圆锥可以无限的延长),根据圆锥曲线的定义可知,所截出来的图形
就是椭圆。且可以看出圆心O在截面上所成的像1 O 并不是椭圆中心M (只有截
面与圆平面平行时才是椭圆中心),所以圆经小孔成像为椭圆得证。
性质 5:圆的某一条切线的切点的像,仍然是椭圆的切线,而且切点的像就是椭
圆的切点。
证明:由性质3 的推广可知,任意两条曲线的交点的像就是他们两个像的交点。
因为圆中切线与圆是相交于一点的,那么像中切线的像与圆的像(椭圆)至
少也会有个交点。假设圆的切线的像不再是像当中椭圆的切线,则切线的像与椭
圆必有两个交点。根据光路可逆原理,我们可以把原来的圆看成是像(椭圆)经
小孔所得到的像。那么对于另外一个交点,它的原像也必然是原像平面上两曲线
的交点,则原像平面中处切点外还有另外一交点,产生矛盾,故假设不成立,所
以结论得证。
模型一:变换矩阵模型
问题一
在标靶上,以某个圆的圆心为原点建立空间直角坐标系,由W W X ,Y , W Z 轴组
成,称其为世界坐标系;在像空间上建立像坐标系,由u、v轴组成;由于摄像
机可以安放在环境中的任何一个位置,我们也建立一个坐标系来描述,由
C C C X 、Y 、Z 轴组成,原点位于光心,称其为光心坐标系。光心坐标系与世界坐
标系间的转换可以同过旋转矩阵R 和平移矩阵T 来实现。如空间一点P 在世界坐
标系和光心坐标系中的坐标分别为( ) ( )T T
W W W C C C X Y Z 、X
Y Z ,于是存在关系:
[ ] [ ] C 1 [ 1] 1 1 ...........(0)
0 1
T T T
C C T W W W W W W
R T
X Y Z X YZ M X Y Z
é ù
= ê ú
=
ë û
其中R 为3×3的正交单位矩阵,0 (0 0 0)T = T, 1 M 为4×4 矩阵,0T和1的
加入只是为了方便以后的计算。空间点p的像在像坐标系的位置与p在光心坐标
系中的关系如图可得:
C , C
C C
x fX y fY
Z Z
= =
其中(x, y)为点p的像在像坐标系的坐标,写成矩阵形式就是:
数学中国论文共享
www.madio.cn
6
0 0 0
0 0 0 ....................................(1)
1 0 0 1 0
1
C
C
C
C
X
x f
Y
Z y f
Z
é ù
é ù é ùê ú
ê ú = ê ú ê ú ê ú ê úê ú
êë úû êë úû ê ú
ë û
所以可以得到:
0 0 0
,
0 0 0 ...............................(2)
0 ,1
1 0 0 1 0
1
W
W
C T
W
X
x f
R T Y
Z y f
Z
é ù
é ù é ù ê ú ê ú ê ú é ù ê ú ê ú = ê ú ê úë û ê ú êë úû êë úû ê
ú
ë û
参照附录中对于旋转矩阵的说明,我们可以得到:
sin
sin cos sin sin cos
cos cos
W
C W Z W W Z
W
X
Z Y T X Y T
Z
b
a b b a b
a b
é ùé ù
= ê- ú ê ú - = - - ê úê ú
êë úû êë úû
所以由:
0
0 0 0
,
0 0 0
0 ,1
1 0 0 1 0
1 1
W W
W W
C T
W W
X X
x f
R T Y Y
Z y f M
Z Z
é ù é ù
é ù é ù ê ú ê ú ê ú ê ú é ù ê ú ê ú ê ú = ê ú ê ú = ë û ê ú ê
ú êë úû êë úû ê ú ê ú
ë û ë û
可得:
11 12 14
21 22 24
( sin sin cos )
( sin sin cos )
W W Z W W
W W Z W W
X Y T u a X a Y a
X Y T v a X a Y a
b a b
b a b
- - = + +
- - = + +
可以通过上式解得:
22 14 14 22 12 24 12 24
22 11 11 22 12 21 21 12
21 14 14 21 11 24 11 24
22 11 1
( sin cos ) ( sin cos )
( sin sin cos ) ( sin sin cos )
( sin ) ( sin )
( sin sin cos
Z Z
W
Z Z
W
X uT a va a a vT a ua a a
ua va aa va ua a a
Y uT a va a a vT a ua a a
ua va a
a b a b
b a b b a b
b b
b a b
+ + - + +
=
- - - - -
- + - -- + -
=
- -1 22 12 21 21 12 a ) - (va sin b -ua sina cos b - a a
)
因为W W X Y 、为圆周上的点,所以在世界坐标系满足:
2 2 2
W W X +Y = r
带入可得到二次曲线方程:
2
22 24 12 14 14 22 12 24
2
21 24 14 11 14 21 11 24
2 2
22 21 11 12 11 22 21 12
[( sin cos ) ( sin cos ) ( )]
[ ( sin ) ( sin ) ( )]
[( sin sin cos ) ( sin cos sin ) ( )]
Z Z
Z Z
Ta a u Ta a v a a a a
Ta a u a T a v a a a a
r a a u a a v a a a a
a b a b
b b
b a b a b b
- - - + +
+- + + + - +
= + - + - -
其代表一椭圆。
注意到此处的C Z 是标靶上圆上点在光心坐标系中C Z 方向上的坐标。对于本
数学中国论文共享
www.madio.cn
7
题来说,因为标靶在光心坐标系中C Z 方向上的坐标应该大于相机的两倍焦距,
即应在1 米附近;而对于标靶上一个圆,其半径仅为12mm,当标靶平面与光心
坐标系中C C C X O Y 平面存在一定夹角q时,那么一个圆上所有点在C Z 上坐标的差
异只是12sinq ,与1m 相比较而言,其误差是比较小的,所以我们可以近似认为
对于标靶上同一个圆上的点,其C Z 是相同的(后面我们将来讨论这种近似所带
来的误差,会发现其误差是非常小的,可见这种近似的合理性),
所以在此情况下,我们认为C Z 是不变的值,于是有:
0 0 0
1 , 0 0 0 .........................(3)
0 ,1
1 0 0 1 0
1
W
W
T
C W
X
x f
R T Y
y f
Z Z
é ù
é ù é ù ê ú ê ú ê ú é ù ê ú ê ú = ê ú ê úë û ê ú ëê ûú ëê ûú ê
ú
ë û
由于在处理像平面上的点时,我们经常用的单位为像素,对像平面坐标为
(x, y)的点,改用像素为单位后,其坐标为:
u xL
v yL
= ìí
î =
其中L 为每毫米的像素单位。于是:
0 0 0 0 0 0 0
1 , 0 0 0 0 0 0 0 ...(5)
0 ,1
1 0 0 11 0 0 1 0 0 1 0
1 1
W W
W W
T
C W W
X X
u L x L f
R TY Y
v L y L f M
Z Z Z
é ù é ù
é ù é ù é ù é ù é ù ê ú ê ú ê ú ê ú ê ú ê ú ê ú é ù ê ú ê ú ê
ú = ê ú ê ú = ê ú ê ú ê ú = ë û ê ú ê ú êë úû êë úû êë úû êë úû êë
úû ê ú ê ú
ë û ë û
其中M 为一3×4 的坐标变换矩阵,而且在相机本身内部参数和相对标靶的
位置不变时,M 是一个常数矩阵。由上式可以看出对于世界坐标系中的任意一
点( , , ) W W W X Y Z ,经坐标变换矩阵M 变换,便可以得到其像在像坐标系中的坐标
(u, v)。
针对题目所给信息,为简化模型,我们可分别随标靶中每个圆分别讨论(因
为对每个圆讨论的方法完全相同,故本文只详细讨论一个圆,其他的可完全类
比)。对某个圆讨论时,取该圆所处坐标系为世界坐标系,坐标原点取在圆心处
(后面会发现这样选取的精妙之处),所得像处于像坐标系,数码相机处于摄像
坐标系。可以很直觉的发现圆周上的点是至关重要的,那么我们先来探讨圆周上
的各点及其所成像的位置。
因为在圆周上,注意到此时的世界坐标系里只需其二维情形,取ZW=0,故其
曲线方程为:
2 2 2
W W X +Y = r
数学中国论文共享
www.madio.cn
8
因为变换矩阵M 为3×4 矩阵,可设
11 12 13 14
21 22 23 24
31 32 33 34
a a a a
M a a a a
a a a a
é ù
= ê ú ê ú
êë úû
所以圆周上任意一点( , , ) W W W X Y Z ,经坐标变换矩阵M 变换后得到像在像坐
标系的坐标:
11 12 13 14
21 22 23 24
31 32 33 34
..............(6)
0 0
1
1 1
W W
W W
X X
u a a a a
Y Y
v M a a a a
a a a a
é ù é ù
é ù ê ú é ù ê ú
ê ú = ê ú = ê ú ê ú ê ú ê ú ê ú ê ú
êë úû ê ú êë úû ê ú
ë û ë û
有:
11 12 14
21 22 24
W W ...................(7)
W W
u a X a Y a
v a X a Y a
= + + ìí
î = + +
与圆的曲线方程联立可得像的曲线方程:
2 2
22 14 12 24 12 14 11 24 2
11 22 21 12 12 21 22 11
a (u a ) a (v a ) a (u a ) a (v a ) r
aa a a a a a a
é - - - ù é - - - ù
ê ú + ê ú = ë - û ë - û
从上式不能很显然的发现什么性质,所以我们将其展开:
2 2 2 2 2 2 2 2
21 22 11 12 11 21 12 22 14 21 22 24 11 21 12 22
2 2 2 2 2 2 2 2
24 11 12 14 11 21 12 22 21 22 14 11 12 24 11 21
2 2
12 22 14 24 12 21 22 11
( ) ( ) 2( ) [2 ( ) ( )]
2[ ( ) ( )] [( ) ( ) 2(
) ( ) ] 0
a a u a a v a a a a uv a a a a aa a a u
a a a a a a a a v a a a a a a a a
a a aa a a a a r
+ + + - + - + + +
- + + + + + + + - +
- - = ......(8)
简化上式为:
2 2
1 2 3 4 5 6 k + k v + k u + k uv + k v + k u =
0..........................(9)
可见这是一个一般二次曲线方程,故其图像为一椭圆(当然也可能会是一条
直线,暂不考虑这种情况)。其中1 k 至6 k 分别为上式的系数。观察(8)式可以发
现这个重要信息:u、v都是14 24 a 、a 和u2、v2、uv的系数组成,而14 24 a 、a 正是世
界坐标系的圆心经转换矩阵得到的像坐标系的坐标。所以其简化式可变为:
2 2
1 24 5 14 4 14 6 24 4 4 5 6 k - (2a k + a k )v - (2a k + a k
)u + k uv + k v + k u = 0.......(10)
我们所需要求的是圆心在像坐标系的像坐标的坐标,由于我们将世界坐标系
的原点放在圆心,故圆心在世界坐标系中的坐标为(0, 0),根据(7)式,在像坐
标系的像坐标的坐标为:
14
24
o
o
u a
v a
= ìí
î =
数学中国论文共享
www.madio.cn
9
于是我们的问题转化为求解a14、a24,而由我们(10)式的分析,我们可以
得到:
24 5 144 2
14 6 244 3
(2 )
.............(11)
(2 )
a k a k k
a k a k k
- + = ìí
î- + =
上式为关于14 24 a ,a 的二元一次方程组,于是问题又转化为求出(2 6) ik £ i £ ,
我们知道(9)式是圆的像在像平面的曲线方程,而i k 为方程的系数,下面我们
来求解这个方程。
我们的数据来源于数码相机拍摄的照片,所以我们首先需要从照片中获取数
据。我们知道,现代的数码照片的存储原理基于点阵。例如在一张灰度照片中存
储了n´m的矩阵P,矩阵中的每个元素P(i, j)记录的是一个灰度值,表示这张
照片在该点的明暗程度(而在一张彩色照片里每个点存储的是一组RGB)。所以
我们的直接数据是一个768´1024的矩阵,下面我们需要在这个矩阵中提取圆的
像的边缘以便求解该圆的像的曲线方程。
传统的图像边缘检测的算法有基于灰度直方图的边缘检测、基于梯度的边缘
检测、Laplacan 边缘算子、Canny 边缘检测算子、模糊推理的边缘检测和
Mallat小波边缘检测算子等,不同的算法有着不一样的应用背景,而在本文中
所需处理的照片经过简单的处理后可以转化为单色的照片,所以用最简单的基于
灰度直方图的边缘检测算法便可得出较理想的效果。
图四 经过黑白处理后的照片
经过黑白处理后,所有点的灰度值只有0(黑色)和255(白色)两种,我
们用下面的算法进行边缘检测(对圆A、B、C、E、D依次编号为圆1、2、3、4、
5):
边缘检测算法:
Step 1 设灰度阀值T,当某点的灰度小于T 时,表示其为黑色,反之则
为白色。在本文中我们取T=128;
数学中国论文共享
www.madio.cn
10
Step 2 为每个圆的像估算一个矩形邻域,使得该矩形邻域可以完全包含
该圆的像且仅包含该圆的像;
Step 3 对第k 个圆的像,执行第4-6 步(k=1,2,3,4,5);
Step 4 对矩形邻域进行逐行扫面,当出现P(i -1, j) ³ T 且P(i, j) <
T
或者P(i, j) < T 且P(i +1, j) ³ T 时将P(i,
j)加入该圆的像的边缘
点集k B ;
Step 5 对矩形邻域进行逐列扫面,当出现P(i, j -1) ³ T 且P(i, j) <
T
或者P(i, j) < T且P(i, j +1) ³ T 时将P(i,
j)加入该圆的像的边缘
点集k B ;
Step 6 除去k B 中重复的点。
我们得到的点集{( , )| ( ,) i } k i i i i B = u v iÎZ +且u v
位与第个圆的像的边缘,这里
, i i u v 的单位为像素,像坐标系以左上角为原点,图片的正下方为u轴,正右方为
v轴。
对于每一个圆,我们将其所有边缘点带入其曲线方程:
2 2
1 2 3 4 5 6 0,1 | | .........(12) i i ii i i k k + v k + u k +
u v k + v k + u k = £ i £
B
我们得到了关于(1 6) ik £ i £ 的线性齐次方程组,而且在这个方程组中,方程
的个数远远大于未知数的个数,而且我们发现若直接求解往往只能得到零解,于
是我们将其化为多元线性回归模型求解。
由上文讨论,像的曲线为椭圆,故平方项系数不为零,故原方程可化为:
1 2 3 4 2 5 2
6 6 6 6 6
0,1 | | ..........(13) i i i i i i k
k v k u k u v k v k u i B
k k k k k
+ + + + + = £ £
令
6
j 1 5
j
k
b j
k= , £
£
,
则
有
:
2 2
1 2 3 4 5,1 | | ...........(14) i i i ii i k -u = b + v b + u
b + u v b + v b £ i £ B
上式是典型的多元线性回归模型,应用最小二乘法可以求解。当我们解出了
1 2 3 4 5 b ,b ,b ,b ,b 后,由(11)式得:
24 5 144 2
14 244 3
(2 )
(2 )
a b a b b
a a b b
- + = ìí
î- + =
数学中国论文共享
www.madio.cn
11
v
u
于是我们得到了圆心的像在uov平面上的坐标:
3 5 2 4
14 2
4 5
2 3 4
24 2
4 5
2
4
.......................(15)
2
4
o
o
u a b b b b
b b
v a b b b
b b
ì - ï = = - ïí
ï - = =
ïî -
最后,我们需要将坐标转换到题目中给定的坐标系中,(x 与v 同向,y
与u 反向)即:
( )
( 512) /
(384 )/ .......................... 16
/
o o
o o
o
x v L
y u L
z f L
= - ìï
= - íï
î = -
问题二
首先,我们用matlab 实现边缘检测算法,我们得到了每个圆的像的边缘点
集k B ,如下图:
图五
然后我们对每个圆的像的边缘点用多元线性回归去求该圆的像的曲线方程
(置信度取0.05),如下表所示:
表1 第1个圆的像的曲线方程参数
参数 参数估计值参数置信区间
b1 140397.20023 [139748.34930 141046.05115]
b2 -640.44473 [-644.04310 -636.84637]
b3 -407.58657 [ -409.32833 -405.84480]
b4 0.08784 [ 0.08246 0.09323]
b5 0.96592 [ 0.96064 0.97120]
R2=0.99999393 , F= 11279037.31254 , p=0
数学中国论文共享
www.madio.cn
12
表2 第2个圆的像的曲线方程参数
参数 参数估计值参数置信区间
b1 227111.11949 [ 225858.66719 228363.57179]
b2 -869.84886 [-875.31270 -864.38502]
b3 -453.55647 [-456.21530 -450.89765]
b4 0.14049 [ 0.13421 0.14677]
b5 0.99545 [ 0.98926 1.00165]
R2=0.999993179 , F=9786431.732757331 , p=0
表3 第3个圆的像的曲线方程参数
参数 参数估计值参数置信区间
b1 521021.19332 [ 517789.67136 524252.71528]
b2 -1431.47741 [ -1441.06500 -1421.88982]
b3 -603.80738 [-608.27782 -599.33694]
b4 0.27703 [ 0.27005 0.28402]
b5 1.07232 [1.06508 1.07957]
R2=0.99999442 , F=11019117.97897741 , p=0
表4 第4个圆的像的曲线方程参数
参数 参数估计值参数置信区间
b1 343623.44295 [ 342627.40480 344619.48111]
b2 -557.42512 [-561.92415 -552.92609]
b3 -1057.92607 [-1059.63303 -1056.21910]
b4 0.19027 [0.18428 0.19625]
b5 0.81131 [0.80593 0.81670]
R2=0.99999918 , F=75335174.25731819 , p=0
表5 第5个圆的像的曲线方程参数
参数 参数估计值参数置信区间
b1 660749.59267 [ 657278.55376 6.64220.63157]
b2 -1228.17779 [-1237.71197 -1218.64362]
b3 -120816690 [-1212.41769 -1203.91612]
b4 0.34664 [0.33934 0.35393]
b5 0.90414 [0.89715 0.91114]
R2=0.99999914 , F=65011077.78295716 , p=0
在上表最后一行中,R2是回归方程的决定系数,F 为F 统计量值,p 为与F
统计量对应的概率值。从数据看来,R2都达到了99.999%以上,表明-u2的99.999%
以上的可以由我们求出来的方程确定,F 值远远超过F 检验的临界值,p 更是远
小于置信度0.05,所以我们求出的曲线方程是高度精确的。
将以上数据代入(15)式,我们便得到了五个圆心的像的坐标:
数学中国论文共享
www.madio.cn
13
表6 圆心像的坐标
圆 o u o v ( ) o u 像素o v(像素)
1 189.61023 322.89682 190 323
2 197.06280 423.00273 197 423
3 213.26380 639.91289 213 640
4 501.87983 284.67960 502 285
5 503.07979 582.75218 503 582
用(16)式将其转换到题目所给的坐标系:
表7 光心坐标系下圆心像的坐标
圆 ( ) o xmm ( ) o ymm ( ) o zmm
1 -50.00000 51.32275 -417.19577
2 -23.54497 49.47089 -417.19577
3 33.86243 45.23809 -417.19577
4 -60.05291 -31.21693 -417.19577
5 18.51851 -31.48148 -417.19577
以上我们就求解出标靶上各圆圆心的像在像坐标系下的具体坐标。
模型二:公切线模型
图六图七
如图所示,为标靶平面与像平面的示意图。图中直线为公切线。物平面内,圆1 O
与O2、O3 的公切点分别为B、D和A、C ,对应像平面中为椭圆1 O '与椭圆
2 3 O '、O '的公切点分别为D'、B'和C'、A',由本文所给出的性质可以很容易得出
D'、B'、C'、A'分别为D、B、C、A的像点。
数学中国论文共享
www.madio.cn
14
图八图九
在标靶平面,显然四边形PQRS为正方形,则PR与QS连线交于圆心1 O ,由引理
中性质3 可知P'Q'R'S '分别为PQRS的像,则P'R'与Q'S '也分别为PR与QS的
像,则圆心1 O 在像平面所对应的像为P'R'与Q'S '的交点1 O '(注意这里1 O '并不
一定就是该椭圆中心)。
两椭圆公切线求解算法如下:
Step 1:由模型一中灰度检测,我们已经得到像平面上椭圆圆周上各点坐标;
Step 2:在两椭圆上各任取一点,连成直线l,固定其中一点(静点),扫描另一
点(动点)所在椭圆的圆周上各点,每扫描到一点在直线l上方(或下方,对应
的是另一条切线),用该点取代动点成为新的动点,此时直线l也变成静点和此新
的动点的连线,扫描一周后停止。
Step 3:固定此时的动点(这时就称之为静点),扫描另一椭圆圆周上各点,当该
点在直线l上方时,取代动点成为新的动点,此时直线l也变成静点和新的动点
的连线,扫描一周后停止;
Step 4:循环step 2 和step 3,直至两椭圆圆周上都不再存在位于直线l上方的点,
此时的l就是两椭圆的公切线。
对该算法而言,收敛速度很快,一次迭代后就几乎接近了切点,所以可以用
来高效地计算切点。
我们按照上述思想编写了程序,所得的切线图如下:
图十
数学中国论文共享
www.madio.cn
15
在我们的程序中,我们还计算了切线之间的交点及其交点连线的交点,即上
文所分析的圆心的像,结果如下:
表8 圆心像的坐标
圆 o u o v ( ) o u 像素o v(像素)
1 189.87734 323.30371 190 323
2 197.51240 423.29359 198 423
3 213.70159 640.08497 214 640
4 502.26962 285.03296 502 285
5 503.31215 582.23246 503 582
表9 光心坐标系下圆心像的坐标
圆 ( ) oxmm ( ) o ymm ( ) o zmm
1 -49.91965 51.35520 -417.19577
2 -23.46730 49.33535 -417.19577
3 33.88491 45.052489 -417.19577
4 -60.04419 -31.28799 -417.19577
5 18.58001 -31.56405 -417.19577
问题三
在这个部分中,我们针对我们在问题1,2 中的方法可能带来的误差及稳定
性提出了一种基于计算机模拟的检验方法。
在我们的前面的推导中,我们知道任意一个点在世界坐标系上的点的坐标和
其像在像坐标系的点的坐标满足(0)式,在模型中我们为了简化计算把c Z 作为
了一个常量来计算(而事实上当拍摄距离较远,标靶上的圆又较小时,这种假设
是合理的),这带来了我们的第一个误差;我们的模型中的另外两个个误差来自
于边缘检测的精确度和用多元线性回归来计算像的曲线方程所带来的误差,而实
际上边缘检测的误差实际上来自于数字图像本身(由于一张图片只能存储有限的
点,所以将具体的像点映射到像素矩阵时会带来误差,但不会超过一个像素),
而从问题二的结果中我们可以得知多元线性回归求曲线方程有着很高的精度,所
以这两个误差相对于我们的第一个误差影响非常的小。
我们知道,(2)式在我们的针孔模型中是精确成立的,如果我们知道了相机
和标靶的确切的位置关系,那么对于标靶上的任意一点,我们总可以根据(2)
式计算出该点的像在像坐标系中的位置,也就是说在这个前提下,圆心的像在像
坐标系中的位置可以精确计算出来,这就是我们这个模型中的理论解。我们的目
的就是分析我们模型的解和理论解之间的关系。
在提出检验方法之前,我们先提出下面的引理。
引理 1:若世界坐标系上按照旋转矩阵3 3 R ´ 进行旋转,然后按向量
数学中国论文共享
www.madio.cn
16
T = (XT ,YT ,ZT )进行平移得到像坐标系,那么原来在世界坐标系中的一点
( , , ) W W W p X Y Z 在光心坐标系中的坐标的Z 轴分量C Z 满足:
C 31 W 32 W 33 W T Z = R X + R Y + R Z + Z
证明:直接由(0)式化简可得。
在我们下面的方法中,我们的世界坐标系W O 的原点在半径为12mm 的圆的圆
心上, W W X OY 平面与圆所在平面平行,光心坐标系C O 的原点及相机的光心,
W W X OY 平面与像平面平行,将光心坐标系C O 沿主光轴平移到主光轴与像平面的
交点即得像坐标系P O ,如下图:
图十一
引理 2:圆心的像在像平面的坐标的理论值为( T , T ,0)
T T
LfX LfY
Z Z
- -
证明:由:
0 0
0 0 0 0 0 0 0 0
1 , 0 1 , 0 0 0 0 0 0 0 0 0
0 ,1 0 0 ,1 0
1 0 0 1 0 0 1 0 0 0 1 0
1 1
o
o T T
C T
u L f Lf
R T R T
v L f Lf
Z Z
é ù é ù
é ù é ù é ù ê ú é ù ê ú ê ú ê ú ê ú é ù ê ú ê ú é ù ê ú ê ú =
ê ú ê ú ê ú = ê ú ê ú ë û ê ú - ë û ê ú ëê ûú ëê úû êë úû ê ú êë úû
ê ú
ë û ë û
得:
T
o
T
T
o
T
u LfX
Z
v LfY
Z
ì = - ïïíï
= -
ïî
数学中国论文共享
www.madio.cn
17
检验方法的提出:
Step 1 选取合理的旋转矩阵R 和平移向量T ,即初始化相机与标靶的相对
位置。
Step 2 在圆周上进行等间隔采样,对于每一个采样点,其在世界坐标系的
坐标可以表示为(12cosq ,12sinq ,0),那么可以得到其像点在像坐标系的坐标为:
0 0
31 32 33
12cos
1 1 12sin
0
1
1 1
W
W
C W W W W T
X
u
Y
v M M
Z Z R X RY R Z Z
q
q
é ù é ù
é ù ê ú ê ú
ê ú = ê ú = ê ú ê ú ê ú + + - ê ú
ëê ûú ê ú ê ú
ë û ë û
Step 3 经过第二步后我们得到了圆的像的边缘点集,这些点是精确的,也
是我们的模型中可以从图片中准确获取的,所以我们将这个边缘点集作为我们模
型的输入,然后利用该边缘点集按照我们的模型求解出圆心的像在像坐标中的坐
标(u, v,0)。
Step 4 计算(u,v,0)与理论值( , ,0) o o u v 的误差D(这里用欧氏距离计算)。
Step 5 重复执行1到4步,统计误差数据。
在用计算机进行模拟的时候,我们特别关心两个问题,一是相机与标靶的距
离对结果的影响;二是像平面与圆平面之间的偏角(由旋转矩阵R决定)对结果
的影响。
接下来我们来研究R和T 的合理选取的问题。
由上文的分析可知,旋转矩阵R中的三个参数a,b,l分别代表y轴向z 轴旋
转的角度,z 轴向x轴旋转的角度,x轴向y轴旋转的角度,例如下图:
图十二
我们几乎不用考虑l ¹ 0的情况,因为这实际上是将相机左右倾斜,这并不
影响拍出的物体的形状和大小,而只会影响拍出的物体在整张照片中的位置。而
数学中国论文共享
www.madio.cn
18
a ¹ 0对应的情景时我们在物体的上(下)前方将镜头对着其进行拍摄,b ¹ 0对
应的情景时我们在物体左(右)前方将镜头对着其进行拍摄,a, b 体取值范围为
[ , ]
2 2
p p
- 。
对 ( , , ) T T T T = X Y Z ,由于圆心在光心坐标系的坐标即( , , ) T T T -X -Y
-Z ,而圆必
须在相机的前方,故0 T -Z > ,即有0 T Z <
。除此之外,还需考虑T 的长度,即
相机与标靶的距离,从我们前一问的数据中,可以估算出这段距离大约在500
毫米左右。
基于以上分析,下面我们用计算机模拟的方式对我们的模型进行检验。
1.相机与标靶的距离对结果的影响
令a = g = 0,b = -0.7,T = (-k,0, -k),1£ k £
1000,我们对每一个k对应的
相机和标靶的相对位置计算出我们模型的解与理论解的误差。
实验结果表明,当1£ k £ 20,即相机与标靶的距离大约小于30 毫米时,存
在较大误差,模型很不稳定,误差图如下:
图十三
特别指出的是,当k在某个数附近时出现了极不稳定的情况,误差也达到了
极大值。
同时,模拟结果也表明,当k ³140时,即相机与模型的距离大约大于200
毫米时,误差始终保持在1 个像素以内,而且随着k 的增加误差D 在不断减小。
数学中国论文共享
www.madio.cn
19
图十四
这个结论表明,在相机距离标靶在200 毫米(即20厘米)以上时有着非常
高的精度以及稳定性,而我们实际拍摄时的物距很少小于20厘米,而且本文中
的物距大约是500 毫米左右,由图可以知道当距离大于500 毫米时,误差仅仅在
0.1 个像素之内,所以可以认为我们的模型是比较理想的。
2.像平面与圆平面之间的偏角对结果的影响。
2.1 a 的变化对结果的影响
从直观上来说,a 的变化对应相机的“上仰”或“下翻”,
2
p
a = ± 对应的相
机状态为镜头水平朝上或水平向下,所以当a 变化时是会带来误差的。
令g = 0, b = -0.7,T = (-350,0, -350), [ , ]
2 2
p p
aÎ - ,在计算机处理时,我
们对[ , ]
2 2
p p
- 这个区间进行等间隔采样已得到a 值,对于每一个a 值计算出我们
模型的解与理论解的误差。
图十五
数学中国论文共享
www.madio.cn
20
由图可以知道当a 在[ , ]
2 2
p p
- 变化时,误差始终在2 个像素以内,而且在
a = 0时取得最小误差。特别是当-0.5 £ a £ 0.5时误差小于一个像素,对我们模
型的结果几乎没有产生影响。
2.2 b 的变化对结果的影响
从直观上来说,在a = 0的前提下,像平面和圆所在平面都在铅垂面,b 的
变化直接对应着像平面和圆所在平面的夹角的变化,所以也会带来误差。我们采
用同样的方法用于检验b 的变化对结果的影响。这里我们令g = a = 0,
T = (-350,0, -350), [ , ]
2 2
p p
bÎ - ,计算出我们模型的解与理论解的误差
图十六
从图中我们发现,b 的变化对结果的产生的误差稍大于a ,但是最大的误差
只有2 个像素左右,但是大约在-1£ b £ 0.5时误差仍然小于1 个像素,可以说
对我们的模型的结果几乎没有影响。
在我们模拟出来的结果中,我们发现除了在b = 0外还有一处的误差达到局
部最小误差且接近于零,这是在我们的模型里没有得出的结论,我们通过大量实
验模拟发现这个结论总是成立,而且在通过对数据的反推得到下面的猜想。
猜想:在本题所建的坐标系中,若世界坐标系的Z 轴过光心C O ,则用我们
的模型可以得到精确的圆心像在像坐标系中的坐标。
3 总结
总体而言,在合理的相机与标靶的相对位置的前提下,我们的模型是很理想
的:边缘提取以及线性回归的误差可以忽略不计;像平面与圆平面之间的偏角对
结果也几乎不会造成影响,相比而言相机与标靶在距离内在小于200 毫米会造成
数学中国论文共享
www.madio.cn
21
很大的误差,而且在小于20毫米时极不稳定,但是在大于200毫米时误差始终
不超过一个像素,而且随着距离的增大误差趋向于零。所以对于合理的相机与标
靶的相对位置,我们的模型始终是有效而精确的。
问题四
用两部相机进行拍摄定位时,就会存在两个相机坐标系,我们记为
C1 : XC1YC1ZC1、C2 : XC2YC2ZC2、,世界坐标系与两者的关系分别为1 1 2 2 R、T和R、T
,
对任意一点,其在世界坐标系与两个相机坐标系下的坐标为W C1 C2 P 、P 、P ,则由
上文可得:
1 1 1
2 2 2
C W
C W
P RP T
P RP T
= +
= +
消去W P 得:
1 1
C1 1 2 C2 1 2 2 1 P = R R - P - R R - T +T
令 1 1
1 2 1 2 2 1 R = R R -、T = R R - T +T ,则可得到两相机坐标系的转换关系:
1 2 .................(12) C C P = RP +T
对两相机分别标定,则可得到相应的1 1 2 2 R、T、R、T ,从而可得两相机坐标
系的转换关系,即得到这两部相机的相对位置。
因为有:
0
1
1
W
W
C
W
X
x
Y
Z y M
Z
é ù
é ù ê ú
ê ú = ê ú ê ú ê ú
êë úû ê ú
ë û
所以对于世界坐标系里一点P在两个光心坐标系的两个像点1 P 2 、P ,有:
1 1 1 1
1 11 12 13 14
1 1 1 1
1 1 21 22 23 24
1 1 1 1
31 32 33 34 1
1
W
W
C
W
X
x a a a a
Y
Z y a a a a
Z
a a a a
é ù
é ù é ù ê ú
ê ú = ê ú ê ú ê ú ê ú ê ú êë úû ëê ûú ê ú ë û
2 2 2 2
2 11 12 13 14
2 2 2 2
2 2 21 22 23 24
2 2 2 2
31 32 33 34 1
1
W
W
C
W
X
x a a a a
Y
Z y a a a a
Z
a a a a
é ù
é ù é ù ê ú
ê ú = ê ú ê ú ê ú ê ú ê ú êë úû êë úû ê ú ë û
式中, 1 1 2 2 (x , y ,1)、(x , y ,1) 分别为1 P 2 、P
在两个光心坐标系的奇次坐标,
数学中国论文共享
www.madio.cn
22
( , , ,1) W W W X Y Z 为P在世界坐标系里的奇次坐标。展开上式可得:
1 1 1 1 1 1 1 1
1 31 11 1 32 12 1 33 13 1 34 14
1 1 1 1 1 1 1 1
1 31 21 1 32 22 1 33 23 1 34 24
( ) ( ) ( )
..........(13)
( ) ( ) ( )
W W W
W W W
xa a X xa a Y xa a Z xa a
ya a X ya a Y ya a Z ya a
ì - + - + - = - ïí
îï - + - + - = -
2 2 2 2 2 2 2 2
2 31 11 2 32 12 2 33 13 2 34 14
2 2 2 2 2 2 2 2
2 31 21 2 32 22 2 33 23 2 34 24
( ) ( ) ( )
..........(14)
( ) ( ) ( )
W W W
W W W
xa a X xa a Y xa a Z xa a
ya a X ya a Y ya a Z ya a
ì - + - + - = - ïí
îï - + - + - = -
方程(a)表示过1 P点在世界坐标系下的直线方程,方程(b)表示过2 P 点在世界坐标
系下的直线方程,那么(a)与(b)联立就可以得到两直线的交点,即得到原来在世
界坐标系下的物点坐标,这样我们就通过双目定位的方法,从两相机中像在光心
坐标系中的参数推导出了物在世界坐标中的参数,达到了定位的目的.如图:
图十七
现在我们再给出另外一种可行的模型:
根据假设,我们认为相机的光心和像屏中心的连线垂直于象平面。
k
A'
A
图十八
数学中国论文共享
www.madio.cn
23
问题一、二中我们已经求得了靶标上圆在像屏上面的坐标,以摄像机的光心为原
点,建立空间中的直角坐标系。x - y平面平行于像平面,Z 轴即为相机的光轴。
可以标出象平面上所有点的坐标,这里设为1 1 1 2 2 2 3 3 3 (a ,b , c )、(a ,b , c
)、(a ,b , c ) ,这三个
点对应在靶标平面上的三个点,这样点A的到原点地向量为1 1 1 (a ,b , c ),由于A,
O,A'在同一条直线上,那么A’坐标可以表示为1 1 1 1 k (a ,b ,c ) ,同理B'C '点的坐标
为2 2 2 2 k (a ,b , c ), 3 3 3 3 k (a ,b ,c
),由于ABC三个点之间的相对位置是固定的,也
即AB、AC、BC的长度是已知的,设为1 2 3 L、L、L ,可以得。
2 2 2 2
1 1 2 2 1 1 2 2 1 1 2 2 |AB| =(k a -k a ) +(k b -k b ) +(k c
-k c ) ;
2 2 2 2
3 3 2 2 3 3 2 2 3 3 2 2 |BC| =(k a -k a ) +(k b -k b ) +(k c
-k c ) ;
2 2 2 2
3 3 1 1 3 3 1 1 3 3 1 1 |AC| =(k a -k a ) +(k b -k b ) +(k c
-k c )
即得到三个方程,未知数为1 2 3 k、k 、k ,下面就可以解得三个方程。
该方程展开后是这种形式,
2 2
1 1 1 2 12 1
2 2
2 2 2 3 2 3 2
2 2
3 3 3 1 31 3
k p kk qk r
k p kk qk r
k p kk qk r
+ + =
+ + =
+ + =
该方程直接求解比较困难,附录中提供了一种求解方案。
这样靶标上的点关于我们建立的坐标系的绝对位置可以求得。这三个点为
A, B,C。
同样的道理我们可以得到靶标上三个点关于第二个相机为原点坐标系的坐标。如
图:
图十九
在求得A 点关于两个坐标系的坐标后,也就是知道了向量1 2 AO AO
uuuur uuuur
、,那么就可
数学中国论文共享
www.madio.cn
24
以知道O2O1 = AO1 - AO2
uuuur uuuur
即两个相机的相对位置。
五、模型的思考
1.离心率问题
观察题中所给的像图,可以很明显的看出不同椭圆的离心率是不相同的,那
么是什么原因导致同一世界坐标系下不同位置的圆经转换后而形成离心率不同
的椭圆的呢?
前文已经给出坐标系转换的变换矩阵M ,如下:
0 0 0
1 , 0 0 0
0 ,1
1 0 0 1 0
1 1
W W
W W
T
C W W
X X
x f
R TY Y
y f M
Z Z Z
é ù é ù
é ù é ù ê ú ê ú ê ú ê ú é ù ê ú ê ú ê ú = ê ú ê ú = ë û ê ú ê
ú êë úû êë úû ê ú ê ú
ë û ë û
观察坐标变换矩阵,对世界坐标系中的不同点,R、T、f 都有相同的值,唯
一的不同就是C Z ,也就说明最终导致离心率不同的因素就是C Z ,即点在相机坐
标系里的C Z 方向上的坐标。设想如题中所给定的标靶平面与光心坐标系的
C C X -Y 平面平行,那么很直观的可以到处最后所成的像都是相同的,即有相同
的离心率。
2.我们分别给出了两种相机定标的方法,并且可以互为印证。变换矩阵模型具有
很强的适应性,对一般的双目定标问题给出了求解的方法。公切线模型主要是在
平面上讨论问题,比较直观,容易理解;
3.本文是建立在小孔成像的基础之上,而实际相机成像是透镜成像,远离图像中
心处,镜头畸变会比较大,从而会给小孔成像模型带来一定误差;
4.矩阵模型中对旋转矩阵R 的求解不能非常清晰明确的给出,需要进一步完善。
数学中国论文共享
www.madio.cn
25
参考文献
[1] 马颂德张正友,计算机视觉,北京:科学出版社,1998
[2] 姜启源谢金星 叶俊,数学模型,北京:高等教育出版社,2003
[3] 李庆扬等,非线性方程组的数值解法,北京:科学出版社,1997
[4] 陈亚浙吴兰成,二阶椭圆型方程与椭圆型方程组,北京:科学出版社,1991
数学中国论文共享
www.madio.cn
26
附录
1.坐标系旋转变换
对于二维的坐标旋转变换,如由xoy绕o 点逆时针旋转q 到新的坐标系
x 'oy ',则有以下变换关系:
'cos sin
'sin 'cos
x x y
y x y
q q
q q
= - ìí
î = +
即:
' cos sin
' sin cos
x x y
y x y
q q
q q
= + ìí =- + î
此时有:
cos sin
sin cos
R
q q
q q
é ù
= ê ú ë- û
且R是正交的单位矩阵。
推广到三维情况,三维坐标系里的旋转我们可以等效为坐标系分别以
x、y、z 轴为转轴旋转,统一规定旋转方向为正面坐标系的逆时针旋转,绕
x、y、z轴旋转的角度分别为a、b、g 。
初始坐标为(x, y, z)经绕x轴旋转a 后得:
1
1
1
cos sin
sin cos
x x
y y z
z y z
a a
a a
= ìï
= + íï =- + î
再绕y轴旋转b 可得:
2 1 1
2 1
2 1 1
sin cos
cos sin
x z x
y y
z z x
b b
b b
= - + ìï
= íï
î = +
再经z轴旋转g 可得:
3 2 2
3 2 2
3 2
cos sin
sin cos
x x y
y x y
z z
g g
g g
= + ìï =- + íï
î =
那么此三维坐标系经旋转最后可得到坐标系XYZ 与原坐标系的坐标转换关
数学中国论文共享
www.madio.cn
27
系为
3
3
3
cos cos (sin sin cos cos sin ) ( cos sin cos sin sin )
( cos sin ) ( sin sin sin cos cos ) (cos sin sin sin cos
)
sin ( sin cos ) cos cos
X x x y z
Y y x y z
Z z x y z
b g a b g a g a b g a g
b g a b g a g a b g a g
b a b a b
= = + + + - + ìï
= = - + - + + + íï
î = = + - +
所以三维情况下的坐标旋转矩阵为:
cos cos sin sin cos cos sin cos sin cos sin sin
cos sin sin sin sin cos cos cos sin sin sin cos
sin sin cos cos cos
R
b g a b g a g a b g a g
b g a b g a g a b g a g
b a b a b
é + - + ù
= ê- - + + ú ê ú
êë - úû
非线性方程组的求解
有映象
F:D Í Rn ®Rn。求解非线性方程组
F(x)=0 ,假定X*ÎD是方程的一个精确解, Xk是X*的一个近似,通过Xk可
以定义仿射映像:
n n
k L:R
®
R
为
k k
k k L (x)=A (x-x )+F(x ) ,其中
n
k AÎ
L(R
)
为
非
奇
异
n
阶
矩
阵
,
显
然
k
k L (x)=F(x ),若用线性方程组
k k
k k L (x)=A (x-x )+F(x )=0的解x=xk+1作为方程的新
近似,即
k 1 k 1 ( k ), 0,1......
k x + = x - A - F x k = 即为非线性方程的线性化迭代法。
对所有的k 都取( ) n
k Aº
A
Î
L
R
非
奇
异
,
于
是 xk +1 = xk - A-1F(xk ), k = 0,1...... ,
成为n维平行弦方法。
应用Newton法解非线性方程组,还可以得到更好的收敛速度和自校正的特
点,但由于时间限制,我们并没有深入分析下去。
数学中国论文共享
www.madio.cn