曲线二乘拟合

2024-07-02

曲线二乘拟合（精选7篇）

曲线二乘拟合篇1

前言

在实验数据的处理、分析时, 数据拟合是经常采用的一种方法。数据拟合的目标是找到能反映变量之间关系的一种表达式, 使其在某种准则下最佳地接近已知数据。其原理有最小二乘法、契比雪夫法等, 且以最小二乘法最为常见和重要。

随着人类认识能力的不断进步以及计算技术的快速发展, 对于变量之间的未知关系, 应用曲线拟合的方法揭示其内在规律具有重要的理论和现实意义。针对最小二乘曲线拟合的有关问题以及相应的MATLAB实现进行探讨。

1 最小二乘曲线拟合

曲线拟合是指:已知n个数据点 (xi, yi) , i=1, 2, ……, n, 其中xi不全相同, 寻求函数f (x;a1, a2, …, am) 的待定参数a1, a2, …, am的一组取值, 使得在这组取值之下, 函数f (x;a1, a2, …, am) 与已知n个数据点整体上最为接近。

最小二乘曲线拟合方法根据已知数据, 首先构造出能够反映含有待定参数的函数f (x;a1, a2, …, am) 与n个数据点 (xi, yi) , i=1, 2……, n偏离程度的函数:

然后应用数学方法求函数J (a1, a2, …, am) 的最小值a1, am2, i…n, am=J (a1, a2, …, am) , 此时a1, a2, …, am的取值就是所求的待定值。这样一组取值使得函数f (x;a1, a2, …, am) 与n个数据点在二次平方和意义下最为接近。

2 最小二乘曲线拟合的MATLAB实现

2.1 线性最小二乘曲线拟合的MATLAB实现

线性最小二乘曲线拟合的含有待定参数的函数形式为:

其中rk (x) 为事先选定的一组已知函数, ak为待定系数, k=1, 2, …, n, m

该方程组可简化为RTRA=RTY, 其中:

当r1 (x) , …, rm (x) 线性无关时, RTR可逆, 方程组RTRA=RTY有唯一解。

根据线性最小二乘拟合理论可得关于待求参数构成的矩阵A的解, 利用MATLAB下命令A=RY, 可以直接求得待求参数, 下面举例说明。

例1下表给出了一组实测数据

已知数据 (xi, yi) 的函数原型为:

试用已知数据进行曲线拟合, 求出待定参数的值。

在MATLAB中输入以下程序:

运行程序, 输出待定参数结果为:

下面画出拟合得到的曲线及已知的数据散点图:

由图1可见, 曲线拟合效果非常好。事实上, 所给数据就是由已知曲线:

2.2 非线性最小二乘曲线拟合

在最小二乘曲线拟合中, 如果待定参数不能全部以线性形式出现, 如指数拟合函数f (x) =a0+a1e-a2·x, 这便是非线性最小二乘曲线拟合。

非线性最小二乘曲线拟合与线性最小二乘曲线拟合的原理没有什么区别, 但是最小二乘的解常常难以通过人的手算实现, 从而制约了该方法的应用。随着计算机技术的进步、专业软件的不断涌现, 这一问题的求解已不再困难。但是, 非线性曲线拟合中初值的选取是一个重要的问题, 目前为止还没有固定的理论或方法给出一般性的结论。

在MATLAB中实现非线性最小二乘曲线拟合有三个常用的命令。

(1) lsqcurvefit () 命令, 其使用格式为x=lsqcurvefit (fun, x0, xdata, ydata) , 其中fun是要拟合的非线性函数, x0是初始参数, xdata, ydata是拟合点的数据, 该函数最终返回系数矩阵。

(2) nlinfit () 命令, 其应用格式为beta=nlinfit (x, y, fun, beta0) , 其中x和y是拟合点数据, fun是被回归 (拟合) 的函数, beat0是初始参数。

(3) lsqnonlin () 命令, 其应用格式为x=lsqnonlin (fun, x0) , 其中fun的定义与lsqnonlin () 函数中fun的定义有差别, x0仍为初始参数向量, 将输出的系数结果放在变量x中。

例2假设已知

并已知该函数原型为y (x) =a1exp (a2x) +0.54exp (a3x) .sin (a4x) , 其中ai为待定系数。

在MATLAB中输入以下命令:

%建立函数原型, 则可以根据他来进行下面的求取系数的计算

运行结果为:

所求得的系数与原式中的系数相近。

如果想进一步提高精度, 则需修改最优化的选项, 函数的调用格式也将随之改变。

在MATLAB中输入以下命令:

%修改精度, 即是修改其限制条件

%两个空矩阵表示系数向量的上下限a', res运行结果为:

下面画出由拟合得到的曲线及已知的数据散点图, 如图2。

3 结论

首先介绍了最小二乘法的原理, 其次对线性最小二乘曲线拟合和非线性最小二乘曲线拟合应用MATLAB进行了具体实现, 使得相关曲线拟合理论更加生动易懂, 这使我们可以在今后的研究和工作中建立曲线模型对相应对象进行曲线拟合, 并应用MATLAB来实现, 从而找到更好更形象的反映变量之间关系的曲线, 也就找到最合适的模型了。

参考文献

[1]薛定宇, 陈阳泉.高等应用数学问题的MATLAB求解[M].北京:清华大学出版社, 2004.

[2]刘琼荪, 何中市, 傅鹂等.数学实验[M].北京:高等教育出版社, 2000.

[3]梁国业, 廖健平.数学建模[M].北京:冶金工业出版社, 2004.

[4]John H.Mathews, Kurtis D.Fink.数值方法[M].北京:电子工业出版社, 2002.

曲线二乘拟合篇2

关于电流互感器(current transformer,CT)铁心的磁化曲线拟合问题,国内外学者提出了多种方法,较常用的有拉格朗日插值法,最小二乘拟合法、分段线性插值法。除此之外新兴的方法还包括基于人工神经网络的曲线拟合和支持向量机,其中人工神经网络又包括BP神经网络和径向基函数神经网络,这两钟方法本质上都是利用网络的自学习功能来自动寻找最优的连接权系数,以达到非线性函数拟合的目的[1,2]。优点是选取的样本容量越大,拟合出来的精度越高,但此种方法所形成的网络结构较复杂,需要调节的连接权系数较多,训练过程也较繁琐,一些重要参数的选择不当甚至会导致整个学习过程陷入局部极值的问题;并且最后得到的拟合公式也很复杂。而传统支持向量机(SVM)虽然没有上述问题,但其训练问题是一个二次规划问题或凸规划问题,当样本数目较大时,其训练速度较慢,占用内存较大[3,4,5]。为此本文提出了一种利用最小二乘支持向量机来实现励磁特性曲线拟合的新算法。该算法不仅具有支持向量机在小样本情况下拟合精度高、泛化能力强的优点,同时还具有计算简单、求解速度快,内存需求少的特点。

1 最小二乘支持向量机原理

最小二乘支持向量机[6](Least Squares Support Vector Machines,LS-SVM),是支持向量机的一种改进,将传统支持向量机中的不等式约束改为等式约束,同时把误差平方和损失函数作为训练集的经验损失,将经验风险由偏差的一次方改为二次方,最终将求解二次规划问题转化为求解线性方程组的问题,避免了不敏感损失函数,大大降低了计算复杂度,提高了求解问题的速度和收敛精度。

1.1 具体算法描述

给定一个有N个训练样本的集合{xk,yk},k=1,2,…,N,其中训练样本n维向量xk∈R n,yk∈R。

首先用一个非线性映射φ(·)把原空间样本从Rn映射到特征空间φ(xi),这样就把低维空间的非线性逼近问题转化为高维空间的线性化逼近问题,在这个高维特征空间中构造最优决策函数:

依据结构风险最小化原则,寻找ω,b就是最小化;

根据统计学理论,函数拟合问题就变为求解如下最优化问题:

其中:|ω|2控制模型的复杂度;γ是正规化参数,控制对超出误差样本的惩罚程度;ω为权矢量;ξ为误差变量;b为偏差量;Remp为误差控制函数,也即不敏感损失函数。常用损失函数有线性损失函数,二次损失函数,Huber损失函数,当选取不同的损失函数,可构成不同形式的支持向量机。本文采用的损失函数为误差函数ξ的二次项。

用拉格朗日法求解这个优化问题:

其中:ak,k=1,2,…,N为拉格朗日乘子。根据最优化理论中的KKT(Karush-Kuhn-Tucker)条件可得到:

可推得:

其中:ai=Cξiωϕ(x i)+b+ξi-yi=0。

定义核函数ϕ(x i,yj)=ϕ(x i)iϕ(x j)是满足条件的对称函数。优化问题转化为求解以下线性方程组解的问题。

根据最小二乘法求出a与b,得到非线性拟合模型:

1.2 关于核函数的选取

LS-SVM的非线性拟合能力都是通过“核映射”的方法来实现的,对于一个具体问题,如果核参数取得不合适,LS-SVM就无法达到预期的拟合效果(数据子空间的维数决定了线性分类面能达到的最小经验误差),并且,核函数的类型[7,8]也应根据不同的情况加以正确选择才能达到理想的拟合效果。一般只要满足Mercer条件的函数都是核函数,常用的有以下几种:

如果选取式(1),那么LS-SVM实现的是一个多项式的向量机,参数q由用户自己取值;式(2),每个基函数的中心对应于一个支持向量,得到的是径向基函数向量机;式(3),实现的是一个两层的多层感知器神经网络。本文采用较常用的径向基函数作为核函数,其中:δ为核宽度。

1.3 数据预处理

为了便于后续处理,需要对训练各样本值按式(10)进行标准归一化处理

归一化后的各样本值的范围在(0,1)之间。

2 具体实现步骤

图1中是用Matlab7.0编程实现LS-SVM的结构框图,括号中是相应用到的函数,详细计算步骤如下:

(1)创建输入输出样本,用Load命令加载数据。

(2)将测得的数据样本按照式(10)进行标准归一化处理。

(3)设定最小二乘支持向量机参数gam和相应核函数sig2,type;其中:gam和sig2是最小二乘支持向量机参数,gam是正则化参数,决定适应误差的最小化和平滑程度,sig2是RBF函数的参数。Type中有两种类型:一类是关于分类的是classfication;另一类是用于函数拟合的是function approximation。

(4)用trainlssvm函数建立回归模型,根据数据样本的输入输出和上一步预先设置好的训练参数,对网络进行训练,得到最小二乘支持向量机的支持向量`alpha`和相应的阀值`b`。

(5)读取预估数据的输入,进行数据预处理,得到实际波形。

(6)将各训练后的点进行拟合,查看拟合情况。

(7)如果上步中出现的拟合波形与实际励磁特性曲线相比不能满足实际要求,可利用模型调整部分的函数进行(3)、(4)两部分建立好的回归模型的相应参数的调整,直到符合要求为止。

3 仿真波形验证

表1列出了某种硅钢材料励磁曲线B-H的部分测量数据值,共9组,将各数据依据式(10)进行标准归一化,之后将其按照以上介绍输入Matlab7.0中,经过多次对回归模型参数的调整实验,最终得到比较理想的拟合波形,如图2所示,此时Matlab7.0命令窗口中部分命令如下:

其中:B,H0为测量值;H1为训练后的样本值;%为两者之间的相对误差。

图2中,黑点为训练后的各样本点,实线为实际的励磁特性曲线。从图中可见,经调整参数训练后的各样本点基本都在实际励磁特性曲线上,拟合效果较理想。同时从表1中列出的训练后的样本值与实际值的偏差程度,发现二者之间相差甚小,拟合精度较高,可以满足实际的需要。

4 结语

本文提出使用最小二乘支持向量机来拟合电流互感器励磁特性曲线,为CT铁心饱和特性的建立提供了一种新的算法途径。实验仿真波形验证了该算法在铁心非线性逼近方面的有效性和准确性。此外,该算法除了可以应用于电流互感器的铁心拟合外,同时还可以进一步推广到其它领域的非线性曲线拟合与回归中。

摘要：针对传统支持向量机在电流互感器铁心励磁特性曲线拟合时样本数目较大出现的训练速度慢、占用内存大的问题,提出了一种新的基于最小二乘支持向量机算法。该算法将实测数据由径向基函数把非线性逼近问题转化为线性逼近问题,依据最小二乘法的思想,利用Matlab7.0求一个线性方程组的解,得到拟合曲线的近似表达式。实验结果表明,新算法训练速度快,误差小、拟合精度高。

关键词：电流互感器,最小二乘支持向量机,非线性,径向基函数,曲线拟合

参考文献

[1]施晓秋.BP网络与CMAC用于非线性曲线拟合[J].温州大学学报,2001(1):46-49.

[2]李贵存,刘万顺.用于磁化曲线拟合的高精度混合型径向基函数神经网络[J].电网技术,2001,25(12).

[3]Vapnik V N.The Nature of Statistical Learning Theory[M].New York:Springer,1995.

[4]Drucker H,Burges C J C,Kaufman L.Support vector regression machines[C].//Advances in Neural Informa-tion Processing Systems.Cambridge:MIT Press,1997:155-161.

[5]赵春晖,陈万海,郭春燕.多类支持向量机方法的研究现状与分析[J].智能系统学报,2007,2(4):11-17.ZHAO Chun-hui,CHEN Wan-hai,GUO Chun-yan.Research and analysis of methods formulticlass supportvector machines[J].CAAI Transa on Intellient Systems2007,2(4):11-17.

[6]阎威武,邵惠鹤.支持向量机和最小二乘支持向量机的比较及应用研究[J].控制与决策,2003,18(3):358-360.

[7]郑小霞,钱锋.高斯核支持向量机分类和模型参数选择研究[J].计算机工程与应用,2006,42(1):77-79.

曲线二乘拟合篇3

非线性最小二乘曲线拟合法, 就是利用非线性最小二乘法的基本思想和一些典型的非线性特征曲线来实现预测的方法。以下首先介绍线性最小二乘法的基本思想。然后, 尝试使用非线性特征曲线拟合法, 包括指数曲线拟合法、饱和指数曲线拟合法、皮尔 (Pearl) 曲线拟合法以及高柏茨 (Gompertz) 曲线拟合法。在历年中国网民的规模数据基础上, 对历年中国网民的规模做一个简单的预测。

一、最小二乘法的基本思想

作为其预测值, 对现有的x所对应的估值为:

二、指数曲线拟合法

大量的研究表明, 世界上许多事物的发展在一定时期内往往呈现出“雪崩”式的快速增长——指数增长。互联网行业同样如此, 伴随技术的进步、经济的起飞, 其增长趋势接近“指数规律”, 用指数函数建立拟合的模型较为满意。

1. 基本思想

其中, A>0, B>1。对式两边取对数, 得lgy=lg A+xlg B

令Y=lgy, a=lg A, b=lg B

这样, 便将指数曲线化为了直线, 于是, 可以套用算式 (6) 求得系数a和b, 再求反对数即获得指数增长模型中的待定系数A和B。对式 (9) 可做类似处理, 求出待定系数A和B。

2. 实例分析

截至2008年底, 中国网民规模达到2.98亿人, 较2007年增长41.9%, 互联网普及率达到22.6%, 略高于全球平均水平 (21.9%) 。继2008年6月中国网民规模超过美国, 成为全球第一之后, 中国的互联网普及再次实现飞跃, 赶上并超过了全球平均水平。 (1)

三、饱和指数曲线拟合法

1. 基本特征

实践证明, 各种事物的生长过程或成长过程均具有一定基本的共性, 主要具有下述两个特征。

(1) 发展的全过程科粗略分为三个阶段:首先是萌芽阶段, 此时, 事物处于初生或初创时期, 进展比较缓慢;其次是发展阶段, 此时在数量上出现明显的增长, 且往往有一个高速发展或突破性进展的事情;最后是饱和阶段, 此时的发展越来越慢, 趋于淘汰。 (2) 发展中存在一个上线或极限, 制约着每一个具体发展的全过程。

基于这两个事物共同的基本特征, 为了尽可能准确地拟合实际中各种各样的具体的事物发展过程, 便提出了有关描述“事物发展”的曲线拟合方法 (又称为生长模型法) , 其中, 饱和指数曲线拟合法以及皮尔 (Pearl) 曲线拟合法和高柏茨 (Gompertz) 曲线拟合法是比较具有代表性的方法。

2. 基本思想

其中, k, a, b为待求参数, 满足k>0, a>0, 0

通常用“三组法”来估计算式 (10) 中的参数, 以下简介之。

将 (2) 中的数据分别带入 (10) , 得

将方程组 (3) 中的算式顺序分成三组, 每组还有n个方程, 然后把各组加起来, 得

于是, 由方程组 (4) 便可以求得k, a, b的值:

3. 实例分析

将网民规模再次进行计算, 探讨中国网民的极限值。

按算式 (14) 计算b, a, k.。得出b=1.2868, 不在饱和指数曲线要求满足k>0, a>0, 0

四、结论

最小二乘法是一种数学优化技术, 它通过最小化误差的平方和找到一组数据的最佳函数匹配。用最简的方法求得一些不可知的真值, 而令误差平方之和为最小。在电子商务各个领域都能得到有效的应用。

参考文献

[1]陈庄, 等.ERP原理与应用教程[M].北京:电子工业出版社, ISBN7-50536-8632-8/TP.

曲线二乘拟合篇4

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面。其中一种简单的情形是:根据相关影长数据确定参照物位置, 这一技术也称为太阳影子定位技术。在本文中, 根据2015 年全国大学生数学建模竞赛A题附件1 的数据[1], 确定出该地点的位置。

2 影长与时间的函数关系

2.1 影响影长的各个参数

太阳高度角是指太阳光的入射方向和地平面之间的夹角, 可以用太阳高度角来估计影子场电影。太阳高度角h与实测地的地理纬度θ, 太阳直射角纬度 φ, 时角 Ω 满足以下函数关系式[2]:

其中, 关于各个位置参数确定如下:

太阳直射点每时都在向西移动, 每小时移过15 度经度, 在地理题的计算中可粗略取每天移动0.25 度纬度[3]。据此求解出太阳直射角纬度与日期之间的关系, 进而通过MATLAB编程, 得到了每天太阳直射点的纬度, 部分表见表1, 在4 月18 日的太阳直射点纬度为 φ=7。

一个天体的时角表示该天体是否通过了当地的子午圈, 其数值则表示了该天体与当地子午圈的角距离, 以小时来计量 (1 时角=15 度) 。根据以上理论, 建立时角的关系表达式[4], 即:

其中Cr表示北京时;LC表示经度校正 (4min/度) , 经度大于120 的经度校正为正, 小于120 的经度校正为负;EQ为时差。

2.2 建立影长随时间变化函数

根据式 (1) , 我们得到了太阳高度角与时角、太阳直射点纬度以及实测地的地理纬度的函数关系式, 进而我们可以通过太阳高度角与影子长度之间的关系, 建立影长随时间变化的函数关系式, 其步骤如下:

根据太阳高度角的定义, 得到太阳高度角与影子长度之间满足以下关系:

其中L为影子长度;H为旗杆高度3m;h为太阳高度角;

根据式 (1) , 将各个参数值代入公式进行整理, 得到太阳高度角与各个参数之间的函数关系:

再根据式 (4) , 将 (5) 转化为:

其中:

3 利用最小二乘拟合确定直杆高度以及所处位置

3.1 模型的建立

已知对应时刻的影子长度 (ti, Li) , 影长与时间的函数关系, 见公式 (6) , 利用最小二乘原理, 建立如下模型:

其中c1, c2, c3的初始值设为1, 0.1, 10

3.2 参数及经度的确定

根据参数的初始值, 将21 组 (ti, Li) 代入拟合方程中进行拟合。在拟合之前, 本文根据c1, c2, c3的表达式做了以下猜想:发现c1, c2关于其自身的自变量的表达式十分简单, 而c3表达式较为复杂, 并且c3在拟合方程中的表达式也较为复杂, 所以关于c3的拟合效果相对于c1, c2来说, 拟合效果较差, 可能得不到精确的解。

而在实际拟合过程也证明了我们的猜想。在实际拟合过程中, 随着初始值的不断变化, c1, c2变化不大或者只在固定的值的范围内小幅度变化, 而c3的变化幅度较大, 且不能找到唯一确定的值, 因此通过求出的纬度以及直杆高度来确定经度。经求解:c1=2.1316m, c2=0.313或0.1287, 即H=2.1316m, φ=18.24N或7.3969N。

将H, φ 中的数据代入 (5) 中, 用MATLAB编程计算, 我们得到21个时间点所对应的经度位置, 我们发现, 这21 个经度变化都十分小, 故我们经度该位置的经度就选取求出的21 个经度的平均值, 得到4 个经度, 进而得到四个固定直杆可能所处的地理位置: (18.24N, 105.9528E) , (18.24N, 37.4972E) , (7.3969N, 107.1988E) , (7.3969N, 36.2512E) 。

4 地理位置的检验

根据测定的关于x, y的坐标, 我们发现随着时间的增加, 影子长度也在不断增加, 说明当地时间已经过了正午时间12:00, 而此时的北京时间为14:42-15:42, 结合北京时间以东八区120 度为标准, 每向西15 度, 当地时间往前1 小时的规律, 我们计算出在79.5 E时, 北京时间为14:52, 当地时间恰好为12:00, 因此我们所确定的经度必须大于79.5E。进而确定可能的地理位置为: (18.24N, 105.9528E) 和 (7.3969N, 107.1988 E) 。

根据地图确定地理位置为 (18.24N, 105.9528E) 的地点为越南, 地理位置为 (7.3969N, 107.1988E) 的地点为泰国湾 (为靠近马来西亚的一个海峡) , 故这个位置不符合我们的标准, 故我们确定的直杆所处位置大致在越南。

5 模型的检验

我们将所确定的地理位置代入关于太阳高度角的式子 (5) 中, 画出14:42-15:42的影子长度的变化曲线l1, 同时拟合出附件1中影子长度随时间变化曲线l2, 两者之间的误差, 可以这21个时刻的误差分别为:0.0004, 0.0005, 0.0005, 0.0006, 0.0006, 0.0006, 0.0006, 0.0006, 0.0006, 0.0006, 0.0005, 0.0005, 0.0005, 0.0005, 0.0004, 0.0004, 0.0004, 0.0004, 0.0004, 0.0004, 0.0004。两条影子长度随时间变化的曲线 (见图1) 。

根据图1, 结合21 个点的误差, 可以发现两条曲线误差最大的为0.0006, 误差非常小, 认为该模型合理。

参考文献

[1]http://www.shumo.com/2015cumcm.html[OL].

[2]http://baike.baidu.com/view/86609.htm[OL].

[3]http://baike.baidu.com/link?url=_a Uu Zt IGSo3DBZgj Wi SGC2Nq H9NNBj Ec XX37JYu5MWs F6M22qio Ndju Y81e Nje P31Llpv JKTBhf_e PWF9g5CXa#1[OL].

[4]王国安, 米鸿涛, 邓天宏, 李亚男, 李兰霞.太阳高度角和日出日落时刻太阳方位角一年变化范围的计算[J].气象与环境科学, 2007, 30:161-163.

曲线二乘拟合篇5

地面三维激光扫描是获取空间数据的新技术, 可快速获取大量点的三维坐标 ( 点云数据) ,并进行三维建模。在点云数据平面拟合过程中,若只考虑观测向量的随机误差,可以构建高斯—马尔科夫 ( G-M) 模型,采用最小二乘法 ( Least Squares,LS) 获取参数最或然值; 但是由于仪器精度、观测条件、观测对象特性等各种因素的影响,扫描仪所采集的点云数据都含有随机误差,拟合过程中仅仅考虑观测向量误差是不严谨的,文献 [1] 最早提出能同时顾及观测向量误差和系数矩阵误差的总体最小二乘法 ( Total Least Squares,TLS) ; 但总体最小二乘法仅是考虑系数矩阵和观测向量均服从独立等精度分布的情况,由于采集的点云数据是不等精度的,采取经典最小二乘方法和总体最小二乘方法求解所得参数估值并非最优解; 文献 [2] 提出了加权最小二乘方法 ( Weighted Total Least Squares,WTLS) ,随后, 国内外学者对加权总体最小二乘进行了深入的理论和应用研究[3~6],并用于点云数据平面拟合[7~10]。本文主要研究了点云数据平面拟合的加权总体最小二乘法及系数矩阵协因数阵的不同确定方法,并与总体最小二乘和加权最小二乘进行比较和分析。

1加权总体最小二乘平面拟合算法

1.1加权总体最小二乘原理点云数据空间平面方程式:

式 ( 1) 中的a,b,c为待求的平面拟合参数。

将式 ( 1) 写成矩阵形式:

式 ( 2) 中,

同时顾及观测向量Z和系数矩阵A的误差,建立EIV模型:

式 ( 3) 中,Z为观测向量,ez为观测向量的随机误差矩阵,A为系数矩阵,EA为系数矩阵的随机误差矩阵; vec表示矩阵按列拉直,顺序从左到右,eA= vec( EA) ,eA∈ R3n ×1; σ02是验前单位权方差; Qe= PZ-1,QA= PA-1,分别表示观测向量的协因数阵和系数矩阵协因数阵。

加权总体最小二乘参数估计准则:

1.2几种确定系数矩阵协因数阵的方法

点云数据是不等精度观测的,应根据各点精度确定其权值。相关文献研究表明点云数据入射角越小,其点位精度越高[9,11],当入射角 θ 越小,相应的余弦值cosθ 越大,因此可以把余弦值cosθ 作为平面拟合的权重。根据点云拟合平面方程式,i点的余弦值解算公式可表达为[9]:

式 ( 6) 中,θi为i点入射角,( a,b,- 1) 为拟合平面法向量,( xi,yi,zi) 为i点的三维坐标。

则i点权值及其协因数为[9]:

1.2.1直积形式构造系数矩阵协因数阵

点云数据平面空间方程式的系数矩阵A为一固定形式,假设点云数据在x,y,z三个方向等精度观测,则有 σx= σy= σz。根据系数矩阵的特点,构造系数矩阵A的权阵:

式 ( 8) 中,P0代表系数矩阵列向量之间的权阵,对角线元素分别为1,1和0,其中1表示系数矩阵A第一列和第二列是等精度观测,第三列为常数,不需要改正,因此P0第三对角元素为0。Px∈Rn ×n代表系数矩阵行向量之间的权阵,Px∈Rn ×n代表观测向量x的权阵。PA∈ R3n ×3n代表系数矩阵A的权阵, 表示 “Kronecker积”。Px、Pz与入射角余弦值cosθ 有关。

权阵对应的协因数阵:

由式 ( 8) 和式 ( 9) 得到系数矩阵协因数阵QA:

该定权方法构造系数矩阵A的协因数阵QA,便于理解与编程,适用于系数矩阵是固定形式时使用。但当Q0奇异时,QA也奇异。系数矩阵A非固定形式时,不能使用该方法确定系数矩阵的协因数阵。

1.2.2Mahboub五条原则构造系数矩阵协因数阵

文献 [5] 提出了五条原则构造系数矩阵协因数阵QA。先把系数矩阵从左到右按列拉直,然后按照原则构造系数矩阵的协因数阵QA。根据点云数据平面拟合系数矩阵A的特点,构造的QA为:

式 ( 11) 中,Qxi= Qei= Pi-1。

该方法不受系数矩阵形式的限制,较好地保持了系数矩阵的结构性,构造出的协因数阵便于编程,而且可以直观地表现系数矩阵协因数阵的具体形式; 但构造协因数阵全过程需要手工操作,操作过程比较繁琐复杂,耗时长。

1.2.3协因数传播定律推导系数矩阵协因数阵

根据系数矩阵的特点,第三列为常数1,其协因数为0,只顾及第一第二列,按列拉直有:

式 ( 12) 中,B1∈ R2n ×2n为过渡的系数矩阵。

由协因数传播定律可把式 ( 12) 系数矩阵A1的协因数阵QA1表达为以下形式:

顾及系数矩阵A的第三列,系数矩阵A的协因数阵QA表达形式如下:

式中,

在系数矩阵没有固定形式时可采用该法推导系数矩阵协因数阵QA; 可先根据实际情况推导协因数阵Q,然后根据式( 12) 确定系数B1,顾及到系数矩阵A第三列为常数,在B1基础上确定系数B,但系数B需要复杂的手工操作,编程实现也比较麻烦。

1.3加权总体最小二乘迭代算法

点云数据拟合的Jazaeri算法具体迭代步骤如下[6]:

( 1) 将最小二乘平面拟合参数解作为起始值。根据式 ( 6) 和式 ( 7) 计算各点的协因数并组成观测向量协因数阵Qe和系数矩阵协因数阵QA;

( 2) 根据加权最小二乘法,计算平面拟合参数初始值:

( 3) 重新计算各点协因数,更新协因数阵,并计算相关值:

( 4 ) 重复步骤 ( 3 ) 直到( ε 为预设的小值) ,迭代结束。

( 5 ) 计算

( 6 ) 计算单位权中误差估计,点到拟合平面的距离( di) 以及平面拟合精度) :

式中,n为观测点个数。

2算例分析

2. 1模拟算例

根据点云数据空间平面方程式,首先给定参数真值a = c = 1,b = 2,拟合平面方程式为z = x + 2y + 1。利用MATLAB软件从该平面中随机选取1000个点 ( 见图1) ,在坐标 ( xi,yi,zi) 上均加上均值为u = 0,方差为 σ2= 1 / ( 100 × cosθi) 的随机误差( 其中i ∈ ( 1,…,1000) 表示点号,cosθi为最小二乘解所确定的i点入射角余弦值) ,得到一组含有随机误差的模拟点云数据。为避免随机误差对结果的影响, 作20次随机实验,获取20组含有随机误差的模拟点云数据。通过MATLAB编程,利用直积形式、五条原则和协因数传播律构造系数矩阵协因数阵的Jazaeri算法、文献 [2] 提出的Schaffrin算法、总体最小二乘法和加权最小二乘法对上述模拟点云数据进行参数估计,并取参数估值的平均值作为每个参数的最终计算结果,计算结果如表1所示。其中,直积形式、五条原则和协因数传播律构造系数矩阵协因数阵的Jazaeri算法分别表示WTLS1、 WTLS2和WTLS3,Schaffrin算法以WTLS4表示。

从表1可以看出,三种WTLS解算的三个参数估值都与真值最为接近,TLS次之,WLS解算的参数估值离真值最远; 三种WTLS的单位权中误差比WLS提高了60. 96% ,比TLS提高了39. 61% ; 平面拟合精度比WLS提高了9. 62% ,与TLS相比, 改善程度有限; 点到拟合平面最大距离也由WLS和TLS的12. 49cm和11. 32cm降为11. 28cm; 在迭代次数方面,三种Jazaeri算法迭代次数平均次数均为21次,Schaffrin算法平均迭代49次,在计算效率上,Jazaeri算法效果较好。

2.2实际点云数据平面拟合

使用RIEGL VZ-400三维激光扫描仪分别对地面、木板扫描,获取两组原始点云数据如图2所示。计算结果如表2、表3所示。

从表2中可以看出,三种WTLS单位权中误差比TLS提高了42. 66% ,与WLS相比改善有限,而WLS比TLS提高了42. 64% ; 由于TLS解算模型未顾及观测点点位精度,其单位权中误差最大,拟合效果最差; 而三种WTLS以各点入射角余弦值为权值,加权解算平面拟合参数,得到比较合理的平面拟合参数解; 在迭代次数上,Mahboub五条原则确定系数矩阵协因数阵的Jazaeri算法迭代次数为2次,其他三种WTLS方法迭代次数均为3次,计算效率上基本一致。

在表3中可以看出,综合三个精度评判指标, 三种WTLS相比WLS均有所改善,与TLS相比单位权中误差略大,在平面拟合精度上却有所改善; 造成这一结果的原因是由于WTLS、TLS和WLS单位权中误差计算方法不同,不能简单以单位权中误差大小判定算法的优劣,应该结合其他精度评判指标进行判定[9]; 顾及三个精度评判指标,三种WTLS与TLS精度相当,这是因为采用的木板表面附有零散建筑水泥斑点,导致斑点处的点云数据受到 “污染”,这部分点云数据入射角余弦值与实际权值有偏差,入射角余弦值已经不能代表该点的点位精度[17]; 在迭代次数上,直积形式、Mahboub五条原则以及协因数传播律确定系数矩阵协因数的Jazaeri算法和Schaffrin算法迭代次数分别为7、9、 14和23次。直积形式确定系数矩阵协因数阵的Jazaeri算法计算效率最高效,Schaffrin算法计算效率最差。

3结论

( 1) Jazaeri算法比Schaffrin算法在计算效率上要更加高效,无论是模拟算例还是在实际点云数据平面拟合上都表现其高效性;

( 2) 直积形式构造系数矩阵协因数阵的方法较简单方便、便于理解以及编程且计算效率快速高效。Mahboub五条原则构造系数矩阵协因数阵在计算效率上相对高效,但在构造协因数阵过程需要手工操作,大大增加了工作量和复杂程度。点云平面拟合以直积形式构造系数矩阵协因数阵较为理想;

曲线二乘拟合篇6

关键词：数字图像相关,移动最小二乘法,应变测量

数字图像相关技术是随着光电技术,视频技术,计算机技术和图象处理技术的不断成熟发展起来的一种非接触式光测力学方法。这种方法最早在20世纪80年代初就被提出并不断发展完善,在实验力学中应用越来越广泛。它通过对物体表面散斑在变形前后的图像进行相关运算来确定物体的位移,再用得到的位移场计算得出应变场。因此精确的位移场计算是得到精确应变场的前提,所以提高测量精度的亚像素位移定位算法被认为是数字相关方法的关键技术之一。为此,Bruck等[1]提出了Newton-Raphson(NR)算法,应用双三次样条插值对亚像素进行重建。Lu等[2]在NR算法中考虑了模板窗口的非均匀变形,引进了二次位移梯度对计算相关系数的影响。芮嘉白等[3]提出了十字搜索算法,将整像素和亚像素定位的二维搜索算法转变成一维搜索,适用于相似系数呈单峰性好的情况,多用于像素级的搜索定位。李善祥等[4]提出了相似系数的曲面拟合算法,具有抗噪声能力强、计算效率高等特点。Peng Zhou等[5]将梯度算法成功应用于实验力学中,该方法易于实现,求解效率高。曲面拟合法和梯度算法不需要对亚像素进行重建,直接进行亚像素位移定位,算法精度高,效率高,在数字图像相关方法得到很好的应用。

潘兵等[6,7]对亚像素位移定位算法的进展做了详细的论述,并对几种主要的亚像素位移测量算法进行了研究对比。文献[6]对主要的几种亚像素定位算法的研究对比表明,NR算法精度较梯度算法高,曲面拟合算法精度最低,且NR算法能同时得到位移应变信息。但NR算法应用了双三次样条插值对亚像素进行重建,计算量大,且在得到应变时受插值误差和图像局部灰度值的影响很大。梯度算法和曲面拟合算法进行亚像素位移定位时误差比NR算法大,但可以经过对该算法得到的位移场进行拟合平滑来减小误差,进而得到比较可靠的应变场。本文尝试应用移动最小二乘法来对梯度算法得到的位移场进行拟合,应用拟合后的位移场计算得到应变场。编制了相应的计算程序,得到了可靠的位移场和合理的应变场。

1 基本原理

数字散斑图像相关技术是根据物体表面的随机散斑粒子在变形前后的概率相关性来确定物体的变形的。数字散斑图像相关技术方法的基本过程是:利用物体表面制作的散斑场或自然散斑场,由CCD采集物体变形前后物体表面的散斑图,根据两幅散斑图像灰度值进行处理,实现表面位移场的测量。这个过程的基本原理是在物体变形前的散斑场中选定一个子集作为参考子集区域,称之为样本子区,然后在物体变形后的散斑场中搜索与参考子集相对应的子集,这个子集称之为目标子区。由概率与统计理论,两个子集相关系数定义为

$\begin{array}{l} C (u, v, \frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}, \frac{\partial v}{\partial x}, \frac{\partial v}{\partial y}) = \\ \frac{\sum \sum [f (x_{i}, y_{j}) - \bar{f}] [g (x_{i}^{*}, y_{j}^{*}) - \bar{g}]}{\sqrt{\sum \sum [f (x_{i}, y_{j}) - \bar{f}]^{2} [g (x_{i}^{*}, y_{j}^{*}) - \bar{g}]^{2}}} (1) \end{array}$

式(1)中,C为相关系数;f(xi,yj)和g(x*i,y*j)分别表示变形前后在(xi,yj)和(x*i,y*j)点的灰度值。 $\bar{f}$ 和 $\bar{g}$ 表示样本子区和目标子区的灰度平均值。

通过试凑位移法,在变形后的图像上移动目标子区,使样本子区和目标子区的相关系数C取得最大值,便可得到物体变形前后的位移值,进而利用差分运算得到应变值。

2 搜索定位及算法

像素级的搜索定位采用改进后的十字搜索算法。具体算法步骤是:首先定位第一个点时对其位移初值u, v赋零直接进行十字搜索,定位后记下第一个点像素级的位移值,然后把它作为第二个点的位移初值进行十字搜索,以后始终把前一点得到的像素级的位移作为后一点的初始位移进行十字搜索,这样就消除了整个视场的刚体位移,提高了效率。亚像素级位移定位尝试应用了梯度算法。

3 位移场拟合分析

应变是位移的导数,因而数字散斑图像相关技术得到的位移场微小误差将导致应变场的较大偏差。故本文提出在DIC计算得到的位移场后,对位移场应用移动最小二乘法来进行拟合,进而对其进行平滑和去噪,从而得到较可靠的应变场。移动最小二乘法[8,9]的拟合函数可以表示为

f(x,y)=pT(x,y)α(x,y) (2)

式(2)中α(x,y)为待求系数,它是坐标的函数;p(x,y)是基函数,现采用二次基函数,即p(x,y)=[1,x,y,x2,xy,y2]T。α(x,y)的取值应使加权离散范式

$J = \sum_{Ι = 1}^{n} w (x - x_{Ι}, y - y_{Ι}) [f (x, y) - f_{Ι} (x_{Ι}, y_{Ι})]^{2} (3)$

取得最小值。所以,令式(3)对α(x,y)的偏导数为零,得故α(x,y)=A-1(x,y)B(x,y)F。则拟合函数表达式

f(x,y)=pT(x,y)A-1(x,y)B(x,y)F (4)

式(4)中:fI(xI,yI)是样点(xI, yI)处的测量值,w(x-xI,y-yI)是样点(xI, yI)处的权函数。

B(x,y)=[w(x-x1,y-y1)p(x1,y1),

w(x-x2,y-y2)p(x2,y2),…,

w(x-xn,y-yn)p(xn,yn)]。

为了得到平滑性好的曲面,采用三次样条权函数,即:

$\begin{array}{l} w (x - x_{Ι}, y - y_{Ι}) = \\ {\begin{matrix} \frac{2}{3} - 4 \bar{s}^{2} + 4 \bar{s}^{3} ‚ & (\bar{s} \leq \frac{1}{2}) \\ \frac{4}{3} - 4 \bar{s} + 4 \bar{s}^{2} - \frac{4}{3} \bar{s}^{3} ‚ & (\frac{1}{2} \leq \bar{s} \leq 1) \\ 0 ‚ & (\bar{s} > 1) \end{matrix} (5) \end{array}$

式(5)中, $\bar{s} = \frac{s}{s_{\max}} ‚ s = \sqrt{(x - x_{Ι})^{2} + (y - y_{Ι})^{2}}$ 。

移动最小二乘法拟合不需要事先确定拟合函数的类型,不需要求解线性方程组,避免了求解方程组时系数矩阵的病态情况,且具有精度高的特点。利用式(4)拟合得到的是去噪和平滑后的位移场,可以提高精确度,减小误差。且拟合的过程,把离散的位移场转变成了全场性的连续位移场。

4 算例

4.1 数值模拟实验

数值模拟生成仿真散斑图是验证和评价数字散斑图像相关技术性能的重要方式,它可以精确控制散斑图的位移和应变信息,还能排除实际图片可能受到的各种噪声干扰。用计算机分别生成了一对既有刚体位移又有均匀变形的图像(如图1;图像大小:256 pixels×256 pixels),位移模式为:u(x)=3.01-0.03(x-20); v(y)=3.01+0.03y;其中x、y的取值范围为:[60,200],用编制的程序计算得到的位移场和应变场如图2、图3。

两张图象行列每间隔10像素取一样点进行搜索定位,然后进行移动最小二乘拟合样点位移值。利用拟合后的位移场进行差分运算得到应变场。如图2和图3。由图2和图3中的(a)、(b)图可以看出拟合后的位移场平滑性好,与预先施加的随轴向坐标线性变化的位移场分布十分吻合。如图2和图3中的(c)图所示,均匀应变的应变场在整个区域是一个常量,从(c)图中可以看出,本文提出对位移场进行移动最小二乘法拟合后再进行应变计算,得到的应变是可靠和精确的。

4.2 复合材料平板拉伸试验

某复合材料板,长300 mm,宽36 mm,厚2.4 mm,铺层顺序为[45/0/-45/90/0/0/45/0/-45/-45]s。用数字散斑图像相关技术来测量它的模量。实验过程:在复合材料平板表面上用黑白漆交替喷涂,直到板表面出现具有明显反差的随机散斑分布,然后在MTS 810 Teststar ±100 KN液压伺服材料试验系统上进行拉伸实验,在实验开始用相机(SONY a200 、1 020万像素;Tokina AT—X (50—250)mm镜头,视场约为20 mm×15 mm)记录下初始散斑图片作为参考图片,然后在两端(沿x方向)每增加3 000 N拉力记录一次变形信息。然后利用编制的程序进行应变计算(应变场的平均值作为复合材料平板的均匀应变值)。有两个试件,标号分别是D26和D24,用数字散斑图像相关技术法来进行计算得到试件D26的表面应变数据,并用最小二乘法进行拟合。D24数据是用传统的应变电测方法测得的数据。

从图5两种测量方法对复合材料板D26和D24的测量数据的对比可以看出,测量得到的数据是吻合的,从而证实编制的程序在实验力学中应用是可以达到实验要求的。

图5(a)是x方向应变与应力的关系图。拟合曲线的斜率就是杨氏模量Ex。故可以得出复合材料板的杨氏模量Ex=61.1 GPa。图5(b)是y方向应变和应力的关系图。拟合曲线的斜率就是杨氏模量Exy,故可以得出复合材料板的Exy=133.3 GPa,得到板的泊松系数 $ν_{x y} = \frac{E_{x}}{E_{x y}}$ =0.458 4。

5 结论

(1) 对十字搜索法根据具体的实验做了修改,简化了算法步骤,效率更高。

(2) 对亚像素定位计算的位移场应用最小移动二乘法进行拟合,实践证明,拟合后的位移能有效减小误差。

(3) 应用拟合后的位移场计算得到了可靠的应变场,使得数字散斑图像相关技术能够应用于实验力学中,并用数字散斑图像相关技术测量了某复合材料板的拉伸模量。

参考文献

[1] Bruck HA,McNeil S R,Sutton MA,et al.Digital image corretion:using Newton-Rapshon method of partial differential correction.Ex-perimental Mechanics,1989;29(3):261—267

[2] Lu H,Cary P D.Deformation measurement by digital image correla-tion:implementation of a second-order displacement gradient.Experi-ment Mechanics,2000,40(4):393—400

[3]芮嘉白,金观昌,徐秉业.一种新的数字散斑相关方法及其应用.力学学报,1994;(09):599—607

[4]李善祥,孙一翎,李景镇,数字散斑相关测量中亚像素位移搜索的曲面拟合研究.光子学报,1999;(07):638—640

[5] Peng Zhou,Letnnet E.Goodson subpixel displacement and deforma-tion gradient measurement using digital image/speckle correlation.Opt Eng,2001;40(8):1613—1620

[6]潘兵,谢惠民,戴福隆.数字图像相关中亚像素位移测量算法的研究.力学学报,2007;36(2):245—252

[7]潘兵,谢惠民,续伯钦,等.数字图像相关中亚像素位移定位算法进展.力学进展,2005;(08):345—351

[8]曾清红,卢德唐.基于移动最小二乘法的曲线曲面拟合.工程图学学报,2004;(1):84—89

曲线二乘拟合篇7

空间直线拟合在工程中具有重要应用[1,2,3,4,5,6,7,8],如在计算机视觉[4]、三坐标测量问题[5]、机器人焊缝中心三维定位[6]和空间管线网络干线的确定[7]中。

由于最小二乘准则下问题求解比较容易,一些文献在该准则下给出确定空间直线的算法。文献[1]将直线理解为两个平面的交集,通过求取两个平面方程,从而确定直线。文献 [2,3]根据最佳平方逼近原理,证明直线必过数据中心,并用逐步迭代给定方向向量,从而确定点向式直线方程。文献 [5]则通过两次降维,先按全最小二乘准则拟合二维平面,将空间中的点投影到该平面上,再将该平面上数据利用全最小二乘准则拟合直线,从而确定空间直线。文献[6]先确定直线经过的点,再求导得到一个非线性方程组,解方程组确定方向向量,从而最终确定空间直线。这些算法比较复杂,本文给出一种基于主成分分析的简便算法,通过计算机仿真验证了算法的有效性。

由于最小一乘准则下确定空间直线是一个不可微优化问题,难度较大。文献 [7]提出了空间中按全加权最小一乘准则下确定直线问题,但没有给出具体的算法,文献[8,9]给出了平面上的求解方法,没有对空间中全加权最小一乘准则下确定直线问题进行讨论。本文给出一种基于无约束优化的搜索算法,并进行了计算机仿真。

1 空间直线拟合问题

空间直线拟合问题可以描述为:在三维空间给定n个点P1(x1,y1,z1),P2(x2,y2,z2),…,Pn(xn,yn,zn),求空间直线L: $\frac{x - x_{0}}{m} = \frac{y - y_{0}}{n} = \frac{z - z_{0}}{p}$ ,其中M(x,y,z)为直线上任意点,M0(x0,y0,z0)为直线上确定的一点,非零向量 $\vec{s} = (m, n, p)$ 为直线的方向向量,使得目标函数:

$\sum_{i = 1}^{n} d^{2} (Ρ_{i}, L) (1)$

为最小,其中 $d (Ρ_{i}, L) = \frac{| \vec{Μ_{0} Ρ_{i}} \times \vec{s} |}{| \vec{s} |}$ 表示点Pi到直线L的距离[10] ,其中|·|表示向量的模。式(1)实际上是按全最小二乘拟合空间直线,等价于文献[1,2,3,4,5,6]中拟合空间直线的准则。

如在上述假设下,改为求目标函数:

$\sum_{i = 1}^{n} d (Ρ_{i}, L) (2)$

最小,则是按全最小一乘拟合空间直线,这样得到的直线不易受奇异点(outliers)的影响。

如在上述假设下,改为求目标函数:

$\sum_{i = 1}^{n} w_{i} d (Ρ_{i}, L) (3)$

最小,其中wi(i=1,2,…,n)为权重,则是按加权全最小一乘拟合空间直线,这样得到的直线可以用于空间中管道的主干线确定问题。当wi=1(i=1,2,…,n)时,式(3)变为式(2)。

2 拟合空间直线的算法

2.1 全最小二乘准则下的拟合算法

将空间中三维分布的点集拟合为一维的直线,其本质就是降维,因此可以按照主成分的思想,对三维分布的空间数据进行主成分分析,提取第一主成分向量作为直线方向向量。文献[3]证明在最小二乘准则下,空间直线必通过数据中心点,因此可以由点向式确定空间直线方程。按此思想给出求问题式(1)中直线的算法如下:

第1步输入三维空间中点坐标的n个数据点矩阵Xn×3,求出均值 $\bar{X} = (\bar{x}, \bar{y}, \bar{z})$ ,拟合直线必过数据中心 $(\bar{x}, \bar{y}, \bar{z})$ ;

第2步求协方差矩阵 $S = \frac{1}{n} (X - \vec{1}_{n} \bar{X})^{Τ} \cdot (X - \vec{1}_{n} \bar{X})$ ,其中1n为n行1列元素全为1的n维列向量;

第3步求出S的最大特征值λ1和对应特征向量 $\vec{v}_{1}$ ,以 $\vec{s} = \vec{v}_{1} = (a, b, c)$ 作为直线方向向量;

第4步由点向式确定直线 $\frac{x - \bar{x}}{a} = \frac{y - \bar{y}}{b} = \frac{z - \bar{z}}{c}$ ,其中(x,y,z)为直线上任意一点。

由主成分意义,点到直线距离平方和为n(λ2+λ3),其中λ2与λ3为S的其余两个特征值,与由式(1)得到的结果相同。

2.2 加权全最小一乘准则下的拟合算法

由于问题式(2)是式(3)的特殊情形,因此讨论式(3)。在问题式(3)中目标函数中含有绝对值,是一个不可微函数,因而求解比较困难。

要在空间中确定一条直线,无论采用确定两个相交平面还是点向式方程都至少需要确定6个未知参数。

将问题式(3)中的目标函数看作是一个含有6个自变量的不可微无约束优化问题,可以采用文献[11]中不依赖于导数的Nelder-Mead的单纯形算法与Hooke-Jeeves的模式搜索法求解。MATLAB中有基于单纯形算法的fminsearch()命令可以直接调用。

本文采用点向式方程确定直线,约定6个自变量按下面方式排列 $\vec{v} = (x_{0}, y_{0}, z_{0}, a, b, c)$ ,给出求问题式(3)中直线的算法如下:

第1步输入三维空间中点坐标的n个数据点矩阵Xn×3,给出初始值,给出 $\vec{v} = (x_{0}, y_{0}, z_{0}, a, b, c)$ 的初始值。初始值可以选为全最小二乘准则下给出的值。

第2步利用初始值和MATLAB中fminsearch()对目标函数 $\sum_{i = 1}^{n} w_{i} d (Ρ_{i}, L)$ 进行迭代求解,给出收敛值,从而确定直线经过的点(x0,y0,z0)和方向向量(a,b,c)。收敛的最终结果(x0,y0,z0)不唯一,方向向量(a,b,c)经标准化并指定方向后唯一。

第3步由点向式确定直线 $\frac{x - x_{0}}{a} = \frac{y - y_{0}}{b} = \frac{z - z_{0}}{c}$ ,其中(x,y,z)为直线上任意一点。

2.3 空间点到拟合直线的投影

为了在计算机仿真中绘制空间中点到拟合直线的垂线,这里给出点到直线的垂足坐标确定方法。

由于空间拟合直线已经求出,故可以确定空间中每一个点Pi(xi,yi,zi)在直线上的垂直投影点P′i(x′i,y′i,z′i)。由空间解析几何知识[10]可知投影点P′i(x′i,y′i,z′i)为直线在经过点Pi(xi,yi,zi)且与直线垂直平面上的垂足,联立直线的参数式方程与直线垂直的平面方程。

${\begin{cases} x = \bar{x} + a t y = \bar{y} + b t z = \bar{z} + c t \\ a (x - x_{i}) + b (y - y_{i}) + c (z - z_{i}) = 0 \end{cases} (4)$

可以解出 $t = \frac{a (x_{i} - \bar{x}) + b (y_{i} - \bar{y}) + c (z_{i} - \bar{z})}{a^{2} + b^{2} + c^{2}}$ ,从而代入参数式方程解出(x,y,z),即为垂足P′i(x′i,y′i,z′i)。

3 计算机仿真

为了检验本文算法的正确性和有效性,给出下面例子。

仿真环境为 Intel酷睿2双核 P8400、2G内存、Windows Vista OS、MATLAB7.5环境。

例1 设直线方程为x=1+t,y=1+t,z=1+t,t为参数。

参数t取值区间为[-6,6],先将该区间等分为30个小区间,t取区间端点,从而可在直线上确定31个点,可以确定31×3的数据矩阵。对直线上每个点进行随机扰动,让x,y,z分别加上随机数种子为rand(‘twister’,10),rand(‘twister’,20),rand(‘twister’,30) 的服从[-1,1]上均匀分布的n维随机列向量。绘出这些点三维散点图,如图1所示。

表1给出不同准则和算法得到的数值结果。其中方向向量已经标准化并约定向上(c>0),加权最小一乘准则中取权重全为1,即最小一乘,优化使用MATLAB自带的fminsearch()。

按最小二乘准则,绘制出拟合直线和点到拟合直线的垂线,如图2所示。

按最小一乘准则绘制出拟合直线和点到拟合直线的垂线,如图3所示。

4 结论

本文讨论了空间中的直线拟合问题,针对全最小二乘和全加权最小一乘准则下,分别给出了基于主成分分析和无约束优化的空间直线拟合方法,这些方法在计算机上实现起来比较简便,计算机仿真验证了算法的有效性。这对使用空间直线拟合方法的科研人员具有一定的参考价值。

空间中曲线拟合比较复杂,未见到统一的方法,然而特殊情况比较容易。如空间中拟合圆周,可以先拟合平面,将三维数据投影到拟合平面上(降维到2),再根据平面法向量的方向角对平面上投影数据旋转使得一个分量全部为0,将旋转后的数据拟合圆周,得到圆周对应的圆心和半径。最后进行旋转逆变换,将球面方程与平面方程联立即得空间圆周方程。

参考文献

[1]袭杨.空间直线拟合的一种方法[J].齐齐哈尔大学学报,2009,25(2):64-68.

[2]许海涛,高彩霞.直线拟合算法[J].电脑知识与技术,2009,5(4):864-867.

[3]杜明芳.空间直线拟合[J].北京印刷学院学报,1996,4(2):27-31.

[4]李畅,张剑清,邹松柏,等.自动识别街道立面凹凸边界[J].计算机工程与应用,2009(15):6-8.

[5]胡学军,王俊亮,王刚.全最小二乘法在三坐标测量中的应用[J].武汉工程大学学报,2008,30(4):112-113.

[6]郭学军,王喜平,刘贺平.机器人焊缝中心三维定位的解析方法[J].数学的实践与认识,2009,39(9):123-127.

[7]梁怡,吴可法.一类非光滑最优化问题的有限步解法[C].中国运筹学会第六届学术交流会论文集(章祥荪等主编),Global-Link出版公司,香港,2000,801-811.

[8]林诒勋,尚松蒲.平面上的点—线选址问题[J].运筹学学报,2002,6(3):61-68.

[9]冯守平,杨桂元.关于加权全最小一乘法[J].应用概率统计,2009,25(2):135-142.

[10]同济大学应用数学系.高等数学(第五版上册)[M].北京:高等教育出版社,2002.

【曲线二乘拟合】推荐阅读：