建模第三次真题训练195队

1、 2016高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。 我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置

2、报名号的话）： 所属学校（请填写完整的全名）： 青岛科技大学 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月日赛区评阅编号（由赛区组委会评阅前进行编号）： 2009高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）： 数码相机定位 摘要 本文通过建立相机线型投影模型实现了数码相机双目标定中物象间坐标的转换，建立针孔线性模型进行相机标定，确立了从物平面到像之间之

3、间一一对应的关系。 首先，本文建立了四个坐标系：世界坐标系、相机坐标系、图像坐标系、像素坐标系。 针对问题一，依据光学的小孔成像原理，考虑到相机与物平面不在同一个平面，以及一些畸变因子的干扰，必须要通过空间内坐标变换在像平面和物平面之间找到一一对应的转换关系，这个转换过程是由世界坐标系相机坐标系理想图像坐标实际图像坐标计算机图像坐标。建立针孔线性模型，运用直接线性变换的方法，对变换过程中的一系列公式求解后，最终总结出一个算法计算靶标上圆的圆心在像平面的像坐标。 针对问题二，实质上是对问题一所建立的模型的求解过程。由于靶标上的圆映射到像平面后会变形，所以首先取利益中确定在像平面内根据变形后圆周上

4、的点找到圆心的方法，然后依据问题一中所建立的模型，根据五个圆的圆心在物平面和像平面的坐标位置，利用直接线性转换法结合最小二乘法计算旋转正交矩阵R和平移矩阵T，从而确立相机与物平面的相对位置。针对问题三，实质上是对问题一所建模型进行的对精度和稳定性模型检验. 为了检验本模型的稳定性，本文采用计算机模拟的方法来破坏像的轮廓并进行求解，比较两个结果就可以判断模型稳定性的好坏。为了简便起见，本文只考虑轮廓缩小的情况，不考虑轮廓放大的情况。对于轮廓缩小的情况，随机减少轮廓上的点来模拟图像轮廓的变化。在轮廓上减少点，利用模拟退火算法求解外公切线交点坐标。 针对问题四，首先每个相机旋转变换矩阵 R 和平移向

5、量 T， 可以得到两相机的变换关系：即相对位置关系，并理论推导了从两相机中像在光心坐标系中的参数得到物在世界坐标中的参数，实现双目定位。 另外在相机的光心和像屏中心的连线垂直于象平面基础上， 我们还给出另外一种合理模型，通过矢量的方法求出物相对于光心坐标系的精确位置， 从而可以得到两相机的相对位置关键词 小孔成像；相机标定；针孔线性模型；直接线性变换法； 1、 问题重述数码相机定位在交通监管（电子警察）等方面有广泛的应用。所谓数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的定位方法是双目定位，即用两部相机来定位。对物体上一个特征点，用两部固定于不同位置的相机摄得物

6、体的像，分别获得该点在两部相机像平面上的坐标。只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。于是对双目定位，精确地确定两部相机的相对位置就是关键，这一过程称为系统标定。标定的一种做法是：在一块平板上画若干个点， 同时用这两部相机照相，分别得到这些点在它们像平面上的像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置。然而，无论在物平面或像平面上我们都无法直接得到没有几何尺寸的“点”。实际的做法是在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点了。而它们的像一般会变形，所以必须从靶标上的这些圆的像中把圆心的像精确

7、地找到，标定就可实现。 现有一个靶标如下，取1个边长为100mm的正方形，分别以四个顶点（对应为A、C、D、E）为圆心，12mm为半径作圆。以AC边上距离A点30mm处的B为圆心，12mm为半径作圆，用一位置固定的数码相机摄得其像。要求我们处理一下几个问题：（1） 建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标, 这里坐标系原点取在该相机的光学中心，x-y平面平行于像平面；（2） 对给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标, 该相机的像距（即光学中心到像平面的距离）是1577个像素单位(1毫米约为3.78个像素单位)，相机分辨率为1024×768；（3）

8、设计一种方法检验模型，并对方法的精度和稳定性进行讨论；（4） 建立用此靶标给出两部固定相机相对位置的数学模型和方法。2、 模型假设与约定（1）假设相机的焦距是固定的，不存在相机自动调焦的情况，固定两相机像素单位长度相同（2）假设所给靶标上圆的相对位置以及投影后得到的像图足够精确；（3）假设相机两个方向上的焦距很小，像距约等于焦距；（4）假设靶标像的中心恰好在光轴上；（5）假设两架相机的相对位置固定且两架相机在空间中的位置相对固定；（6）假设不考虑光线经过透镜时的畸变，假设光线是以直线传播的方式到达光学中心；（7）假设数码相机中图像平面与光轴垂直；（8）假设忽略边缘图像颜色不确定带来的误差；三、

9、符号说明及名词解释序号符号符号说明4、 问题一 模型建立、求解及方案设计4.1问题分析 要求建立确定靶标上圆的圆心在像平面的坐标的模型。由于存在畸变，所以像平面的图形并不是规则的形状，单纯从像平面的图形出发显然不能找到这个圆心在像平面的坐标。题目中给出了靶标的尺寸等信息，但是靶标像上的点却很难与靶标上的点建立起对应关系。因此我们需要将模型简化，简化成一个小孔成像的模型，建立物平面到像之间之间一一对应的关系。由于靶标上的圆映射到像平面后会变形，所以首先取利益中确定在像平面内根据变形后圆周上的点找到圆心的方法，根据五个圆的圆心在物平面和像平面的坐标位置，利用直接线性变换法结合最小二乘法计算旋转正交

10、矩阵R和平移矩阵T，从而确立相机与物平面的相对位置。4.2模型建立与求解（一）模型准备：1：有关小孔成像原理的介绍：通过查阅有关资料，我们知道对于小孔成像模型有以下几个性质（性质的证明略）： 图一 图二（1） ：直线经小孔成像后仍然为直线（见图一）。（2） ：线段中点经小孔成像后所得的点不一定还是像线段的中点（见图一）。（3） ：若两直线相交，则两直线交点的像仍然是两直线像的交点（见图二）。（4） ：圆经小孔成像为椭圆（见图三）。（5）：圆的某一条切线的切点的像，仍然是椭圆的切线，切点的像也是椭圆的切点（见图三）。 图三2：有关坐标系的定义和介绍（1）世界坐标系：世界坐标系是用于描述物体空间位

11、置的一个基准绝对坐标系，表示为。相机所拍摄的内容还取决于相机在世界坐标系中的位姿，这个位姿参数就是相机的外参数。在本题中为方便求解，设定世界坐标系的原点为正方形靶标的中心，其Z轴垂直于靶标平面，X轴、Y同正方形相邻边分别平行。（2）相机坐标系 相机坐标系是以相机的光学中心为原点，以光轴为z轴的坐标系，对应于题中所定义的笛卡尔坐标系，具体表示为：。（3）图像坐标系 成像平面坐标系：以透镜光轴和成像平面的交点为原点，轴和轴分别平行于相机坐标系的轴和轴，有效焦距就是光学中心到成像平面的距离。点在理想小孔透视变换模型下的图像坐标为，考虑到透镜存在畸变，它的实际成像坐标为；在理想模型中，不考虑相机镜头畸

12、变，和重合。（4）像素坐标系相机采集的数字图像在计算机内可以存储为数组，数组中的每一个元素称为像素，像素的值即为图像点的亮度或灰暗度，若为彩色图像，则图像的像素亮度将有红、绿、蓝三种颜色的亮度表示。在图像上定义直角坐标系，每一像素的坐标，分别是该像素在数组中的列数和行数。所以，是以像素为单位的图像坐标系坐标。3：关于模型数据的预处理根据问题二的要求，我们知道相机的分辨率为1024×768，为了符合题目的要求，我们通过相关的修图软件修改图片的大小为1024×768像素，去除其它颜色的干扰，转化成单色图，即图像中只有黑白两种颜色，然后将该图以1024行、768列的矩阵存入计算机

13、，利用matlab中的edge函数求出各圆的边缘点，并将边缘点的位置坐标保存保存下来，供求解模型使用（求解程序见附录二）。边缘点的图像见图四：图四 边缘点的分布（2） ：模型的建立在拍摄照片的过程中，外界景物经历了从世界坐标到相机坐标再到平面像坐标的变换，为了建立模型并求解，首先要了解世界坐标系、相机坐标系、平面像坐标系之间的映射流程及其变换。从上面的叙述中我们已经知道世界坐标到平面坐标的映射流程为：世界坐标相机坐标平面像坐标。下面说明各个坐标系之间的变换关系（图五）。图五 各个坐标系之间的关系1.世界坐标系到相机坐标系的变换世界坐标系中的坐标转换为相机坐标系中的坐标包含平移和旋转两种变换，用

14、矩阵型时表示为： （4.1）其中，R为旋转正交矩阵，T为平移矩阵。R决定了相机相对于世界坐标系的方向；T决定了相机相对于世界坐标系的位置。 （4.2）2.相机坐标系和成像平面坐标系的针孔透视变化在理想情况下，根据针孔相机的几何结构，用透视投影变换将相机坐标系下的三维坐标变换为图像坐标系下的而为坐标，有几何关系（式4.3）： (4.3)式中，为相机焦距。3.成像平面坐标系到计算机坐标系之间的转换由相机获得的被测物体的图像，最终有图像采集卡变成数字图像并输入计算机。从相机图像坐标到计算机图像坐标之间的转换为： （4.4）式中，为相机坐标实际成像坐标，和分别为水平方向和垂直方向的像素长，为计算机图像

15、的光学中心坐标。在理想情况下，和重合，得到： （4.5）综合以上的变换过程，相机的理想变换模型为： （4.6）上式写成矩阵形式即为： （4.7）（二）相机标定方法直接线性变换方法（DLT）将像点和物点的成像几何关系在齐次坐标下写成透视投影矩阵的形式： （4.8）其中，为第i个点在世界坐标系中的其次坐标，为第i个点在图像坐标系中的图像坐标，为投影矩阵R的第i行第j列元素，s为未知尺度因子。将上述矩阵方程展开： （4.9）整理式（4.7）（4.9）得： （4.10）再次整理得： （4.11）将上式写成矩阵形式： （4.12） 将上式简写成那么DLT参数的最小二乘解为： （4.12）知道不同坐标系之

16、间的变换后，我们就能将世界坐标转化为相机像平面上的坐标。根据已有数据，以E的圆心为原点，建立空间直角坐标系如图六所示：zyx5.2.1%wdx.mbound_zb;dqp=0.3;wdx_zb=zb_bound;total_result=;for cs=1:10 zb_bound=wdx_zb; for i=1:5 for j=length(zb_boundi(:,1):-1:1 ranum=rand(1); if ranum<dqp zb_boundi(j,:)=; end end end qj; total_result=total_result,result;end %qj.mqx

17、;res_l=fsolve(fu,1500 1500 1500 1500);load -ascii c.txt;ca=c;oa=res_l(1)*ca(1,:);ob=res_l(2)*ca(2,:);oc=res_l(3)*ca(3,:);od=res_l(4)*ca(4,:);ab=ob-oa;ac=oc-oa;cd=od-oc;bd=od-ob;a3a=(12/124)*ac+(12/124)*ab;a3b=(42/124)*ac+(12/124)*ab;a3c=(112/124)*ac+(12/124)*ab;a3d=(112/124)*ac+(112/124)*ab;a3e=(12/

18、124)*ac+(112/124)*ab;oA=oa+a3a;oB=oa+a3b;oC=oa+a3c;oD=oa+a3d;oE=oa+a3e;A=xj2wl(oA);B=xj2wl(oB);C=xj2wl(oC);D=xj2wl(oD);E=xj2wl(oE);A;B;C;D;E;result=round(wl2xs(A);round(wl2xs(B);round(wl2xs(C);round(wl2xs(D);round(wl2xs(E)%qx.m%求切线主程序cd 'qx'sa;%best=143 285 539 240 177 685 537 613;zb0=best(2

19、) best(1);best(4) best(3);best(6) best(5);best(8) best(7);zb_x=;for i=1:length(zb0(:,1)zb_x=zb_x;wl2xj(xs2wl(zb0(i,:);endsave c.txt zb_x -asciiclear zb0;%fu.mfunction f=fu(l)load -ascii c.txt;oa=l(1)*c(1,:);ob=l(2)*c(2,:);oc=l(3)*c(3,:);od=l(4)*c(4,:);ab=ob-oa;ac=oc-oa;cd=od-oc;bd=od-ob;f(1)=dot(ab,

20、ac);f(2)=dot(ac,cd);f(3)=dot(cd,bd);f(4)=dot(ab,bd);f(5)=norm(ab)-468.72;f(6)=norm(ac)-468.72;f(7)=norm(cd)-468.72;f(8)=norm(bd)-468.72;return;%sa.m%模拟退火算法r=143,285,539,240,177,685,537,613;%初始解zb_y=zb_bound1;zb_bound2;zb_bound3;zb_bound4;zb_bound5;Tm=5;T=Tm;%初始温度Tn=1;%最低温度E=;%能量dE=;%能量差k=0;%迭代次数l=le

21、ngth(r);L0=3;L1=1;%随机扰动范围随温度变化，最大和最小范围best=r;zx=mbFunction(r);while T>=Tnk=k+1;h=0;while h=0hh=0;while hh=0r1=r+randomDisturb(l,Tm,T,Tn,L0,L1);if isMatch(r1,zb_y)hh=1;endend%扰动产生一个新可行解E=mbFunction(r1);dE=E-mbFunction(r);if dE<=0r=r1;h=1;elsep=exp(-dE/T);if p>=rand(1)r=r1;h=1;endh=1;endif mo

22、d(k,10)=0T=0.9*T;T;endif E<zxzx=E;best=r;endendendbestmbFunction(best)clear E;clear T;clear Tn;clear dE;clear h;clear hh;clear k;clear l;clear p;clear r1;clear L0;clearL1;clear Tm;clear r;clear zb_y;clear zx;cd '.' 图六 对靶标建立空间直角坐标系对于图六所示的坐标系，我们根据题目中给的已知条件得到在所建的坐标系上，物平面的五个圆的圆心坐标如表一所示：表一 物平面

23、的空间坐标ABCDE坐标（0,0,100）（0,30,100）（0,100,100）（0,100,0）（0,0,0）在得知物平面上圆心坐标的情况下，我们就可以通过方程4.1求出像平面上的圆心坐标，但是R和T这两个参数未知，我们可以根据两图中圆心的位置关系求出R,T的值，从而确定相机与物平面的关系。通过图四可以明显的看到，像平面上的圆发生了畸变，说明相机和物平面不平行，而要确定靶标上的圆的圆心在像平面上的坐标，只能通过畸变后的图像利用有关成像的性质寻该圆的圆心。我们利用小孔成像的性质3和性质4确定了采用最值连线求交点的方法计算靶标在像平面上的圆心，具体的做法如下：为了适应matlab的计算模式，

24、我们以左上角为圆心，横向为x轴，纵向为y轴建立坐标系（见图七），在图中找出x值最大值和最小值点A,C，y值得最大值最小值点B,D，连接AC,BD，它们的交点O即为靶标在像平面上的圆心（可以通过几何方法加以证明，证明略）。图七 椭圆圆心的确定下面通过解方程的方法确定交点O：设，解方程4.13可以得到O点的坐标。 （4.13）为了避免在图中搜索最值的时候出现多个圆的重合现象，我们采用对每一个圆的边缘值单独搜索的方法进行计算，最终求出各个圆的最值点的坐标见表二：各个圆的边缘点集的坐标的计算见附录二；表二 各个圆边缘点的最值坐标圆12345x最大坐标（227，315）（233,409）（248,635

25、）（533,280）（531,570）x最小坐标（150,320）（158,425）（177,639）（468,279）（471,582）y最大坐标（180，362）（184，460）（205,674）（489,321）（496,616）y最小坐标（182,284）（197,385）（210,606）（499,248）（499,551）通过求出的边缘点的最值坐标结合方程4.13,我们能够求出每个椭圆的中心坐标，也即是在像平面上的圆心坐标（见表三）。表三 像平面上的圆心坐标圆圆心坐标1（181.13，317.98）2（191.30,417.90）3（207.70,637.27）4（494.70,

26、279.41）5（497.82,576.64）现在我们已经求出了实际靶标上的圆心在像平面上的坐标，下面我们就能计算R,T参数，从而确定方程4.1，然后通过靶标上的坐标计算出不考虑畸变因子的情况下的理想坐标。通过已经建立的物平面和像平面坐标之间的关系模型，借助matlab软件编程求解可以得到在理论上像平面上的圆心坐标（计算的程序见附录二），像平面的圆心坐标的理论值见表二： 五、模型的检验5.2 稳定性检验 由于种种原因，像的轮廓上的点很多情况下是难以确定的，而这些难以确定的点常常会影响模型求解的结果。一个好的模型不仅要精确，而且在轮廓出现变化时，模型的结果不能变化太大，即要求模型的稳定性要好。

27、为了检验本模型的稳定性，我们采用计算机模拟的方法来破坏像的轮廓并进行求解，比较两个结果就可以判断模型稳定性的好坏。为了简便起见，我们只考虑轮廓缩小的情况，不考虑轮廓放大的情况。对于轮廓缩小的情况，我们可以随机减少轮廓上的点来模拟图像轮廓的变化。在轮廓上减少点，利用模拟退火算法求解外公切线交点坐标时，其最优解目标函数值只会比减少钱更小，因此不妨令模拟退火算法的初始点为以上我们模型确定的点 。 令 p 为轮廓上点的保留概率，则轮廓上的点被删除的概率为 1-p，我们在 3 个不同的 p 值下利用以上模型计算圆心的坐标。 每个 p 值做了 10 次模拟 （模拟程序见附录5.2.1 ），模拟的结果如下：

28、P总偏移距离平均偏移距离单个点最大偏移距离0.910.0210.850.110.713.4140.26821.414其中每个点的偏移距离即新的圆心与前面模型求解出的圆心之间的距离； 总偏移距离即 10 在一个 p 下 10 次模拟每个点偏移距离的综合；每个 p 下做 10 次模拟，每次模拟有 5 个坐标，因此每个 p 下有 50 个点，平均偏移距离既为总偏移距离除以 50。表中各数据单位为像素。 从以上数据可以看出，在轮廓发生改变后，即使在轮廓损失了近 30% 的信息量时，圆心的平均偏移距离也只有 0.2682，不到一个像素，因此该求解圆心的模型具有很好的稳定性。 六、问题四 模型建立、求解及

29、方案设计 6.1问题分析： 问题要求建立用此靶标给出两部固定相机相对位置的数学模型和方法。用两部相机进行拍摄定位时，就会存在两个相机坐标系，所以我们考虑建立两个世界坐标系可得到两相机坐标系的转换关系，就可以得到这两部相机的相对位置。每个相机旋转变换矩阵 R 和平移向量 T， 可以得到两相机的变换关系：即相对位置关系，并理论推导了从两相机中像在光心坐标系中的参数得到物在世界坐标中的参数，我们考虑得到原来在世界坐标系下的物点坐标，这样我们就可以通过双目定位的方法，从两相机中像在光心坐标系中的参数推导出了物在世界坐标中的参数，达到了定位的目的6.2 模型建立与求解： 用两部相机进行拍摄定位时，就会存

30、在两个相机坐标系，我们记为、，世界坐标系与两者的关系分别为和,对任意一点，其在世界坐标系与两个相机坐标系下的坐标为，则由上文可得： 消去得： 令 ,则可得到两相机坐标系的转换关系: 对两相机分别标定，则可得到相应的从而可得两相机坐标系的转换关系，即得到这两部相机的相对位置。因为有: 所以对于世界坐标系里一点P在两个光心坐标系的两个像点,有: 式 中 ,、 分 别 为 在 两 个 光 心 坐 标 系 的 奇 次 坐 标 ，为P在世界坐标系里的奇次坐标。展开上式可得：方程（a） 表示过点在世界坐标系下的直线方程， 方程(b)表示过点在世界坐标系下的直线方程，那么(a) 与(b)联立就可以得到两直线

31、的交点， 即得到原来在世界坐标系下的物点坐标，这样我们就通过双目定位的方法， 从两相机中像在光心坐标系中的参数推导出了物在世界坐标中的参数，达到了定位的目的如图： 现在我们再给出另外一种可行的模型：根据假设，我们认为相机的光心和像屏中心的连线垂直于象平面。问题一、二中我们已经求得了靶标上圆在像屏上面的坐标， 以摄像机的光心为原点，建立空间中的直角坐标系。平面平行于像平面，Z 轴即为相机的光轴。可以标出象平面上所有点的坐标，这里设为,这三个点对应在靶标平面上的三个点，这样点 A 的到原点地向量为由于 A，O，在同一条直线上，那么 坐标可以表示为,同理 点的坐标为由于 ABC 三个点之间的相对位置是固定的，也即 AB 、AC 、BC 的长度是已知的，设为，可以得。 即得到三个方程，未知数为，下面就可以解得三个方程。该方程展开后