二乘拟合

2024-08-04

二乘拟合（共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

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面。其中一种简单的情形是:根据相关影长数据确定参照物位置, 这一技术也称为太阳影子定位技术。在本文中, 根据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.

二乘拟合篇3

地面三维激光扫描是获取空间数据的新技术, 可快速获取大量点的三维坐标 ( 点云数据) ,并进行三维建模。在点云数据平面拟合过程中,若只考虑观测向量的随机误差,可以构建高斯—马尔科夫 ( 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五条原则构造系数矩阵协因数阵在计算效率上相对高效,但在构造协因数阵过程需要手工操作,大大增加了工作量和复杂程度。点云平面拟合以直积形式构造系数矩阵协因数阵较为理想;

二乘拟合篇4

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

数字图像相关技术是随着光电技术,视频技术,计算机技术和图象处理技术的不断成熟发展起来的一种非接触式光测力学方法。这种方法最早在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

二乘拟合篇5

利用GPS定位技术测定的GPS高程是基于WGS-84参考椭球的大地高,而工程中所采用的高程一般是基于似大地水准面的正常高,它们两者之间相差一个高程异常,因此如何更好地确定高程异常是确定正常高的关键。GPS高程拟合方法主要有地球重力场模型法和数学模型拟合法。地球重力场模型法需要精度高、分布良好的重力数据和地形数据,这在实际工程中往往很难满足。常采用的数学模型拟合法有多项式拟合法、多面函数法、平面函数拟合法、二次曲面拟合法等方法[1]。由于这些拟合法对似大地水准面作了人为的假设,存在一定的模型误差。近年来一些学者提出了基于神经网络的GPS高程拟合方法[2~4]、支持向量机法[5~7]等,神经网络算法虽然简单,可塑等特点,但它存在收敛速度慢、过学习,容易陷入局部最优、有时找不到全局最小值等缺点。Vapnik[8]在统计学理论结构的风险最小化原理上首次提出了支持向量机(Suppor Vector Machine,SVM)。该方法不仅能较好地解决了小样本、非线性、高维数、局部极小点等实际问题,还克服了如神经网络等一般机器学习中存在的过学习、局部优化和样本数量要求多等问题。LS-SVM是支持向量机的一种改进,比SVM具有更快的求解速度,所需计算资源少等优点。

本文将最小二乘支持向量机(LS-SVM)引入GPS高程拟合领域中,通过实例数据进行建模分析,并与传统的多项式拟合、BP神经网络、GA-BP神经网络进行比较,结果表明LS-SVM方法拟合的精度较高。

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

最小二乘支持向量机是Suykens和Vandewalb在1999年提出的[9],其原理是将最小二乘估计引入支持向量机中。它是将传统支持向量机中的不等式约束改为等式约束,且将误差平方和(Sum Squares Error)损失函数作为训练集的经验损失,这样就把解二次规划问题转化为求解线性方程组求解问题,提高求解问题的速度和收敛精度,同时使误差平方项达到最小化的计算过程。因此,LS-SVM能够解决SVM在样本数目多时运算速度慢、抗噪能力差、复杂等问题。给定一个有M个训练样本的集合(xi,yi),i=1,2…M,其中训练样本m维向量,xi∈Rm,输出数据yi∈R。按照结构风险最小化原则,函数拟合问题可转换为下列函数约束优化问题[7,9]:

约束条件为:

式中:W∈Rm为权矢量;g(xi)是将x从输入空间映射到特征空间的函数;ξi∈R为i个样本点的训练误差;C为正则化参数,b为偏置量。

同时定义拉格朗日函数,将式(2)转化为无约束目标函数:

式(3)中:a=(a1,…aM)为拉格朗日乘子,ξ=(ξ1,…ξM),i=1,…M。

分别对(3)式中的W、b、ξ求导,并令其结果为0,同时结合式(2)中的约束条件,则得下列优化条件:

式中:i=1,…M。

然后,消元去掉ξi和W,得到式(8):

式(8)中:Ω=[g(x1)g(x2)…g(xM)]T,珗a=(a1,…aM)T,Y=[y1y2…ym]T,l为M×1单位列向量。非线性函数拟合的ΩΩT内积运算,可用满足Mercer条件的任何对称的核函数对应于特征空间的点积的核函数K(xk,xj)代替。这样,可以得到LS-SVM的非线性预测模型:

目前常用的核函数有以下三类,多项式核:K(x,x)=(x·y+1)P;高斯核(径向基函数核):K(x,x)=exp(-x-y2/(2σ2))(σ为核函数宽度);线性核:K(x,x)=x·y。由于数值限制条件和参数较少以及优秀的局部逼近特性,通常径向基核函数作为LS-SVM核函数的首选,因此本文采用径向基核函数进行分析计算。

由上述原理可知,LS-SVM仅需确定核函数的宽度σ和正则化参数C,而不需要不敏感损失函数的ε值,这既简化了计算,又便于实际应用。

2 GPS高程拟合精度评定

对拟合后的结果要进行精度评定,评价的指标较多,有方差、标准差等,鉴于样本中既有拟合点又有检核点,本文对高程拟合精度采用内符合精度指标μ1和外符合精度指标μ2进行评价:

其中:v1为拟合高程异常的残差,n1为拟合计算点的个数,v2为拟合外推高程异常的残差,n2为拟合外推计算点的个数。内、外符合精度越小,表明拟合和预测的精度越高;反之,则表明拟合和预测的精度越差。

3 实例分析

某沿江地形平缓区域的GPS控制网有无粗差且同精度的水准高程点17个,平均边长约1km,区域面积约10km2,按国家GPS网B级要求实测,采用二等水准联测各GPS点,即17点每个点都获取了平面位置和高程异常[10],选择其中的8个平均分布的点作为学习样本点,其余的9个点作为测试样本点。

采用LS-SVM模型进行GPS高程拟合实验计算过程:

(1)归一化处理。在进行数据拟合之前,由于所用数据数值较大,为了避免其对拟合结果产生影响,需要对数据进行归一化处理。

(2)采用基于交叉核实(cross validate)的格网搜寻法,并从大范围至小范围进行逐步搜索,确定LS-SVM参数C和σ。

(3)初始化模型并进行训练学习和测试,其中分别采用了tunelssvm,initlssvm,trainlssvm,simlssvm四个函数。

随后,分别采用多项式拟合、BP神经网络、GA-BP神经网络方法进行GPS高程拟合,将拟合的结果与LS-SVM模型拟合结果进行比较,见表1和图1。

从表1可看出,支持向量机的内、外符合精度分别为0.71cm和0.99cm,均为最小,其拟合精度最高,多项式的拟合的精度最差,GA-BP神经网络比BP神经网络的拟合精度稍高。从图1可以看出,LS-SVM模型拟合的残差变化比较平稳,多项式拟合、BP神经网络和GA-BP神经网络方法三种拟合方法拟合的残差变化上下起伏较大。

4 小结

本文探讨了LS-SVM在GPS高程拟合的运用,然后通过实例说明,在恰当选取核函数和核参数的前提下,LS-SVM在GPS高程拟合中的优越性,拟合精度比多项式拟合、BP神经网络和GA-BP神经网络方法都要高,为GPS高程拟合提供了一种新途径。

摘要：本文论述了最小二乘支持向量机(LS-SVM)的算法,提出了基于最小二乘支持向量机进行GPS高程拟合的方法,并在MATLAB中编制了相应的LS-SVM程序,建立了相应的GPS高程拟合模型。以实例数据讨论了LS-SVM的GPS高程拟合的分析方法,通过与多项式拟合、BP神经网络拟合、GA-BP神经网络拟合的结果比较,可知LS-SVM的拟合精度较高。

关键词：GPS,LS-SVM,神经网络,高程拟合

参考文献

[1]乔仰文等.GPS高程转换的若干问题的研究[J].测绘通报,1999,(11):17～19.

[2]吴良才,胡振琪.基于神经网络的GPS高程转换方法[J].工程勘察,2004,(2):49～51.

[3]尹爱明,张楚.转换GPS高程的BP神经网络方法研究[J].测绘科学,2008,33(6):78～80.

[4]孔令杰,李捷斌,陈伟.基于遗传算法的BP神经网络在高程拟合中的应用[J].工程勘察,2010,(3):61～64.

[5]谢波,刘连旺.支持向量机在GPS高程异常中的应用[J].测绘科学,2011,36(1):172～174.

[6]周理含.最小二乘支持向量机在GPS高程转换中的应用[J].工程地球物理学报,2010,7(2):243～247.

[7]黄磊等.粒子群最小二乘支持向量机在GPS高程拟合中的应用[J].测绘科学,2010,35(5):172～174.

[8]Vapnik.The Nature of Statistical Learning Theory[M].NewYork:Springer-Verlag,1995.

[9]方瑞明.支持向量机理论及其应用分析[M].北京:中国电力出版社,2007.

二乘拟合篇6

影响皮带秤称重系统精度要素之一就是张力, 利用张力补偿可以提高称重长期稳定性和精度。张力越准确, 称重就越精确。因此, 准确测得张力值十分重要。另外, 张力不稳定是由机械的不平衡与磨损和电机响应运转等方面造成的[1]。本文通过检测皮带垂度并对垂度数据进行算法处理, 从而对张力的不稳定因素进行补偿, 最终根据研究得到的垂度与张力数学模型来计算出张力。

本文先理论推导出张力随流量变化的数学模型, 其次, 通过离散垂度测量来建立垂度与流量的数据表, 利用最小二乘曲线拟合法对垂度与流量曲线拟合, 实现对数据的综合分析, 为研究物料流量与皮带垂度的关系提供科学的参考依据[2]。然后对拟合结果进行分析, 最终得到垂度、皮带倾角等参数的张力数学模型。

1 现场输送带载料模型分析

现场的输送带安装示意图如图1所示, 在具体研究中简化为图2所示模型。

其中, O1为托辊A, B之间输送带上的垂点, AB倾角为β1, 其中虚线是输送带传输方向, 实心直线为水平面, 根据实际测试, O1垂点对应的输送带挠度即垂度v<3%L。现取一小段输送带O1D分析, 建立xo1 y坐标系, 如图3所示。

取O1为坐标原点, 且以过O1点的切线为X轴, 与X轴垂直作Y轴, 其中, T1, TD分别为O1和D点处的张力, P为均布载荷;D点坐标为 (x, y) , 载荷重心坐标E (xe, ye) , 因输送带挠度较小, 曲线O1D近似等于x。因此, O1D上负荷可视为x P, 且, 同时, 根据条件可知。使用平衡条件可得:

式中, MD为关于D点的力矩, P为皮带及物料的质量分布简称流量, β1地面与水平面间的夹角。经整理得: (2)

又已知O1为AB垂点, 且为中间, 因此可得x=L/2, 代入上式得:

式中, L为两托辊间的距离。其中, L与β1均为已知量。当P与y两个变量检测出来, 对应O1点的张力也就随之求出。

2 基于最小二乘法的垂度-流量曲线拟合

2.1 最小二乘法曲线拟合原理

假设流量fi与垂度xi, 其中i=0, 1, 2, …, n。拟合关系为n次多项式, 函数为, 其中a0, a1, …, an为下面方程组的解:

如果方程组中的系数矩阵为非奇异, 则最优解唯一, 这里, 是Rm+1上向量的加权内积, 即, , 通过该方程使用Matlab软件求解出系数a0, a1, …, an, 就可以解得n次多项式[3]。

由于皮带传输物料的重量不均匀, 只有通过垂度值来反推出流量, 从而进一步估算载物的重量, 辅助称重系统判定称重单元计算结果的正确性与故障源。因此, 研究物料流量固定, 通过大量数据采样得到垂度值, 经均值数据处理, 得到一组垂度与流量的数据, 重复该方法获取多组不同流量下检测到的垂度值, 来准确反映流量与垂度的关系。

2.2 流量与垂度实验

实验条件:1) 标定初始值104mm;2) 输送带预紧力为840kg;3) 皮带宽为1.2m;4) 料斗下料给定值设为100t/h。

本实验通过改变电机驱动输送带运转的频率来改变落在输送带的物料重量即流量, 并检测获取相应的垂度, 得到表1所示的9组流量范围在14.62~46.3kg/m内的垂度与流量的关系数据。

2.3 垂度-流量特性研究

根据表1数据采用Matlab对拟合模型y=a0x0+a1x+…+anxn进行多项式拟合, 由最小二乘法确定多项式的系数a0, a1, …, an, 通过拟合发现高于5次多项式曲线就会严重失真, 因此只考虑5次以内拟合曲线, 通过Matlab计算获得各多项式曲线对应的残差如表2所示。

由表2可以观察到5次多项式拟合曲线残差为0.80778, 拟合效果最好, 其拟合曲线如图4所示。

由图可以观察出数据与曲线拟合度非常好, 因此选取5次多项式拟合曲线为最优拟合曲线, 其5次多项式空间为, 且曲线方程[4]如下:

另外, 根据皮带本身性能和受力情况分析可以设想曲线y=f (x) 是双曲型, 它与实际数据规律大致符合, 可以通过变量代换, 化为线性参数的数学模型。拟合数据。其中 (i=1, 2, …, 9) 。其中ti, zi分别由原始数据xi, yi据变量代换公式计算出来, 转换后的数据表如表3所示。

根据最小二乘法曲线拟合基本原理, 并且简化求解方程组得

将以上计算值代入简化方程组得如下方程组:

计算得a1=0.33822, a0=-0.0078494从而拟合曲线为

3 最小二乘法曲线拟合结果分析

3.1 拟合模型误差比较

根据实测垂度数据, 分别利用5次和双曲型拟合曲线计算出输送带上的物料流量, 进而计算出理论流量与计算流量之间的绝对误差, 误差平方[5]。具体计算结果如表4所示。

由表4算得5次与双曲型拟合曲线误差平方总和分别为0.6683, 1.2104。5次多项式拟合曲线误差平方更小, 拟合效果更好。

由上表分析研究选择5次多项式拟合曲线作为最优曲线拟合结果, 其函数如下:

3.2 张力数学模型

根据5次多项式拟合曲线测得的流量值, 将其带入张力公式 (其中β1=6°, L=1.2m) , 张力公式简化为, P=y+10 (其中y为检测流量, 10为皮带重量为一常量其单位为kg/m) 。再将垂度数据代入上式就可算得张力值, 如表5所示。

从表5中可以观察到垂度为4.5~10.3mm范围内张力随着垂度变大而减小, 张力的极小值产生于垂度为10.3mm~11.5mm之间, 张力自极小值之后会随着垂度的增加而增加, 由以上规律可判定张力与皮带性能有很紧密的关系, 可以在以后的实验中作进一步研究。最终形成的张力与垂度的数学模型为:

若将参数β1, L代入公式则可以转化为

式中a为皮带重量分布 (kg/m) ;x为测得的垂度值 (mm) ;T 1为张力 (kg) ;β1为皮带与水平地面之间的夹角;L为两托辊之间的距离 (m) 。

由公式 (7) 可以根据检测到的垂度推算出张力。

4 结论

1) 通过实验获取了流量与垂度关系的对应数据;利用最小二乘曲线拟合法研究垂度与流量的实验数据, 得到垂度与流量关系的拟合曲线;通过实验验证分析得出5次多项式拟合曲线为最优垂度与流量拟合曲线。2) 根据最优垂度与流量拟合曲线及张力与流量数学关系建立了皮带张力与垂度的数学模型。

参考文献

[1]李南, 刘秋楠.有关高速传带设备张力控制的应用研究[J].制造业自动化, 2010.2, 32 (2) :132-133.

[2]谢先伟.最小二乘法在煤矿钻机测试中的应用研究[J].制造业自动化, 2011.5, 33 (9) :59-60.

[3]徐跃良.数值分析.成都:西南交通大学出版社, 2005, 08:59-60.

[4]吕喜明, 李明远.最小二乘曲线拟合的MATLAB实现[J].内蒙古民族大学学报 (自然科学版) , 2009.3, 24 (2) :126-127.

二乘拟合篇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.

【二乘拟合】推荐阅读：