重要性-绩效分析矩阵方法

重要性-绩效分析矩阵方法（完整版）实用资料(可以直接使用，可编辑完整版实用资料，欢迎下载）

重要性-绩效分析矩阵方法重要性-绩效分析矩阵方法（完整版）实用资料(可以直接使用，可编辑完整版实用资料，欢迎下载）第一步：基本数据收集表一：流程重要性—绩效信息收集表衡量指标分值（1-5分）ABCDEFGHI……备注流程重要性=∑(y1,y2,y3)/3y1流程的增值（=∑(A,…I)/9）分值从间接支持、直接支持、增值依次给1-5分y2流程的独特性（=∑(A,…I)/9）分值从一般的、独特的依次给1-5分y3流程类型（=∑(A,…I)/9）分值从支持、战术性的、战略性的依次给1-5分流程绩效性=∑(X1,X2…X5)/5X1质量指标（=∑(A,…I)/9）从产品的符合性、稳定性、过程性能指标来进行打分X2成本指标（=∑(A,…I)/9）从资源的消耗量来打分X3按时交付成果（=∑(A,…I)/9）主要从业务流程周期、关键点响应速度、成果交付时间来打分X4服务指标（=∑(A,…I)/9）主要通过顾客满意度调查来打分X5员工（=∑(A,…I)/9）主要从员工对流程的了解程度来进行打分第二步：重要性—绩效矩阵分析矩阵分析在信号处理中的应用班级：2021级专业硕士班在二十一世纪，人们普遍认为对人类社会发展最有影响的科学领域将是信息科学和生命科学。事实上，我们确实强烈感受到信息科学的飞速发展对我们日常生活方面的影响，信息科学的发展让世界变成了地球村。概括起来说，信息科学研究的是作为信息载体的信号的获取、存储、传输和处理。可见信号处理是信息科学的核心研究内容之一。信号处理在理论上涉及的范围极其广泛，并不断有新的分支出现。从所处理的信号的性质上来看，可分为确定性信号处理和随机信号处理。确定性信号处理研究的确定性信号的分析、线性滤波、重构，反卷积（线性失真补偿）等。除了大家比较熟悉的线性滤波器设计与实现理论、信号分析的各种快速变换算法等之外，还包括信号重构理论、多抽样率信号处理、小波分析等较新的学科分支。确定性信号处理时随机信号处理的重要理论基础之一。本文就是研究基于确定性信号处理的矩阵分析应用。平稳随机过程的功率谱估计方法可以分为经典谱估计方法和现代谱估计方法。经典谱估计的基本方法包括19世纪末由Schuster提出的周期图法和1949年Tukey根据维纳—辛钦定理提出的自相关法即BT法。现代谱估计方法主要是针对经典谱估计分辨率低和估计质量差提出的，因此也称作高分辨率谱估计。现代谱估计是基于随机信号的参数化模型表示的方法。在现代谱估计中，对未能得到的样本数据或未能估计出来的自相关函数，不是简单地当作零处理，而是与所得到的样本数据服从同一模型。根据谱表示定理，从功率谱等价的角度，规则过程可以用AR、MA或ARMA模型来描述。根据我们所学的信号与系统、数字信号处理等课程中系统的幅频特性与系统零、极点之间的关系，可以明确知道AR模型适合于具有尖峰但没有深谷的谱，MA模型适合于具有深谷但无尖峰的谱，ARMA模型是较通用的，对两种极端情况能够表示的模型。通过建立平稳随机过程的AR模型来获得功率谱估计的方法叫做AR谱估计。历史上，最早提出AR模型并用于时间序列建模的是Yule，他在1927年采用AR建模的方法研究太阳黑子的活动周期。Walker在1931年用最小二乘法建立了自回归模型参数与自相关函数关系的Yule-Walker方程。已知AR过程的前p+1个延迟的自相关函数，就可以解一组线性方程组来确定AR参数。Yule－Walker方程用矩阵形式表示：对于平稳过程，我们知道，自相关矩阵是非负定的，也就是说可能是正定的，也可能是半正定的。但对平稳非可预测过程，可以证明是正定的，即自相矩阵是满秩的。因此，由Yule－Walker方程可以唯一确定出AR系数。再由：可以确定白噪声的方差。一般将这两个方程用矩阵形式联立表示为：无论是p阶自相关矩阵，还是p+1阶自相关矩阵，它们在结构上都有特殊的性质，根据自相关函数的共轭对称性，它们是共轭对称的，即实平稳过程的自相关矩阵是对称的Toeplitz矩阵，而复平稳过程的自相关矩阵是共轭对称的Toeplitz矩阵。正是利用了自相关矩阵的共轭对称Toeplitz结构，Levinson提出了一种快速求解Yule—Walker方程的递推算法并由Durbin进行了改进，这种算法称为Levinson—Durbin算法。对于实际谱估计问题，过程的自相关函数是不可能精确知道的。一般只能得到过程的有限长度观测样本。最简单的AR模型参数估计方法是，在Yule—Walker方程中，用样本自相关函数估计代替自相关函数真值，并求解所得到的Yule－Walker近似方程，也就是说自相关法。具体地说，就是求解以下方程组：以及显然自相关矩阵估计是共轭对称的且具有Toeplitz结构，同时可以验证：因此，是正定的。所以，该方程组可用Levinson递推算结果求解，并且可以保证估计出来的AR模型的稳定性。另一种方法是协方差法。在自相关法中，预测误差功率的估计包含了一些不适当的项。协方差法与自相关法的唯一差别就是在估计预测误差功率时不取那些不适当的项。即前向线性预测误差功率估计为：其中：，，。由前向线性预测误差功率的极小化来估计AR模型的参数，即：为求这个二次函数极小化问题，求关于的导数，并令其等于零，可以得到：即：具体写出来就是：求解这个方程就可以得到。这个方程称为AR参数估计的协方差方程。对应的最小线性预测误差功率估计为：由于，因此，矩阵是共轭对称的，但不具备Toeplitz结构，方程的求解不再可以用Levinson算法。如果x(n)为非可预测过程，矩阵是正定的，因此，可以采用运算量较小的平方根法（Cholesky分解法）来求解。同时，利用由Toeplitz矩阵的乘积构成这一特点，还提出了其它更有效的求解方法。协方差法的缺点是理论上不能保证估计出的AR参数对应于一个稳定的AR模型。但其谱分辨率要优于自相关法。对于AR参数估计和AR功率谱估计的统计性能分析，是件很复杂、很困难的事。AR参数估计方法，AR参数估计值和AR功率谱估计都是样本数据的高度非线性函数，因此，要获得它们的偏差、方差等的一般解析表达式，几乎是不可能的。目前所能得到的仅是假设样本长度N很大情况下的一些近似结果，也即大样本渐近统计特性。在这种情况下，前面所介绍的各种AR参数估计方法都是十分接近的，也即它们有相同的渐近统计特性。设服从AR(p)模型，即其中的零点全部位于复平面的单位圆内，为零均值的独立同分布随机过程即纯白噪声过程，，则对自相关法估计得到的AR参数与，当样本长度时具有以下性质：(1)是的渐近无偏估计，具有一阶收敛(2)是的强一致估计，且(3)的分布收敛于Gauss分布，即如果服从高斯分布，则(4)如果服从高斯分布，则是的渐近有效估计，即矩阵与数值分析学院专业班级学号姓名电子信息与电气工程学部生物医学工程刘江涛1：考虑计算给定向量的范数；输入向量x=(x1,x2,,xn)T，输出x,x2,x∞，请编制一个通用程序，并用你编制的程序计算如下向量的范数：1⎫⎛11Tx=1,,,,⎪,y=(1,2,,n)n⎭⎝23对n=10,100,1000甚至更大的n计算其范数，你会发现什么结果？你能否修改你的程序使得计算结果相对精确呢？通用求范数程序：functionNORM(x)y1=sum(abs(x));y2=(sum(x.^2))^(1/2);y3=max(abs(x));fprintf('1-范数=%g；2-范数=%g；inf-范数=%g\n',y1,y2,y3);例题的运行程序：functionxianglaing(n)x=[];y=[];fori=1:nx(i)=1/i;y(i)=i;enddisp('x的范数：');NORM(x');disp('')disp('y的范数：');NORM(y');运行结果如下表：T根据上述的两个表的运行结果，我们可以得知无论n的值如何变化，对于x∞=1恒成立；y∞=n恒成立，其1-范数与2-范数随着n的增大而增大，但是其变化越来越小，这是因为计算在进行数值计算时有误差存在，对于表达式(1)当n很大时1却很n小，会出现“大数吃小数的现象”；修改方案：当n很大时我们避免用n做除数，因为当n非常大时1→0成立；所以在求解其范数时我们从小数开始相加，无穷个非常n小的数值相加也可能是个很大的数，从而可以避免两个数相加时出现“大数吃小数”的现象；2：考虑y=f(x)=ln(1+x)，其中定义f(0)=1，此时f(x)是连续函数，用此公x式计算当x∈[-10-15,10-15]时的函数值，画出图像。另一方面，考虑下面算法：d=1+x;ifd=1theny=1elsey=lnd/(d-1)endif用此算法计算x∈[-10-15,10-15]时的函数值，画出图像，比较一下发生了什么？程序：x=-10^(-15):10^(-20):10^(-15);if(x==0)f=1;elsef=log(1+x)/x;endfigure(1)plot(x,f);d=1+x;ifd==1y=1;elsey=log(d)/(d-1);endfigure(2)Plot(x,y);有图可知，直接用公式f(x)=ln(1+x)计算x∈[-10-15,10-15]的函数值时，除了在xx=0出的值为1，其他的值都是无限趋近于1；而利用算法二算出的结果全为1；出现这这情况的原因是x的取值非常接近于0，在用公式d=1+x求d得过程中出现了大数吃小数的情况，所以在用计算机计算时d=1恒成立，从而使y=1恒成立；3：首先编写一个利用秦九韶算法计算一个多项式在定点的函数值的通用程序，你的程序包括输入多项式的系数以及定点，输出函数值，利用你编写的程序计算f(x)=(x-2)9=x9-18x8+144x7-672x6+2021x5-4032x4+5376x3-4608x2+2304x-512在x=2邻域附近的值，画出p(x)在x∈[1.95,20.5]上地图像。秦九韶算法的通用程序：%A为多项式的以升幂排列的系数，x为初始值functionp=qinjiushao(A,x)a=A;[~,n]=size(a);n=n-1;S=[];S(n+1)=a(n+1);fork=n:-1:1S(k)=x.*S(k+1)+a(k);endp=S(1);利用上述程序计算p(x)在x=2邻域附近的值具体见下表：当x∈[1.95,20.5]时，p(x)的图像如下：画图程序如下：functionhuatu(A,x)[~,n]=size(x);fori=1:ny(i)=qinjiushao(A,x(i));endplot(x,y);程序运行如下：>>x=1.95:0.01:20.5;>>A=[-5122304-46085376-40322021-672144-181];>>huatu(A,x)11p(x)x4：编制计算机给定矩阵A的LU分解和PLU分解的通用程序，然后利用你编写的程序完成下面两个计算任务：考虑⎡1⎢-1⎢A=⎢⎢⎢-1⎢⎣-101⎤⎥⎥01⎥∈Rn⨯n⎥11⎥-11⎥⎦0-1-1自己取定x∈Rn，并计算b=Ax。然后用你编制的不选主元的Gauss消去法求解ˆ。对n从5到30估计计算解的精度。该方程组，记你计算出的解为x（2）对n从5到30计算出其逆矩阵。LU分解的通用程序：functionLU(A)[m,n]=size(A);L=zeros(m,n);U=zeros(m,n);forj=1:nU(1,j)=A(1,j);L(j,j)=1;endforj=2:mL(j,1)=A(j,1)/U(1,1);endfori=2:nforj=i:nsum=0;fork=1:i-1sum=sum+L(i,k)*U(k,j);endU(i,j)=A(i,j)-sum;endforj=i+1:nsum=0;fork=1:i-1sum=sum+L(j,k)*U(k,i);endL(j,i)=(A(j,i)-sum)/U(i,i);endenddisp('L=');disp(L);disp('U=');disp(U);PLU分解的通用程序：functionPLU(A)[~,n]=size(A);Ip=1:n;fork=1:n-1[~,r]=max(abs(A(k:n,k)));r=r+(k-1);ifr>kA([k,r],:)=A([r,k],:);Ip([k,r])=Ip([r,k]);endforp=k+1:nmu=A(p,k)/A(k,k);A(p,k)=mu;A(p,k+1:n)=A(p,k+1:n)-mu*A(k,k+1:n);endendp=eye(n,n);P=zeros(n,n);fori=1:nP(i,1:n)=p(Ip(i),1:n);endL=tril(A,-1)+eye(n);U=triu(A);disp('P=')disp(P);disp('L=')disp(L);disp('U=');disp(U);我选取b=[123n]Tn∈[5,30],n∈N*，则我们可以计算出其精确解2n-1-12n-2-1-1x=[-n-1-n-22222n-1T]2n-1n∈[5,30]n∈N*实现求解的程序如下：functionINV(n)A=ones(n);fori=2:nforj=1:i-1A(i,j)=-1*A(i,j);endforj=i:n-1A(i-1,j)=0;endendfori=1:nb(i)=i;X(i)=-(2^(n-i)-1)/(2^(n-i));endX(n)=(2^n-1)/(2^(n-1));[L,U]=LU(A);x1=inv(U)*inv(L)*b'disp(x1);我们利用上述的程序计算出的结果如下表：利用我们求出的精确解与用程序求出的近似解求其误差，再利用matlab编程实现时其误差为零。（2）对于题目中的A的求逆的通用程序：functionINV(n)A=ones(n);fori=2:nforj=1:i-1A(i,j)=-1*A(i,j);endforj=i:n-1A(i-1,j)=0;endendB=inv(A);disp('A的逆为：');disp(B);n=5时：-0.2500-0.1250-0.0625-0.0625⎤⎡0.5000⎢⎥00.5000-0.2500-0.1250-0.1250⎢⎥-1⎥A=⎢000.5000-0.2500-0.2500⎢⎥0000.5000-0.5000⎢⎥⎢⎥0.25000.12500.06250.0625⎣0.5000⎦n=6时：-0.2500-0.1250-0.0625-0.0313-0.0313⎤⎡0.5000⎢0⎥0.5000-0.2500-0.1250-0.0625-0.0625⎢⎥⎢0⎥00.5000-0.2500-0.1250-0.1250A-1=⎢⎥000.5000-0.2500-0.2500⎢0⎥⎢0⎥0000.5000-0.5000⎢⎥0.50000.25000.12500.06250.03130.0313⎢⎥⎣⎦n=7时：-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156⎤⎡0.5000⎢0⎥0.5000-0.2500-0.1250-0.0625-0.0313-0.0313⎢⎥⎢0⎥00.5000-0.2500-0.1250-0.0625-0.0625⎢⎥A-1=⎢0000.5000-0.2500-0.1250-0.1250⎥⎢00000.5000-0.2500-0.2500⎥⎢⎥000000.5000-0.5000⎢⎥⎢0.50000.25000.12500.06250.03130.01560.0156⎥⎣⎦n=8时：-0.2500-0.1250-0.0625-0.0313-0.0156-0.0078-0.0078⎤⎡0.5000⎢0⎥0.5000-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156⎢⎥⎢0⎥00.5000-0.2500-0.1250-0.0625-0.0313-0.0313⎢⎥0000.5000-0.2500-0.1250-0.0625-0.0625⎥A-1=⎢⎢00000.5000-0.2500-0.1250-0.1250⎥⎢⎥000000.5000-0.2500-0.2500⎢⎥⎢0000000.5000-0.5000⎥⎢⎥0.25000.12500.06250.03130.01560.00780.0078⎥⎢⎣0.5000⎦.......在此不再一一列举，都可以用上述程序算出；5：编制计算对称正定阵的Cholesky分解的通用程序，并利用你编制的程序计算Ax=b，其中A=(aij)∈Rn⨯n,aij=1，b可以有你自己取定，对n从10到20验i+j-1证程序的可靠性。Cholesky分解求L的通用程序：%LT代表L的转置function[L,LT]=Cholesky(A)[n,m]=size(A);L=zeros(n,n);forj=1:msum=0;fork=1:j-1sum=sum+(L(j,k))^2;endL(j,j)=(A(j,j)-sum)^(1/2);fori=j+1:nsum=0;fork=1:j-1sum=sum+L(i,k)*L(j,k);endL(i,j)=(A(i,j)-sum)/L(j,j);endendL=L;LT=L';求解Axb的方程解的程序如下：functioncholesky_qiu_jie(n,b)A=zeros(n);fori=1:nforj=1:nA(i,j)=1/(i+j-1);endend[L,LT]=Cholesky(A);Y=inv(LT)*inv(L)*b;disp(Y');>>n=10;>>b=[1000000000]';>>cholesky_qiu_jie(n,b)1.0e+06*Columns1through5Columns6through10-6.3052258749891819.607833253959356-8.7498890229612894.374900125373453-0.923581794866949>>n=11;>>b=[10000000000]';>>cholesky_qiu_jie(n,b)1.0e+07*Columns1through5Columns6through10Column110.385801129112326n=12;>>b=[111111111111]';>>cholesky_qiu_jie(n,b)1.0e+08*Columns1through5Columns6through10Columns11through12n=13;>>b=[1111111111111]';>>cholesky_qiu_jie(n,b)1.0e+10*Columns1through5Columns6through10-0.0955644689000000.338540141400000-0.795937993200000Columns11through13>>n=14;>>b=[11111111111110]';>>cholesky_qiu_jie(n,b)1.0e+15*Columns1through5Columns6through10Columns11through14-0.862636505012200>>n=15;>>b=[111111111111101]';>>cholesky_qiu_jie(n,b)1.0e+16*Columns1through5Columns6through10Columns11through15>>n=16;>>b=[1111111111111010]';>>cholesky_qiu_jie(n,b)1.0e+16*Columns1through5Columns6through10Columns11through15Column16>>n=17;>>b=[11111111111110101]';>>cholesky_qiu_jie(n,b)1.0e+17*Columns1through5Columns6through10-0.034029369522343Columns11through15Columns16through17>>n=18;>>b=[111111111111101012]';>>cholesky_qiu_jie(n,b)1.0e+16*Columns1through5Columns6through10Columns11through15Columns16through18-3.4602807804032691.887745959501761-0.436632060441746>>n=19;>>b=[1111111111111010120]';>>cholesky_qiu_jie(n,b)1.0e+17*Columns1through5Columns6through10Columns11through15Columns16through191.944997967101147-0.384572610629205n=20;b=[11111111111110101200]';>>cholesky_qiu_jie(n,b)1.0e+17*Columns1through5Columns6through10Columns11through15Columns16through203.295476378249488A-1利用cholesky分解求出的解的精确度高于直接，因为当n逐渐增大的过程中，越来越接近奇异矩阵，使得计算结果的误差增大，而使用cholesky分解可以避免这种现象的产生，是计算结果更加精确。6：（1）编制程序House(x)，其作用是对输入的向量x，输出单位向量u使得(I-2uuT)x=x2e1。编制Householder变换阵H=I-2uuT∈Rn⨯n乘以A∈Rn⨯m的程序HA，注意，你的程序并不显式的计算出H。考虑矩阵⎛1-1A=-2-0⎝232224⎫⎪2⎪eπ⎪⎪-37⎪75/2⎪⎭3用你编制的程序计算H使得HA的第一列为αe1，并将HA的结果显示出来。编制House(x)，其作用是对输入的向量x,输出单位向量U使得(I-2uuT)x=x2e1：House(x)通用程序：functionHouse(x)[n,~]=size(x);e1=zeros(n,1);e1(1)=1;w=x-norm(x)*e1;U=w/sqrt((w'*w));disp('U=');disp(U);编制Householder变换阵H=I-2uuT∈Rn⨯n乘以A∈Rn⨯n的程序HA，注意，你的程序并不显示的计算出H。Householder变换阵H通用程序：functionHouseholder(A)[n,~]=size(A);e1=zeros(n,1);e1(1)=1;x=A(:,1);w=x-norm(x)*e1;U=w/sqrt((w'*w));H=eye(n)-2*(U*U');HA=H*A;disp('HA=');disp(HA);考虑矩阵⎛1-1A=-2-0⎝24⎫⎪32⎪2eπ⎪⎪2-37⎪272⎪⎭3用你编制的程序计算H使得HA的第一列为αe1的形式，并将HA的结果显示：HA的结果显示的通用程序：functionHouseholder(A)[n,~]=size(A);e1=zeros(n,1);e1(1)=1;x=A(:,1);w=x-norm(x)*e1;U=w/sqrt((w'*w));H=eye(n)-2*(U*U');HA=H*A;disp('HA=');disp(HA);程序运行结果：>>A=[1234;-13sqrt(2)sqrt(3);-22exp(1)pi;-sqrt(10)2-37;0275/2];>>Householder(A)HA=4.0000-2.83111.4090-6.53780.00001.38960.8839-1.78050.0000-1.22081.6576-3.88360.0000-3.0925-4.6770-4.107802.00007.00002.50007：用jacobi和Gauss-Scidel迭代求解下面的方程组，输出每一步的误差xk-x*；⎧5x1-x2-3x3=-2⎪⎨-x1+2x2+4x3=1⎪-3x+4x+15x=10123⎩Jacobi迭代通用程序：functionjacobi(A,b,x0,ep)ifnargin==3ep=1.0e-6;elseifnargin<3errorreturnendendD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;whilenorm(y-x0)>=epfprintf('%d步误差\n',n)Error=y-x0;disp(Error)x0=y;y=B*x0+f;n=n+1;enddisp('近似解')disp(y)在命令窗口输入以下命令就可以输出每一步的误差及最优近似解：>>A=[5-1-3;-124;-3415];>>b=[-2110]';>>x0=[000]';>>ep=10e-8;>>jacobi(A,b,x0,ep)近似解Gauss-Scidel迭代法通用程序：functionGaussseidel(A,b,x0,ep)ifnargin==3ep=1.0e-6;elseifnargin<3errorreturnendendD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=(D-L)\U;f=(D-L)\b;y=B*x0+f;n=1;whilenorm(y-x0)>=epfprintf('%d步误差\n',n)Error=y-x0;disp(Error)x0=y;y=B*x0+f;n=n+1;enddisp('x的近似最优解：')disp(y);在命令窗口输入以下命令就可以输出每一步的误差及最优近似解：>>A=[5-1-3;-124;-3415];>>b=[-2110]';>>x0=[000]';>>ep=10e-8;>>Gaussseidel(A,b,x0,ep)x的近似最优解：-1.8032785851349118：取不同的初值用Newton迭代法以及弦截法求解x3+2x2+10x-100=0的实根，列表或者画图说明收敛性；Newton迭代法通用程序：%使用说明f为符号表达式，x0为初始值functionNewton(f,x0)symsx;x=x0;x1=x0-eval(f/diff(f));while(norm(x1-x0)>0.000001)x=x1;xk=x1-eval(f/diff(f));x0=x1;x1=xk;disp(xk);enddisp('xk=')disp(xk)弦截法通用程序：%使用说明f为符号表达式，x0,x1为初始值functionXuanjiefa(f,x0,x1)symsx;x=x0;f1=eval(f);x=x1;f2=eval(f);x2=x1-(f2/(f2-f1))*(x1-x0);while(norm(x2-x1)>0.000001)x=x1;f1=eval(f);x=x2;f2=eval(f);xk=x2-(f2/(f2-f1))*(x2-x1);x1=x2;x2=xk;enddisp('xk=')disp(xk)Newton法：弦截法：从上述两个表的结果我们可以对照得出，在不同的初始值的条件下运行程序都具有收敛性，但是不同的初始值有不同的收敛速度；9：用二分法求解方程excosx+2=0在区间[0,4π]上所有根。functioner_fen_fa(f,a,b)symsx;e=0.0000001;while((b-a)>e)m=(b+a)/2;x=a;f1=eval(f);x=m;f2=eval(f);iff1*f2>=0a=m;elseb=m;endendc=(a+m)/2;disp(c);根据函数的特点将其定义域分为4等分分别为[0,π]，[π,2π]，[2π,3π]，[3π,4π]分别运行程序：>>symsx>>f=exp(x)*cos(x)+2;>>a=0;>>b=pi;>>er_fen_fa(f,a,b近似零点：1.8807>>symsx>>f=exp(x)*cos(x)+2;a=pi;>>b=2*pi;>>er_fen_fa(f,a,b)近似零点：4.6941>>symsx>>f=exp(x)*cos(x)+2;>>a=2*pi;>>b=3*pi;>>er_fen_fa(f,a,b)近似零点：7.8548>>symsx>>f=exp(x)*cos(x)+2;>>a=3*pi;>>b=4*pi;>>er_fen_fa(f,a,b)近似零点：10.9955;x2=4.6941;x2=7.8548;x4=10.9955由上述程序可知其全部解为x1=1.880710：考虑函数f(x)=sin(πx),x∈[0,1]。用等距节点作f(x)的Newton插值，求出插值多项式以及f(x)的图像，观察收敛性。Newton插值求插值多项式的通用程序：functionP=Nweton_xhazhi(X,f)symsx;p=0;[~,n]=size(X);fori=1:n-1F=0;w=1;sum=1;forl=1:i+1;w=w*(x-X(l));endforj=1:i+1x=X(j);F=F+f(j)/(eval(diff(w)));endfort=1:isymsx;sum=sum*(x-X(t));endp=p+F*sum;endP=p+f(1);Disp(‘P=’);用等距节点作f(x)的Newton插值多项式的通用程序：X=0:0.1:1;f=sin(pi*X);O=Nweton_xhazhi(X,f);Y=simple(O);disp(Y);plot(X,f)P=画图程序：X=0:0.1:1;x1=0:0.0001:1;f1=sin(x1*pi);symsxx0=0:0.001:1;f=sin(pi*X);O=Nweton_xhazhi(X,f);x=x0;Y=eval(O);plot(x0,Y,'r');holdonplot(x1,f1);由上图可知：红色图形是利用插值多项式画出的，蓝色图像是利用原函数画出的，将两者对比我们可以得出结论：该插值多项式具有收敛性11：对函数f(x)=1,x∈[-5,5]，取不同的节点数n,用等距节点作Lagrange1+x2插值，观察Runge现象。Lagrange插值通用程序：%X为节点，y为节点对应的函数值，x0为插值节点functionLagrange(X,y,x0)symsx;[~,n]=size(X);sum=0;fori=1:nsum=sum+y(i)*L(X,i);endx=x0;p=eval(sum);plot(x0,p,'r-')holdonplot(X,y,'b-');functionl=L(X,k)symsx;sum1=1;sum2=1;[~,n]=size(X);if(k==1)fori=k+1:nsum1=sum1*(x-X(i));sum2=sum2*(X(k)-X(i));l=sum1/sum2;endelsefori=1:k-1sum1=sum1*(x-X(i));sum2=sum2*(X(k)-X(i));endfori=k+1:nsum1=sum1*(x-X(i));sum2=sum2*(X(k)-X(i));endl=sum1/sum2;end取不同的节点数n,用等距节点作Lagrange插值，观察Runge现象:当n=10;X=-5:0.1:5;f=1./(1+X.^2);x0=-5:5:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)当n=20;X=-5:0.1:5;f=1./(1+X.^2);x0=-5:0.5:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)当n=50;X=-5:0.1:5;f=1./(1+X.^2);x0=-5:0.2:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)n=100;X=-5:0.1:5;f=1./(1+X.^2);x0=-5:0.1:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)-5-4-3-2-1012345由以上几个图形我们可以得知，插值节点越多时，插值多项式的收敛性越好；12：令f(x)=e3xcos(πx)，考虑积分⎰2π0f(x)dx。区间分为50，100，200，500，1000等，分别用复合梯形以及复合Simpson积分公式计算积分值，将数值积分的结果与精确值比较，列表说明误差的收敛性。复合梯形公式通用程序：%f为被积函数，n为积分区间等分的数目，a,b分别为积分下上限functiontixing(f,n,a,b)symsx;sum=0;xk=a:(b-a)/n:b;fori=2:n-1x=xk(i);sum=sum+eval(f);symsx;endx=a;f1=eval(f);x=b;f2=eval(f);T=(b-a)/(2*n)*(f1+f2+2*sum);disp(T)复合sipson公式：functionSimpson(f,n,a,b)symsx;sum=0;sum1=0;xk=a:(b-a)/n:b;fori=1:n-1x=(xk(i)+xk(i+1))/2;sum=sum+eval(f);symsx;endfori=2:n-1x=xk(i);sum1=sum1+eval(f);symsx;endx=a;f1=eval(f);x=b;f2=eval(f);T=(b-a)/(6*n)*(f1+f2+2*sum1+4*sum);disp(T)利用matlab内置函数int()求出该积分式的真是值为：>>symsx>>int(exp(3*x)*cos(pi*x),0,2*pi)ans=3.5232e+07利用符合simpson公式计算如下表：利用复合梯形公式计算如下表：由上述两个表格的值我们可以得出，其误差随着n的增大近似值与真实值之间的误差逐渐减少，说明了误差具有收敛性。13：分别用2点,3点以及5点的Gauss型积分公式计算如下定积分：(1)⎰1x2-x2-1dx(2)⎰2πsinxdxx解：利用Gauss型求积公式得出：2点Gauss积分公式求解时，n=1;查找上表得x0,x1,A0,A1，把他们都代入，可得⎰1x2-x2-1dx≈(0.5773502692)2-(0.5773502692)2+(-0.5773502692)2-(-0.5773502692)2=0.81653点Gauss积分公式求解时，n=2；查找上表可得x0,x1,x2,A0,A1,A2，把他们都代入，可得：⎰1x2-x2-1dx≈0.5555555556⨯(0.7745966692)2-(0.7745966692)222+0.5555555556⨯(-0.7745966692)-(-0.7745966692)=1.05415点Gauss积分公式求解时，n=4;查找上表可得x0,x1,x2,x3,x4,A0,A1,A2,A3,A4，把他们都代入，可得：⎰1x22-1-x0.53846931012dx≈2⨯0.2369268851⨯=1.24952+2⨯0.4786286705⨯-0.53846931012⎰2π02πsinx1sin(π(t+1))sinx2π-02π+0dx，令x=t+=πt+π则⎰dx=⎰，于是利用0-1x22xt+1上表计算Gauss积分2点Gauss积分公式求解时，n=1;查找上表得x0,x1,A0,A1，把他们都代入，可得1sin(π(t+1))sinxsin(π(1+0.5773502692))dx==+⎰0x⎰-1t+11+0.5773502692sin(π(1-0.5773502692))=1.68121-0.57735026922π3点Gauss积分公式求解时，n=2；查找上表可得x0,x1,x2,A0,A1,A2，把他们都代入，可得：1sin(π(t+1))sinxsin(π(1+0.7745966692))dx==0.5555555556⨯⎰0x⎰-1t+11+0.7745966692sin(π(1-0.7745966692))+0.5555555556⨯=1.39951-0..77459666922π5点Gauss积分公式求解时，n=4;查找上表可得x0,x1,x2,x3,x4,A0,A1,A2,A3,A4，把他们都代入，可得：1sin(π(t+1))sinx51⨯⎰0xdx=⎰-1t+1dt=0.2369268859))sin(π(1-0.9061798459))⎫⎛sin(π(1+0.90617984+⎪+1+0.90617984591-0.9061798459⎝⎭2π01))sin(π(1-0.5384693101))⎫⎛sin(π(1+0.538469310.4786286705⨯+⎪1+0.53846931011-0.5384693101⎝⎭=1.418114：考虑微分方程初值问题：1⎧dx=(tx-x2)⎪2⎨dt(t+1)⎪x(0)=2⎩分别用Euler法，改进的Euler法，Runge-Kutta法求解该方程。分别去步长为0.1，0.01，0.001，计算到x(1),画图说明结果。Euler法：%f为被积表达式，h为步长，ab为t的取值范围，x0为初始值functionEuler(f,h,a,b,x0)symstxN=abs(b-a)/h;x=x0;Y=zeros(1,N+1);T=zeros(1,N+1);Y(1)=x0;forn=0:Nt=a+n*h;T(n+1)=t;if(n<N)un=x+h*eval(f);x=un;Y(n+2)=un;endenddisp(un);plot(T,Y,'*-')N=0.1;x(1)=1.072393608200700改进的Euler法：%f为被积表达式，h为步长，ab为t的取值范围，x0为初始值functionGai_Jin_Euler(f,h,a,b,x0)symstxN=abs(b-a)/h;x=x0;Y=zeros(1,N+1);T=zeros(1,N+1);Y(1)=x0;forn=0:Nt=a+n*h;T(n+1)=t;f1=eval(f);t=a+(n+1)*h;un1=x+h*f1;x=un1;f2=eval(f);ifn<Nx=x0;un=x+h/2*(f1+f2);Y(n+2)=un;x0=un;endx=un;enddisp(un);plot(T,Y,'*-');Runge-Kutta法：矩阵分析在-------机械振动中的应用摘要：随着科学技术的迅速发展，古典的线性代数知识已不能满足现代科技的需要，矩阵的理论和方法业已成为现代科技领域必不可少的工具。诸如数值分析、优化理论、微分方程、概率统计、控制论、力学、电子学、网络等学科领域都与矩阵理论有着密切的联系，甚至在经济管理、金融、保险、社会科学等领域，矩阵理论和方法也有着十分重要的应用。本文采用了矩阵论中所学的矩阵相似变换、矩阵正交化及特征方程等相关知识，对多自由度系统的自振动的运动微分方程进行了研究分析，引入正则坐标并采用坐标变化法求得了振动系统的自由响应。关键词：多自由度系统，正则坐标，自由响应一、引言20世纪60年代，随着计算机技术的进步，航空航天技术和综合自动化的发展需要，对于复杂的机械结构特性分析也越来越重要。而对于像航天器等复杂的机械结构需要用更多的自由度来描述，多自由度系统的振动方程式二阶常微分方程组。建立系统方程是振动分析的前提，但随着自由度的增多，所建立的系统运动微分方程也越来越复杂，对于离散系统运用牛顿第二定律的方式来对方程进行求解也越来越困难，为此发展了柔度系数法和刚度系数法，而拉尔朗日方程是建立系统控制方程的最通用方法，他使用功、能和广义力等物理量，得到了完全刻画系统的最少方程。本文只考虑阻尼矩阵能够被无阻尼振形矩阵对角化的情形，分析其基本理论方程，并用实例进行论证求解。二、多自由度系统的自由振动理论本文主要对多自由度系统的自由振动进行求解，在介绍多自由度系统的振动之前，先介绍单自由度无阻尼的自由振动以便了解机械振动理论的基本原理。1.单自由度无阻尼系统的自由振动图1单自由度无阻尼系统对于单自由度系统而言，当系统受到激励时，根据牛顿第二定律，可以列出的运动微分方程为：(1.1其中，m为物体的质量；k为弹簧的刚度；为物体的加速度；x为弹簧的伸缩量。该方程是一个二阶齐次线性常系数微分方程。这为之后的多自由度系统的运动分析提供了理论基础。2.多自由度无阻尼系统的自由振动多自由度系统和单自由度系统的振动特性是有区别的。单自由度系统受初始扰动后，按系统的固有频率作简谐振动。多自由度系统有多个固有频率，当系统按某一个固有频率作自由振动时，各独立坐标在振动过程中相互关系是固定的，这个关系叫振幅比，也叫作主振型或模态。主振型是多自由度系统以及弹性体振动的重要特征。本文主要目的是通过无阻尼自由振动系统来介绍多自由系统的固有频率和振型，它们是多自由振动系统的重要特征。在无阻尼情况下，系统的自由振动微分方程可以表达为：(1.2在单自由度系统中,我们得到无阻尼自由振动解为正弦函数或余弦函数,不失一般性。对于多自由度系统振动解可设为：(1.3列向量和ω均为待定复常数。若系统是振动的，则解必为实数。将式(1.3代入(1.2，得到下列代数齐次方程组：(1.4上面的方程组存在非零解的充分必要条件是系数行列式为零，即：(1.5式(1.5为系统的特征方程，具体写出为：(1.6上式左端的行列式展开后是关于的n次代数多项式：(1.7称为特征多项式，由式(1.6或(1.7可解出n个称为特征值或特征根，将其按升序排列为：显然特征值仅取决于系统本身的刚度和质量参数。这n个特征值在大多数情况下互不相等且不为零，重根的零根说明系统有刚体运动。有零根和情况本书不再讨论，有兴趣的读者可参考相关的线性代数和振动理论书籍。在求得特征值后．把某一个代回式(1.4，可求对应的列向量。由于式(1.4的系数矩阵不满秩，在没有重根和零根情况下只有(n-1个是独立的，故只能求出列向量中各元素、、…的比例关系。我们去掉其中不独立的某一式(例如最后一式，并将剩下的n-1个方程式中某一相同的项(如项移到等式右边，可得代数方程组：我们去掉其中不独立的某一式(例如最后一式，并将剩下的n-1个方程式中某一相同的项(如项移到等式右边，可得代数方程组：(1.8解上面的方程，可得到用表达的解、…，显然都与的值成比例。我们可将这些比例常数用表示，并补充，可得列向量，则有：(1.9列向量是确定的常数，反映列向量中各数的比例关系，叫作特征向量。同比例放大或减小特征向量并不改变其比例关系，所以应用时常根据需要来放大或减小特征向量。不失一般性，我们可在式(1.9中用待定复常数取代，式(1.9可写为：(1.10这样，当成比例变化时，有相应的变化，对应不同的特征值，可得到不同的特征向量。对应于n个特征值可得n个特征向量…，且每一个特征向量都满足式(1.4。对于一个振动系统，特征值就是系统的固有频率，特征值相对应的特征向量就是系统的振形。显然，对应于n个固有频率可得n个振形…。我们将在后面论述。显然，将及代入式(1.3，可得n组满足方程(1.2的解，将这些解相加，可得多自由度系统自由振动的一般解为：(1.11其中2n个待定常数由系统运动的初始位移和初始速度确定。如果系统在某一特殊的初始条件下，使得待定常数中只有≠0，则式(1.11所表示的系统运动方程只保留第k项:(1.12多自由度系统振动一般解的方程可表达为：(1.13这时整个系统按圆频率、振幅比作同步简谐运动。振幅分别为，振幅之间都保持固定不变的比值。因此特征向量完全确定了系统按固有频率振动时的形态，所以特征向量就是按相应固有频率振动时的振型向量，对应的特征向量称为它的第阶主振型或主模态，相应的振动叫主振动。在振动过程中，一般还会产生其它阶主振动。对于一个n自由度系统，一般可以找到n个固有频率，以及相应的n个主振型。我们把各阶主振型组成的矩阵叫做振型矩阵：(1.14三、三自由度系统自由响应求解三自由度的弹簧-质量系统如图11所示，设t=0时。求振系的自由响应。图2三自由度无阻尼系统解：第一步，建立振动微分方程，由刚度法可建立该振系的微分方程第二步，求固有频率和振型。系统的,，故系统矩阵将[S]代入振型方程得故频率方程为由上式解得三个特征值为对应的固有频率为将代入振型方程a消去公因子，并令=1，则有由上式解得，对，做同样的处理，得到相应的振型为第三步，求振型矩阵与正则矩阵。振型可知，振型矩阵即可确定为求正则振型矩阵，需先求出各阶主质量再求出各阶正则振型由正则振型即可构成正则振型矩阵第四步，用正则坐标变换可得到用正则坐标表示的独立方程(i=1,2,3第五步，把初始条件变换到正则坐标上，若将式子两端左乘则有因,即第六步，求振系在正则坐标下的响应。而方程的一般解为代入正则坐标表示的初始条件，，第七步，把正则坐标的响应再变回到物理坐标系下。利用坐标变换式得四、结论本文在研究多自由度系统的自由振动时，将模型简化成无阻尼系统，并使用了特征方程，矩阵逆变换等相关知识进行求解。得出了多自由系统在激励下的自由响应。其实在实际问题中，系统几乎都是有阻尼的，此时，所列出的运动微分方程也更加复杂，所需要用到的矩阵论的知识也更多。在计算机发展和普及的前提下，矩阵论理论的重要性越来越明显，应用也越来越广泛。当然，研究梁单元的振动情况只是矩阵论理论应用领域的一个小方面。但是，这足以说明用矩阵论力量和方法可以方便地解决现代工程技术中的各种问题，它表述简洁，便于进行研究，已经越来越成为从事科学研究和工程设计科技人员的首选工具。参考文献[1]李新.何传江.矩阵理论及其应用[M].重庆大学出版社，2005年8月.[2]李有堂.机械系统动力学[M].国防工业出版社，2021[3]陈天福，冯贤贵.材料力学[M].重庆大学出版社，2006万方数据万方数据万方数据万方数据5.2数值仿真实验(1实验方法模拟透视成像摄像机,其内参数为K=被观测目标是200个3D点,随机分布在边长为360单位的立方体中,立方体中心距离摄像机光心600个单位。两台摄像机相对运动参数为:平移向量t=[40,o,o]7、旋转轴方向为Eo,0,1]7、旋转角度为10。。根据摄像机成像模型,将3D点投影到两个成像平面上,然后叠加零均值高斯噪声,其方差o(noisesigma从0到1变化,步长0.1,同时错误数据所占百分比r(outlierratio从o%到50%变化,步长10%。(2实验结果将几种典型方法计算结果的误差分布曲面显示在图3中。图3基础矩阵计算的仿真实验5.3真实图像实验(1实验方法图4真实图像对和初始特征点匹配首先.采用改进的Harris角点检测算法分别在左右两幅图像中检测特征点,然后,使用零均值相关系数(ZNCC匹配方法获得初始匹配;最后,将匹配结果作为输入数据计算两幅图像的基础矩阵。(2实验结果这里给出3组图像的实验结果。其中,Mars为室外自然景物(512X512,图4(a,Inria为室外人造建筑(512×256,图4(b,Desk为室内人造物(800X600,图4(c。图4中显示了初始匹配特征点及其视差向量。图中还显示了通过LMedS+MEst方法得到的部分极线。为便于比较,将数值仿真实验结果同真实图像的实验结果一并在表1中列出。表1基础矩阵计算的实验结果1-LinEig・2一LinLS,3-IterEig,4一herLS。5-GradEig。6-GradLS,7-NonParam.8-M—Est,9-RANSAC,10一LMedS,1卜M-Est+LMedSd=10o.00001000置:o.9931.0160.989r=u.J422=10‘.吉1.885f!?’36.55435.50232.07232.08144.45532.08131.3672.142o.990I.8861.839Man13.11313.20512.68612.88712.69112.88712.9000.8540.7590.872o.692Inria5.9795.9688.9317.4258.7717.4255.960Desk31.56829.88233.03830.41731.27739.4】629671l】9501178366613636分析与结论以仿真实验以及真实图像实验结果为依据,本文对不同的F矩阵计算方法有如下分析和结论。(1线性方法:如果特征点定位较精确并且没有错误匹配,那么可以得到很好的结果。但是,这类方法对错误匹配的鲁棒性非常差。(2非线性迭代方法:在一定程度上能够处理定位噪声的影响,但是实验结果表明改善效果不明显,同时,当存在错误匹配时,效率很低。(3参数空间优化方法:比线性方法和迭代方法好,但是同样不能处理错误匹配问题。(4鲁棒方法:能够同时处理数据噪声和错误匹配。当存在错误匹配时,M.Kstimator方法的性能有所降低,特别是对于Desk图像,由于初始匹配错误率较大导致聊s较大。这一现象证实了M—Estimator对定位噪声是鲁棒的,但是对错误数据的鲁棒性较差。相比之下,RANSAC和LMedS体现出较好的品质。(5如果使用LMedS方法预先剔除错误匹配,然后使用M—Estimator,可以得到最好的结果。(6同线性最dx_-乘优化方法比较,基于特征分析技术的优化方法效果更好。(下转第289页・247・●印√叫j筋弱●OO∞O∞OO陌一万方数据Ｇｒａｐｈｉｃｓ，１９８４，１８（３）：卜１０２２３—２３１［５］ＰｒｕｓｉｎｋｉｅｗｉｃｚｏｆＰ，ＬｉｎｄｅｎｍｅｙｅｒＡ・ｅｔａ１．ＴｈｅＡｌｇｏｒｉｔｈｍｉｃＢｅａｕｔｙ［１１３ＨａｎｒａｈａｎＰ，ＬａｗｓｏｎＪ．Ａｌａｎｇｕａｇｅｆｏｒｓｈａｄｉｎｇａｎｄｌｉｇｈｔｉｎｇｃａｌ—Ｐｌａｎｔｓ［Ｍ］．Ｓｐｒｉｎｇｅｒ－Ｖｅｒｌａｇ．１９９０ｅｕｌａｔｉｏｎｓ［Ｊ］．Ｃｏｍｐｕｔｅｒ［１２］ＰｈａｒｒＧｒａｐｈｉｃｓ，１９９０，２４（４）：２８９—２９８［６］王辉．基于Ｌ系统的虚拟园林观赏树木生长建模研究［Ｄ］．长沙：中南林业科技大学，２００６：７－９Ｍ，ＨａｎｒａｈａｎＰ．Ｇｅｏｍｅｔｒｙｃａｃｈｉｎｇｆｏｒｒａｙ—ｔｒａｄｎｇｄｉｓ一Ｗｏｒｋｓｈｏｐｐｌａｃｅｍｅｎｔｍａｐｓ［Ｃ］∥Ｐｒｏｃ．ｏｆｔｈｅ７ｔｈＥｕｒｏｇｒａｐｈｉｃｓ［７３［８３ＨｅａｍＤ．ＢａｋｅｒＭＰ．计算机图形学［Ｍ］．蔡士杰。吴春熔，孙正ｏｎＲｅｎｄｅｒｉｎｇ．Ｐｏｒｔｏ．Ｐｏｒｔｕｇａｌ．１９９６：３１—４０兴，译．电子工业出版社，２００２ＦａｒｍＧ．ＣｕｒｖｅｓａｎｄＳｕｒｆａｃｅｓｆｏｒＣｏｍｐｔｅｒＡｉｄｅｄＧｅｏｍｅｔｒｉｃＤｅ－ｓｉｇｎ：ＡＰｒａｃｔｉｃａｌＧｕｉｄｅ，４ｔｈＰｒｅｓｓ・１９９６［１３］ＳｅｈａｕｆｔｅｒＧ．Ｐｒｉｇｌｉｎｇｅｒ札Ｅｆｆｉｄｅｎｔｉｍａｇｅｄｉｓｐｌａｃｅｍｅｎｔｍａｐｐｉｎｇｂｙｗａｒｐｉｎｇ［Ｃ－］｝｝Ｐｒｏｃｏｆｔｈｅ１０ｔｈＥｕｒｏｇｒａｐｈｉｃｓＷｏｒｋｓｈｏｐＥｄｉｔｉｏｎ［Ｍ］．Ｂｏｓｔｏｎ：ＡｃａｄｅｍｉｃｏｎＲｅｎｄｅｒｉｎｇ．Ｇｒａｎａｄａ，Ｓｐａｉｎ，Ｅｕｒｏｇｒａｐｈｉｃｓ，１９９９：１７５—１８６［１４］付恺。李春霞，杨克俭，等．ＢｕｍｐＭａｐｐｉｎｇ原理及在ＯｐｅｎＧＬ下的实现口］．交通与计算机，２００４，２２（２）：５４—５７［１５］杨刚．全景图拼接算法的设计与实现口］．重庆工学院学报；自然Ｇｒａｐｈｉｃｓ，１９８４，１８（３）：［９］李钢，刘华明．基于ＮＵＲＢＳ的扫描曲面造型方法的研究［Ｊ］．机械研究与应用，２０００，１３（３）：５－６［１０］ＣｏｏｋＲＬＳｈａｄｅｔｒｅｅｓＥＪ］．Ｃｏｍｐｕｔｅｒ科学版，２００７，２１（１９）：１０７—１１０（上接第２４７页）结束语Ｆ矩阵是许多计算机视觉应用中的重要参数，其计算的准确性决定了后续的处理步骤能否成功。本文将常见的Ｆ矩阵估计方法划分为线性法、非线性法和鲁棒法３大类，共计１１种，通过仿真数据和真实图像实验对各自的性能进行了评估。实际应用中，由于视觉系统所采用的图像特征定位及匹配方法各有差异，可根据具体情况合理选择本文介绍的Ｆ矩阵计算方法，以实现高精度的３Ｄ信息获取。Ｅ１１］ＦａｎｇｅｒａｓｆｅｒｅｎｃｅＤ，ＬｕｏｎｇＱＴ，ＭａｙｂａｎｋＳＪ．Ｃａｍｅｒａｓｅｌｆ－ｃａｌｉｂｒａｔｉｏｎ：ｔｈｅｏｒｙａｎｄｏｎｅｘｐｅｒｉｍｅｎｔ［Ｃ］？ｊＰｒｏｃｅｅｄｉｎｇｓｏｆｔｈｅＥｕｒｏｐｅａｎＣＯｎ－ＣｏｍｐｕｔｅｒＶｉｓｉｏｎ．ＳａｎｔａＭａｒｇｈｅｆｉｔａＬ．１９９２：３２１—３３４［１２３ＬｏｕｒａｋｉｓＭＩＡ。ＤｅｆｉｃｈｅＲＣａｍｅｒａｓｅｌｆ－ｃａｌｉｂｒａｔｉｏｎｕｓｉｎｇｔｈｅｓｉｎｇｕｌａｒｖａｌｕｅｄｅｃｏｍｐｏｓｉｔｉｏｎｏｆｔｈｅｆｕｎｄａｍｅｎｔａｌｔｏｍａｔｒｉｘ：Ｆｒｏｍｐｏｉｎｔｃｏｒｒｅｓｐｏｎｄｅｎｃｅｓ３Ｄ１９９９ｍｅａｓｕｒｅｍｅｎｔｓ［Ｒ］．３７４８．ＩＮＲＩＡＳｏｐｈｉａ－Ａｎｔｉｐｏｌｉｓ，Ａｕｇｕｓｔ［１３３ＬｏｕｒａｋｉｓＭＩＡ。ＤｅｆｉｃｈｅＲＣａｍｅｒａｓｅｌｆ—ｃａｌｉｂｒａｔｉｏｎｕｓｉｎｇｔｈｅＫｒｕｐｐａｅｑｕａｔｉｏｎｓａｎｄｔｈｅＳＶＤｏｆｔｈｅｆｕｎｄａｍｅｎｔａｌｍａｔｒｉｘ：Ｔｈｅｃａｓｅ参考文献［１］马颂德，张正友．计算机视觉一计算理论与算法基础［Ｍ］．北京：科学出版社，１９９８［２３ＬｕｏｎｇＱ－Ｔ，ＦａｕｇｅｒａｓＯ１３．Ｔｈｅｆｕｎｄａｍｅｎｔａｌｍａｔｒｉｘ：Ｔｈｅｏｒｙ．ａｌ—ｇｏｒｉｔｈｍｓａｎｄｓｔａｂｉｌｉｔｙｐｕｔｅｒｏｆｖａｒｙｉｎｇｉｎｔｒｉｎｓｉｃｐａｒａｍｅｔｅｒｓ［Ｒ］．３９１１．ＩＮＲＩＡＳｏｐｈｉａ－Ａｎｔｉｐｏｌｉｓ，Ｍａｒｓ．２０００［１４］Ｈａｒｔｌｅｙ［１５］ＨａｒｔｌｅｙＲｌ。Ｚｉｓｓｅｒｍａｎ八Ｍｕｌｔｉｐｌｅｖｉｅｗｇｅｏｍｅｔｒｙｉｎｃｏｍｐｕｔｅｒｖｉｓｉｏｎ［Ｍ］．ＣａｍｂｒｉｄｇｅＵｎｉｖｅｒｓｉｔｙＰｒｅｓｓ，２０００ＲＬＩｎｄｅｆｅｎｓｅｏｆｔｈｅ８－ｐｏｉｎｔａｌｇｏｒｉｔｈｍ［Ｃ］∥Ｐｒｏｃｅｅ－ｏｎａｎａｌｙｓｉｓ［Ｊ］．ＩｎｔｅｒｎａｔｉｏｎａｌＪｏｕｒｎａｌｏｆＣｏｍ—ｄｉｎｇｓｏｆｔｈｅ５ｔｈＰｒｏｃ．ＩｎｔｅｍａｔｉｏｎａｌＣＯｎｆｅｒｅｎｃｅＣｏｍｐｕｔｅｒＶｉ－Ｖｉｓｉｏｎ，１９９６，Ｉ（１７）：４３—７６ｓｉｏｒｌ．Ｂｏｓｔｏｎ，／ＶＩＡ：ＩＥＥＥ１０７０ＣｏｍｐｕｔｅｒＳｏｃｉｅｔｙＰｒｅｓｓ，１９９５：１０６４—［３３Ｌｏｎｇｕｅｔ－ＨｉｇｇｉｎｓＨｃＡｃｏｍｐｕｔｅｒａｌｇｏｒｉｔｈｍｆｏｒｒｅｃｏｎｓｔｒｕｃｔｉｎｇａｓｃｅｎｅｆｒｏｍｔｗｏｐｒｏｉｅｅｔｉｏｎｓ［Ｊ］．Ｎａｔｕｒｅ，１９８１，２９３：１３３—１３５ｓｔｒｕｃｔｕｒｅ［１６］ＦａｕｇｅｒａｓｐｏｉｎｔＱ，ＭｏｕｒｒａｉｎＲＯｎｔｈｅｇｅｏｍｅｔｒｙａｎｄａｌｇｅｂｒａｏｆｔｈｅｎ［４］ＭｅｓｔｅｒＲｔｉｏｎＯｎｔｈｅｍａｔｈｅｍａｔｉｃａｌｏｆｄｉｒｅｃｔｉｏｎａｎｄｏｎｎｌｏ—ａｎｄｌｉｎｅｃｏｒｒｅｓｐｏｎｄｅｎｃｅｓｂｅｔｗｅｅｎｉｍａｇｅｓ［Ｃ］ｆ｝Ｐｒｏｃｅｅ－ｅｓｔｉｍａｔｉｏｎ［ｑｆ｝３ｒｄＩｎｔｅｒｎａｔｉｏｎａｌＳｙｍｐｏｓｉｕｍＰｈｙｓｉｃｓｉｎｄｉｎｇｓｏｆｔｈｅ５ｔｈＰｒｏｃＩｎｔｅｒｎａｔｉｏｎａｌＣｏｎｆｅｒｅｎｃｅｏｎＣｏｍｐｕｔｅｒＶｉ—ＳｉｇｎａｌａｎｄＩｍａｇｅＰｒｏｃｅｓｓｉｎｇ（ＰＳＩＰ）．Ｇｒｅｎｏｂｌｅ。Ｆｒａｎｃｅ，Ｊａｎｕａｒｙ２００３ｓｉｏｎ．Ｂｏｓｔｏｎ，ＭＡ：ＩＥＥＥＣｏｍｐｕｔｅｒＳｏｄｅｔｙＰｒｅｓｓ，１９９５：９５１－９５６［１７］Ｚｈａｎｇｚ，ＤｅｒｉｃｈｅＲ．ＦａｕｇｅｒａｓＯ，ｅｔａ１．Ａｒｏｂｕｓｔｔｅｃｈｎｉｑｕｅｆｏｒ［５３ＰｏｌｌｅｆｅｙｓｓｐｉｔｅＭ，ｅｔａ１．Ｓｅｌｆ－ｃａｌｉｂｒａｔｉｏｎａｎｄｍｅｔｒｉｃｒｅｃｏｎｓｔｒｕｃｔｉｏｎｉｎｍａｔｃｈｉｎｇｔｗｏｕｎｃａｌｉｂｒａｔｅｄｉｍａｇｅｓｔｈｒｏｕｇｈｔｈｅｒｅｃｏｖｅｒｙｏｆｔｈｅｏｆｖａ叫ｌａｇａｎｄｕｎｋｎｏｖｃｎｉｎｔｅｒｎａｌｃａｍｅｒａｌｐａｒａｍｅｔｅｒｓ［Ｊ］．ｕｎｋｎｏｗｎｅｐｉｐｏｌａｒｇｅｏｍｅｔｒｙ［Ｊ］．ＡｒｔｉｆｉｃｉａｌＩｎｔｅｌｌｉｇｅｎｃｅＪｏｕｒｎａｌ，ＩｎｔｅｍａｔｉｏｎａｌＪｏｕｒｎａｌｏｆＣｏｍｐｕｔｅｒＶｉｓｉｏｎ，１９９９，３２（１）：７－２０１９９５，７８（１０）：８７－１１９［６３ＺｈａｎｇＺ，ＤｅｆｉｃｈｅＲ，ＦａｕｇｅｒａｓＯ，ｅｔａＬＡｒｏｂｕｓｔｔｅｃｈｎｉｑｕｅｆｏｒ［１８３Ｚｈａｎｇｚ．Ｄｅｔｅｒｍｉｎｉｎｇｔｈｅｅｐｉｐｏｌａｒｇｅｏｍｅｔｒｙａｎｄｉｔｓｕｎｃｅｒｔａｉｎｔｙ：Ａｍａｔｃｈｉｎｇ１１３９ｔｗｏｕｎｃａｌｉｂｒａｔｅｄｉｍａｇｅｓ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｏｎｒｅｖｉｅｗ［Ｊ］．ＩｎｔｅｒｎａｔｉｏｎａｌＰａｔｔｅｒｎＡｎａｌｙｓｉｓａｎｄＭａｃｈｉｎｅＪｏｕｒｎａｌｏｆＣｏｍｐｕｔｅｒＶｉｓｉｏｎ，１９９８，２７Ｉｎｔｅｌｌｉｇｅｎｃｅ，１９９５，８２（２）：１１２９—（２）：１６１－１９５［１９］ＺｈａｎｇＺ．ＰａｒａｍｅｔｅｒＥｓｔｉｍａｔｉｏｎＴｅｃｈｎｉｑｕｅｓ：Ａｐｌｉｃａｔｉｏｎｔｏｔｕｔｏｒｉａｌｗｉｔｈａｐ－［７３ＨａｒｔｌｅｙＲ，ＧｕｐｔａＲ。ＣｈａｎｇＴ－Ｓｔｅｒｅｏｆｒｏｍｕｎｃａｌｉｂｒａｔｅｄｃａｍｅｒａｓｃｏｎｉｃｆｉｔｔｉｎｇ［Ｊ］．ＩｍａｇｅａｎｄＶｉｓｉｏｎＣｏｍｐｕｔｉｎｇＪｏｕｒ—［ｃ］∥ＰｒｏｃｅｅｄｉｎｇｓＶｉｓｉｏｎａｎｄｏｆｔｈｅＩｎｔｅｒｎａｔｉｏｎａｌＣｏｎｆｅｒｅｎｃｅｏｎＣｏｍｐｕｔｅｒ［２０３ｈａｌ，１９９７，１５（１）：５９—７６ＦｉｓｃｈｌｅｒＰａｔｔｅｒｎＲｅｅｏｇｎｉｔｉｏｒＬＩＥＥＥ．ＵｒｂａｒｍＣｈａｍｐａｉｇｎ，ＩＬＭＡ，ＢｏｌｌｅｓＲｃＲａｎｄｏｍｓａｍｐｌｅｃｏｎｓｅｎｓｕｓ：Ａｐａｒａ—［８３ＭｏｈｒＲ，ＢｏｕｆａｍａＢ，ＢｒａｎｄＰ．Ａｃｃｕｒａｔｅｐｒｏｊｅｃｔｉｖｅｒｅｃｏｎｓｔｒｕｃｔｉｏｎｄｉｇｍｆｏｒｍｏｄｅｌｆｉｔｔｉｎｇｗｉｔｈａｐｐｌｉｃａｔｉｏｎｓ／ｏｉｍａｇｅａｎａｌｙｓｉｓａｎｄ［ｃ］∥ＭｕｎｄｙＪ．ＺｉｓｓｅｒｍａｎＡ，ｅｄｓ．ＡｐｐｌｉｃａｔｉｏｎｓｏｆＩｎｖａｒｉａｎｃｅｉｎａｕｔｏｍａｔｅｄｃａｒｔｏｇｒａｐｈｙ［Ｊ］．Ｃｏｍｍ．Ａｓｓｏｃ．Ｃｏｍｐ．Ｍａｃｋ，１９８１，２４（６）：３８１－３９５ＣｏｍｐｕｔｅｒＶｉｓｉｏｎ．ＶｏＩ．８２５ｏｆＬｅｃｔｕｒｅＮｏｔｅｓｉｎｅｅ，Ｂｅｒｌｉｎ：ＳｐｒｉｎｇｅｒＶｅｒｌａｇ，１９９３：２５７—２７６ｃｃ肌ｐｕｔ日Ｓｄｅｎ－［２１］Ｔｏｒｒｅｓ。ＰＨ＆Ｍｏｔｉｏｎｓｅｇｍｅｎｔａｔｉｏｎａｎｄｏｕｔｌｉｅｒｄｅｔｅｃｔｉｏｎ［Ｄ］．［９３Ｚｈａｎｇ乙Ａｎｅｗｍｕｌｔｉｓｔａｇｅａｐｐｒｏａｃｈｔｉｍａｔｉｏｎ：ＦｒｏｍｆｕｎｄａｍｅｎｔａｌｔｏｍｏｔｉｏｎａｎｄｓｔｒｕｃｔｕｒｅＤｅｐｔ．ｏｆ［２幻ＺｕｌｉａｎｉＥｎｇｉｎｅｅｒｉｎｇＳｃｉｅｎｃｅ，ＵｎｉｖｅｒｓｉｔｙｏｆＯｘｆｏｒｄ，１９９５Ｍ．ＫｅｒｍｅｙＣｉｔｓｅｓｓｅｎｔｉａｌｐａｒａｍｅｔｅｒｓｔｏＥｕｃｌｉｄｅａｎｍｏｔｉｏｎｖｉａＳ，ＭａｎｊｕｎａｔｈＢｓＴｈｅＭｕｈｉＲＡＮＳＡＣａｌ—ｍａｔｒｉｘ［Ｒ］．２９１０．１ＮＲｌＡｓｏｐｈｉａ－Ａｎｔｉｐｏｌｉｓ。Ｆｒａｎｃｅａｇｏｒｉｔｈｍａｎｄａｐｐｌｉｃａｔｉｏｎｔｏｄｅｔｅｃｔｐｌａｎａｒｏｎｈｏｍｏｇｒａｐｈｉｅｓ［Ｃ］ｆ｝［１０３ＭａｙｈａｎｋＳＪ，Ｐａｕｇｅｒａｓｏ．Ａｔｈｅｏｒｙｏｆｓｅｌｆ－ｃａｌｉｂｒａｔｉｏｎｏｆｍｏｖ—ＩＥＥＥＩｎｔｅｒｎａｔｉｏｎａｌＣｏｎｆｅｒｅｎｃｅ２００５Ｉｍａｇｅｉｎｇｃａｍｅｒａ［Ｊ］．ＩｎｔｅｒｎａｔｉｏｎａｌＪｏｕｒｎａｌｏｆ８（２）：１２３—１５１ＣｏｍｐｕｔｅｒＰｒｏｃｅｓｓｉｎｇ．ＳｅｐｔｅｍｂｅｒＶｉｓｉｏｎ．１９９２，万方数据・２８９・视觉基础矩阵估计方法的性能比较与分析作者单位：蔡涛，段善旭，李德华，CAITao，DUANShan-xu，LIDe-hua蔡涛,段善旭,CAITao,DUANShan-xu(华中科技大学电气与电子工程学院,武汉,430074，李德华,LIDe-hua(华中科技大学图像识别与人工智能研究所,图像信息处理与智能控制教育部重点实验室,武汉430074计算机科学COMPUTERSCIENCE2021,36(12刊名：英文刊名：年，卷