最小二乘准则

最小二乘准则（共7篇）

最小二乘准则 篇1

0 引 言

空间直线拟合在工程中具有重要应用[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)为直线上确定的一点,非零向量$\overrightarrow{s}=(m,n,p)$为直线的方向向量,使得目标函数:

${\displaystyle \sum _{i=1}^{n}d}{}^2({\rm P}{}_i,L)(1)$

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

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

${\displaystyle \sum _{i=1}^{n}d}({\rm P}{}_i,L)(2)$

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

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

${\displaystyle \sum _{i=1}^{n}w}{}_{i}d({\rm P}{}_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,求出均值$\overline{X}=(\overline{x},\overline{y},\overline{z})$,拟合直线必过数据中心$(\overline{x},\overline{y},\overline{z})$;

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

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

第4步 由点向式确定直线$\frac{x-\overline{x}}{a}=\frac{y-\overline{y}}{b}=\frac{z-\overline{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个自变量按下面方式排列$\overrightarrow{v}=(x{}_0,y{}_0,z{}_0,a,b,c)$,给出求问题式(3)中直线的算法如下:

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

第2步 利用初始值和MATLAB中fminsearch()对目标函数${\displaystyle \sum _{i=1}^{n}w}{}_{i}d({\rm P}{}_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{array}{l}x=\overline{x}+aty=\overline{y}+btz=\overline{z}+ct\\ a(x-x{}_i)+b(y-y{}_i)+c(z-z{}_i)=0\end{array}(4)$

可以解出$t=\frac{a(x{}_i-\overline{x})+b(y{}_i-\overline{y})+c(z{}_i-\overline{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.

[11]唐焕文,秦学志.最优化方法[M].大连:大连理工大学出版社,2004.

最小二乘准则 篇2

近年来基于导频的信道估计算法得到了广泛研究,R.Ne等[1]研究了慢时变信道环境下最优导频的放置方案,指出慢时变或非时变信道环境下,导频应该等间隔分布。良好的信道估计算法需要同时具有很好的估计精度和良好的跟踪能力。从实现准则来分,信道估计方法可分为最小均方误差(MMSE)[2,3]、最大似然估计(MLE)[4]、最小平方误差(LS)[5]等。MMSE算法[2]性能最好,可以得到精确的信道估计,但是MMSE算法需要知道信道的统计信息,并且具有很高的计算复杂度。线性最小均方误差(LMMSE)[3]信道估计器,是根据最佳低阶秩理论推导出的一种低复杂度的信道估计器,该估计器在信噪比(SNR)较低时性能受限,并且仍然需要已知信道的统计特性和信噪比。最大似然(ML)[4]信道估计算法,无须知道信道的频率相关特性,但是需要已知信道最大时延扩展和信噪比,且对有效信道长度具有敏感性。因此实际系统中常采用LS[5]算法,LS算法计算简单,不需要任何先验信道信息,但是LS算法的性能较差。为了提高LS算法的性能,[6]提出了基于离散傅立叶变换(DFT)的变换域信道估计算法,基于DFT的信道估计算法计算简单,不需要任何信道先验信息,性能优于LS算法。

本文首先分析了OFDM系统中子载波间干扰产生的原因和特性[7],接着分析了基于最小二乘LS的信道估计技术的性能,然后提出了同时考虑信道噪声和载波间干扰的OFDM信道估计算法。根据Tikhonov正则化原理[8],导出了总体最小二乘正则化解法的计算公式。该算法根据子载波间干扰特性,避免了一般总体最小二乘[9]方法带来的设计矩阵结构上的缺陷,并通过迭代计算进一步去除载波间干扰和噪声对系统性能的影响。理论分析和仿真结果表明该算法不仅可以有效的提高系统的误码率性能和信道估计的精度,而且对有限长度的CP仍然能保证通信质量,从而提高了频率利用率。

1 OFDM系统模型

1.1 信道模型

信道为瑞利多径衰落信道,且假设在一个OFDM符号间隔T内不变。在第n个OFDM符号间隔内,离散信道脉冲响应函数可表示为:

其中L为多径延迟扩展用采样时间归一化后的时间单位数,hl为第l条路径信号的衰落因子且为零均值的复高斯随机变量;τl第l条路径的多径时延,fDl为第l条路径的多谱勒频移,它会引起接收信号的载波间干扰(ICI);δ(t)为Kronecker delta函数。对于最大有效延迟数为L的信道,采用抽头系数为复系数的FIR滤波器作为信道的模型,其冲激响应的系数h=[h0,h1,…,hL]T。

1.2 OFDM信号模型

基于导频插入的典型的OFDM系统基带模型如图1所示,基于导频信道估计的OFDM系统首先在映射器中根据系统的调制方式对二进制信息进行符号映射分组并调制,串并变换后以固定周期在所有子载波上插入导频,形成频域序列X[k]。长度为N的序列X[k]被调制到工作子载波上,经过逆离散傅里叶变换(IFFT)模块被转化为时域序列x[n]。

为了消除码间干扰ISI,每帧OFDM符号后要插入保护间隔,保护间隔的长度Ng,Ng应大于等于信道冲击相应的最大延时,保护间隔可以由循环前缀组成,插入循环前缀后,此时OFDM系统的输出信号表示为:

发送信号通过带有加性白高斯噪声(AWGN)的频率选择性信道进行传输。在接收端,串并转换除去CP后,FFT之前的时域接收信号向量yt=[yi,0,yi,1,…,yi,N-1]T可以表示为:

式中G,A和B都是N×N矩阵,它们的组成如下面三个式子所示:

G是循环矩阵;A的最后g列为零。(4)式右边第二项和第三项分别代表ICI和ISI。如果g≥L,A和B为零矩阵,这时不存在ICI和ISI,均衡可以简化为简单的频域单抽头均衡。在CP不足时,残余ISI导致的当前符号内的ICI和前一符号带来的ISI会影响系统性能。这时,需要采取措施来消除这些干扰。

经过FFT变换,接收信号的频域表达式为:

式中Y(k)表示接收端第k个子载波的接收信号值,H(k)表示第k个子载波的信道频率特性值,w(k)表示频域AWGN,I(k)表示第k个子载波上干扰项。

2 最小二乘信道估计性能分析

LS信道估计是最简单的一种信道估计方法。设发送的导频信号X=[X0,X1,…,XP-1],从接收信号的导频位置获得导频符号Y=[Y0,Y1,…,YP-1],依据最小二乘准则:得到导频位置的信道特性:

在OFDM系统中,LS估计的MSE性能取决于导频所占用的频点及功率,当不考虑虚载波时,采用[10]中等间隔等功率的最优导频时,LS估计的MSE为:

式中σ2表示加性高斯白噪声信号功率,Ep为分配给导频的总功率。从式(4),(8)和(10)可以看到,LS信道估计结果仅考虑噪声而没有抑制干扰。如果能够采取一定的方法减小LS信道估计结果中噪声和干扰,那么LS信道估计的MSE可以得到降低。

3 正则总体最小二乘信道估计算法

针对上述ICI和ISI问题,接收端信号经过FFT后可以表示为:

式中W表示导频信号产生的数据矩阵,包括了I-CI和ISI干扰项后,不在为常数对角阵;E表示信道的不确定性对数据矩阵扰动;e表示噪声。

总体最小二乘模型的正则化准则[7]为:

其中α是正则化参数。

构造拉格朗日目标函数:

λ为NP×1拉格朗日乘算子矢量。对式(14)中Φ(E,e,λ)求偏导,并令偏导数等于0得:

求解上面(15)和(16)式得:

将(18)和(19)式代入约束方程得:

由(17)和(19)式得:

将(20)代入(21)化简得:

其中:

下面给出利用(22)、(23)式迭代计算的正则总体最小二乘信道估计算法,具体计算步骤如下:

步骤1用最小二乘算法估计信道的冲激响应,再用信道冲激响应的估计值分别构造出G,A和B;

步骤2初始化数据矩阵W和H(0),按照式(23)计算中间变量β;

步骤3按式(22)、(23)计算第1次迭代结果H(1)和β(1);

步骤4重复步骤3求第i次迭代结果:H(i)和β(i);

步骤5判断||H(i)-H(i-1)||<δ(δ为阈值),若成立则停止迭代;反之,重复步骤4。

上面迭代过程可以得到很精确的信道估计,其性能可由其均方误差来衡量:

4 性能仿真

为了评估文中提出的综合考虑了噪声和干扰的正则总体最小二乘RTLS信道估计算法的有效性,并与LS算法作比较。进行了系统的仿真,仿真参数如下:OFDM系统载频为2.0 GHz,子载波数N=128,16QAM调制,子载波间隔为:7.8125k Hz。各径时延依次为:0,2,4,8和12μs,延迟功率谱服从指数分布。先考察循环前缀充足CP=16时的仿真,即不存在存在I-CI和ISI干扰时改进算法有效性。下面给出仿真图如图2所示:

图2给出了两种算法的MSE随信噪比SNR变化的曲线。因为LS算法的性能受限于噪声,从图2可以看出,在低信噪比时,LS算法抑制噪声的能力比较差所以其MSE性能也差;而RTLS的性能有显著的提高。随着信噪比的增加,RTLS算法与LS算法性能越来越接近,因为随着噪声噪声功率的降低,噪声的影响越来越小。

系统BER随信噪比SNR变化的曲线如图3所示。从图3的BER比较图中也可得出与图2相同的结果,RTLS算法更能抑制噪声。在信噪比低于10d B时,RTLS算法有超过5d B的性能提高,信噪比越低性能提高越明显。

为了验证RTLS算法对干扰的稳健性,在两种循环前缀不足的情况下,系统BER随信噪比SNR变化的比较曲线如图4所示。从图4的比较结果可以看出,LS算法受ICI干扰影响较大,CP长度越小影响越大。因为RTLS总体考虑了噪声和干扰,ICI干扰对其影响不明显。随着SNR的增加,干扰是降低估计精度的主要原因,采用干扰抑制RTLS算法性能比LS算法有约6 d B的性能提高。理论上总体最小二乘算法比LS算法性能有10~15d B的性能提高[7],所以RTLS可大大提高估计的精度。

5 结束语

信道估计是实现OFDM系统性能的关键。由于总体最小二乘算法可以同时兼顾噪声和干扰,本文通过对频率选择性衰落信道下的OFDM系统ISI和ICI产生机理分析,避免了一般总体最小二乘算法设计矩阵结构上的缺陷。根据Tikhonov正则化原理,提出了总体最小二乘正则化信道估计算法。仿真结果都验证了该算法的有效性,尤其在存在干扰时改进明显,比普通最小二乘算法估计性能提高5d B以上。因此对噪声和干扰综合考虑可以提升信道估计算法鲁棒性。

参考文献

[1]NEGI R,CIOFFI J.Pilot tone selection for ehannel estimation in amobile OFDM system[J].IEEE Trans oll Consumer Electronics.1998,44(3):1122-1128.

[2]Coleri Sinem,Ergen Mustafa,Puf Anuj,Bahai Ahmad.Channelestimation techniques based on pilot arrangement in OFDM systems[J].IEEE Trans.on Broadcasting,2002,48(3):223-229.

[3]Y.Li,L.J.Cimini,and N.R.Sollenberger,Robust channel estimation forOFDM systems with rapid dispersive fading channels[J].IEEE Trans.Commun.,1998,46(7):902-915.

[4]Deneire L,Vandenameele P,Van der Perre L,et al.AlowcomplexityMLchannel estimator for OFDM[J].IEEE Transactions on Communications,2003,51(2):135-140.

[5]J.-C.Lin,Least-squarechannelestimationformobileOFDMcommunicationon time-varying frequency-selective fading channels[J].IEEE Trans.Veh.Technol.,2008,57(6):3538-3550.

[6]ZhaoY,HuangA.A novel channel estimation method for OFDMmobilecommunication systems based on pilot signals and transform-domainprocessing[C].ProceedingsofIEEE47thVehicularTechnologyConference.IEEE,1997:2089-2093.

[7]ZhongWei,MaoZhigang.EfficientTimeDomainResidualISICancellationfor OFDMBased WLANSystems[J].IEEE,2006,52(2):321-326.

[8]AmirBeck,AharonBen-Tal.OntheSolutionoftheTikhonovRegularizationoftheTotalLeastSquaresProblem[J].SIAMJ.OPTIM,2006,17(1):98-118.

[9]S.Chandrasekaran,MGu,A.h.Sayed,and K.E.Schubert.The degeneratebounded errors-in-variables model[J].SIAM J.Matrix.Anal.Appl.,2001,23(1):138-166.

轴承沟道形状误差的最小二乘评定 篇3

轴承主要用于确定旋转件与固定件相对运动,是起支承或导向作用的典型零部件,在机器及各类机械中占有重要地位。因此,研究轴承沟道形状误差评定对提高轴承的精度及整个机器或机械的性能有重要的意义。目前多使用圆度、线轮廓度来评价轴承沟道形状误差,即用一条线轮廓和底圆形成的几何要素的误差对轴承沟道的形状误差进行评定,而对轴承沟道整体形状误差的评定国内外很少报道。本文依据最小二乘形状误差的定义,提出了轴承沟道形状误差的最小二乘评定。

1 最小二乘算法过程及步骤

1.1 提取测量点坐标

将工件水平放置在工作台上,利用三坐标测量机对轴承沟道的每条线轮廓进行采样取点,从中取得的测量点为Pij(xij,yij,zij)(i表示测量线轮廓数,i=1,2,…,M;j表示每条线轮廓上的测量点数,j=1,2,…,N)。

1.2 求空间每条线轮廓的中心点

通过三坐标测量机对轴承沟道每条线轮廓进行测量并得到采样点,用最小二乘法对每条线轮廓求其半径及圆心的位置,但轴承沟道轮廓是非整圆,不能直接用整圆的最小二乘法公式,需要坐标变换求空间每条线轮廓的中心点。

1.2.1 坐标变换

在空间坐标系中,直接用最小二乘法求每条线轮廓的中心点比较复杂,为了便于计算,将空间内的线轮廓绕z轴转θi角度至yoz平面。其中:

undefined。

空间线轮廓上各点坐标通过坐标(xij,yij,zij)变换至yoz平面上的新坐标为:undefined。

1.2.2 求所有平面线轮廓的圆心

设平面线轮廓的最小二乘拟合方程为:

(y-a′i)2+(z-b′i)2=rundefined。 (1)

其圆心为(a′i,b′i),半径为ri,平面线轮廓最小二乘拟合方程见图1。对式(1)展开并整理得:

y2+z2-2ai′y-2bi′z+(ai′2+bi′2-rundefined)=0 。 (2)

令c=ai′2+bi′2-rundefined,由式(2)得:

y2+z2-2ai′y-2bi′z+c=0 。 (3)

以每条线轮廓上测点坐标(y′ij,z′ij)代入式(3),方程式明显不等于零,而为:

Δij=y′2ij+z′ij2-2a′iy′ij-2b′iz′ij+c 。 (4)

其中:Δij为实际线轮廓上各点与理论线轮廓上对应点的函数偏差值。设:

当偏差值为最小值时,理论线轮廓逼近实际线轮廓效果最佳。令:

依据求最小值的方法,对系数a′i、b′i、c求偏导数后得出:

其中:undefined。

再由式(3)得出:undefined,在求出待定系数a′i,b′i,c和ri后,可以得出平面线轮廓上所有中心点O′i(a′i,b′i)。

1.2.3 坐标变换

将上一步得到的平面线轮廓的中心点代入下式可以得到空间线轮廓上的中心点Oi(ai,bi,ci):

1.3 最小二乘法拟合空间圆

将每条线轮廓的中心点用最小二乘法拟合一个空间圆,但空间圆没有特定的方程, 只能由空间圆和空间平面的方程联立来表示。距差与z轴的关系如图2所示,由于拟合的空间圆不在一个平面内,其圆心(中心点)为z向的一系列点,r′1,r′2,r′0分别为被测轮廓的测点到中心点的距离,通过式r′21 -r′20 =Δz2,证明空间坐标系中被测轮廓的测点与中心点的距离r′相对于z轴的距差(Δz)大得多,故Δz对r′的影响可以忽略不计。因此可得出中心点在空间平面上,而空间圆也在空间平面上且平行于xoy平面。

设空间圆的方程为:

(x-a0)2+(y-b0)2+(z-c0)2=R2。

其中:圆心为(a0,b0,c0),半径为R。展开并整理得:

x2+y2+z2-2a0x-2b0y-2c0z+aundefined+bundefined+cundefined=R2 。 (10)

设A=-2a0,B=-2b0,C=-2c0,D=aundefined+bundefined+c20-R2。由式(6)得:

x2+y2+z2+Ax+By+Cz+D=0 。 (11)

将最小二乘法拟合每条线轮廓得到的中心点坐标Oi(ai,bi,ci)代入式(11),方程式明显不等于零,而为:

δi=aundefined+bundefined+cundefined+Aai+Bbi+Cci+D 。 (12)

其中:δi为空间圆上各点的函数偏差值,则:

当undefined为最小值时,空间圆拟合的效果最佳。令:

undefined。 (14)

依据求最小值的方法,对系数A、B、C和D求偏导数后得出关于方程的系数A、B、C和D的线性方程如下:

解得A、B、C、D后,代入计算得:

undefined。

即可得到圆心坐标(a0,b0,c0)和圆半径R。

1.4 求解空间圆与每条线轮廓平面的交点坐标

空间圆与线轮廓平面的方程为:

如式(17)所示,每条线轮廓所在平面与z轴平行且共面,那么该平面的方向数C等于零;最小二乘拟合的空间圆与xoy平面平行且共面,所以每条线轮廓所在平面与空间圆相交,可以得到所有的交点坐标值Li(xi,yi,zi)。而交点坐标值Li(xi,yi,zi)既在平面内又在空间圆内,并且交点与所对应线轮廓的测点在一个象限内。

1.5 轴承沟道形状误差

轴承沟道形状误差模型见图3。

由图3可知,每条线轮廓上的测量点Pij(xij,yij,zij)至所对应的每个交点Li(xi,yi,zi)的距离为:

由轴承沟道形状误差的定义可以知道,每条线轮廓上各测量点至所对应交点的距离中,存在最大值、最小值和极差(最大值与最小值的差值)。因而可以得到M个最大值max(dij)、最小值min(dij)和极差Δdijcha。

以最大值、最小值所在的线轮廓绕z轴转360°形成圆环区域,比较所有的极差值,其中最大值为max(Δdijcha),即得到包容整个轴承沟道的最小二乘形状误差值TLSM为:

TLSM=max{Δdijcha} 。 (19)

2 结论

采用最小二乘算法进行轴承沟道形状误差评定,计算简单,完全满足最小二乘法的评定标准。另外,评定方法不仅可以给出轴承沟道形状误差,还能适用于现代机械的精密设计和精密测量。这种评定方法简单正确,收敛速度快,易于计算机程序实现,具有较好的推广应用价值。

参考文献

[1]王伦.轴承基本知识[M].北京:机械工业出版社,2002.

[2]《数学手册》编写组.数学手册[M].北京:高等教育出版社,2004.

[3]沈世德,徐辛伯.最小二乘圆弧法在图像分析中的应用[J].机械设计与制造,1999,10(5):46-47.

最小二乘准则 篇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

卫星导航定位在现代军事和经济领域扮演着不可或缺的重要角色, 卫星导航系统除了提供导航功能外, 还必须具有在系统不能使用时及时向用户发出告警的能力, 这种能力叫做系统的完好性。完好性和定位精度是卫星导航中两个重要指标。当卫星导航用于航空、武器系统等生命相关场合时, 完好性成为决定性指标。因此, 卫星导航系统的完好性问题已成为卫星导航领域重点解决的问题, 并将存在于卫星导航整个生命周期。

引起接收机伪距观测量偏差的主要因素有较大的卫星钟漂、导航电文的不正确上传、各种欺骗与干扰以及卫星组成部分的故障等, 这里统称为卫星故障, 在GPS系统中这种故障发生的概率大约为10～4/小时。对这种故障, 导航系统本身做出告警反应的时间较长, 如GPS系统的反应时间为15分钟到数小时, 对于生命相关场合显然太慢, 不能满足需求。因此, 卫星故障的快速检测只有在用户端进行, 因而出现了接收机自主完好性检测技术 (Receiver Autonomous Integrity Monitoring, 英文简写RAIM) 。本文针对导航接收机常用的位置解算方法, 即最小二乘法, 设计了基于最小二乘残差法的RAIM检测方案, 并对方案进行了仿真分析。

1 基于最小二乘残差法的RAIM检测方案

导航接收机位置解算常用最小二乘法, 针对最小二乘位置解算方法, 设计了相应的RAIM检测方案, 即基于最小二乘残差法的RAIM检测方案。基于最小二乘残差法的RAIM方案除了对当前可视卫星的几何结构有很强的依赖性, 还对可视卫星数目提出了很高要求。如果接收机要完成故障检测, 最少需要5颗可视卫星;如果接收机要完成故障排除, 则最少需要6颗卫星。

1.1 最小二乘位置解算法[5,6]

对于伪距定位方程:

undefined. (1)

在X0= (x0, y0, z0, B0) ′ (B=cΔt, 看作一个未知数, 以米为单位) 处进行泰勒级数展开[4,5], 略去二次项及高次项得到:

undefined

令undefined

则:undefined. (3)

undefined. (4)

undefined. (5)

将式 (3) (4) (5) 代入式 (2) , 方程可整理为:

ρi-ρ0=αi (x-x0) +βi (y-y0) +γi (z-z0) + (B-B0)

=αiΔx+βiΔy+γiΔz+ΔB. (6)

测量四颗卫星的伪距, 可获得四个类似式 (6) 的方程, 联立这四个方程得到一个四元线性方程组:

令

undefined

ΔX=ΔX ΔY ΔZ ΔBT.

Y=T.

则式 (7) 变为: GΔX=Y. (8)

解此方程得: ΔX=G-1Y. (9)

若测量伪距的卫星数n大于4, 则G为n×4的矩阵, 解此方程得:

ΔX= (GTG) -1GTY. (10)

所以 X=X0+ΔX. (11)

再以新的X作为X0, 继续上述解算, 迭代数次后, 当满足设定的误差时, 迭代结束, 得到所求的解, 实践经验表明, 迭代次数为5左右。

1.2 故障检测与识别

基于最小二乘位置解算原理, 利用最小二乘迭代解算过程的中间结果, 就可以完成RAIM故障检测与识别。引入RAIM故障检测与识别的最小二乘位置解算过程如图1所示, 其中完好性处理部分的详细流程如图2所示。

完好性处理部分首先对观测量Y进行故障检测后, 观测量Y才能进入最小二乘迭代运算过程。完好性处理主要完成可视卫星数目的判断、水平定位误差保护级HPL、检验统计量及其限值的计算, 并对故障进行检测、识别等。具体过程如下:在RAIM检测之前, 有两个重要的限值需要计算, 它们是水平定位误差保护限值HAL和故障检测门限T。水平定位误差保护限值和故障检测门限可通过文献中给出的算法计算得到, 或者可通过文献提出的基于最小风险准则来确定。

在RAIM检测开始时, 首先根据获取的观测信息, 判断可视卫星数目是否达到故障检测与识别要求, 判断可视卫星数目的流程如图3所示。判断可视卫星的主要依据是计算过境卫星的仰角, 如果仰角大于5度, 认为可视。理论上卫星导航系统中最多可视卫星为12颗, 所以判断可视卫星时最多进行12次循环判断。

如果前面判断的可视卫星数目符合要求, 对水平定位误差保护级HPL进行计算, 并与HAL比较, 看是否超限, 以保证RAIM检测的有效性。当HPL小于HAL时, 对检验统计量进行计算。如果检验统计量小于检测门限值T, 则认为无故障, 否则, 进入故障识别与隔离阶段。如果可视卫星数目只有5颗, 那么只能完成故障的检测, 并对检测到的故障向用户发出告警;如果可视卫星数目大于等于6颗, 则可以对检测到的故障进行隔离。最后, 输出完好性处理结果。

2 方案仿真

为了验证基于最小二乘残差法的RAIM方案对故障检测与识别的性能, 利用STK5.0软件产生星座数据, 采用MATLAB软件对该方案进行仿真。由于GPS系统已有较为准确的误差统计数据, 因此仿真采用的卫星星座是GPS卫星星座。利用了在GPS系统时间1 Jan 2010 13:15:00.00至1 Jan 2010 13:17:00.00之间, 在某地区所有的可视卫星, 它们分别是GPS星座系统第1、2、3、6、9、11、15、20颗卫星。仿真条件如下:故障检测误警率为0.005, 由此计算得到故障检测门限值为10.343 m。设观测噪声均值为0, 标准差为5 m。在观测时段内第30秒时刻, 对其中一颗卫星对应的伪距加入一个100 m突变, 模拟此时卫星出现故障, 相应的伪距出现异常。

在上面仿真条件下, 仿真结果如图4所示。横轴为观测时间, 单位为秒。左边纵轴表示检验统计量, 单位为米, 右边纵轴为观测伪距残差, 单位为米。从仿真结果图上可以看出, 在30秒时刻观测伪距发生异常时, 相应的伪距残差出现突变, 从而导致了检验统计量的突变。图中所示的检验统计量在30秒时刻已经超过了故障检测门限值, 因此, 此时的故障完全能够被检测出来。仿真结果表明, 本文所设计的基于最小二乘残差法的RAIM故障检测方案是可行的, 有效的。

3 结束语

完好的系统故障检测性能、快速的故障反馈能力是现代卫星导航定位系统接收机自主完好性检测的发展趋势。本文提出的基于接收机自主完好性检测方案具有告警时问短、运算量小的特点。接收机位置解算如果采用最小二乘法, 就可以很容易地将该方案嵌入其中。我国二代卫星导航定位系统正在建设当中, 建议在系统建设中考虑好系统的完好性检测扩展能力, 在接收机的设计中加入自主完好性检测功能。

摘要：为了使导航接收机的自主完好性检测算法与最小二乘位置解算法结合起来, 提供准确可靠的定位定时服务, 对最小二乘位置解算法和基于最小二乘残差法的接收机自主完好性检测算法进行了研究。针对导航接收机最小二乘位置解算法, 提出了基于最小二乘残差法的卫星导航接收机自主完好性检测方案, 并对方案进行了仿真分析, 仿真结果表明方案是可行的, 有效的。

关键词：导航接收机,自主完好性,最小二乘法

参考文献

[1]Brown R G.Global Positioning System Theory and Appli-cations[M].American:American Institute of Aeronauticsand Astronantics, 1996:143-164.

[2]Xie G, Pullen S, Luo M, et al.Integrity Design and Up-dated Test Results for the Stanford LAAS Integrity MonitorTestbed[C].Albuquerque, NM:Proceedings of ION 57thAnnual Meeting and CIGTF 20th Biennial Guidance TestSymposium, 2001:681-693.

[3]冯庆堂, 沈林成, 常文森.一种监测GPS完好性的新方法[J].导航与控制, 2004, 3 (1) :4-5.

[4]Bradford W, Parkinson, Penina Axelrad.AutonomousGPS Integrity Monitoring Using the Pseudorange Residual[J].Journal of the Institute of Navigation, 1988, 35 (2) :456-460.

[5]王惠南.GPS导航原理与应用[M].北京:科学出版社, 2003.

[6]刘基余.GPS卫星导航定位原理与方法[M].北京:科学出版社, 2003.

[7]黄晓瑞, 田巍, 李波.GPS接收机的自主完善性监测算法研究[J].遥测遥控, 2003, 24 (1) :1-3.

最小二乘准则 篇6

考虑下列二阶抛物问题

$\{\begin{array}{l}u{}_t(x,t)-\nabla \cdot (a(x)\nabla u(x,t))=f(x,t),\\ (x,t)\in \Omega \times J,\\ u(x,t)=0,(x,t)\in \Gamma \times J,\\ u(x,0)=u{}_0(x),x\in \Omega ‚\end{array}(1)$

(1)式中J=[0,T],Ω⊂Rn为有界凸域,边界Γ是Lipschitz连续的。假定a(x)有界,即存在常数α*和α*,使得

0≤α*≤a(x)≤α*,∀x∈Ω (2)

对上述问题人们提出了许多数值模拟格式,如有限元格式[1],混合有限元格式[2],扩展混合有限元格式[3],最小二乘混合有限元格式[4]等。现将最小二乘思想和扩展混合有限元方法相结合,提出并分析了最小二乘扩展混合有限元格式。该格式继承了最小二乘和扩展混合元方法的优点,即格式能同时高精度逼近未知函数,未知函数的梯度以及流体的通量,所形成的有限元方程组是对称正定的,有限元空间无须满足LBB相容性条件。数值算例说明了该格式的有效性。

1 最小二乘扩展混合元格式

定义空间Λ=(L2(Ω))n,V=H${}_0^1$(Ω),W=H(div;Ω)={q∈(L2(Ω))n,∇·q∈L2(Ω)},相应范数为$\parallel \mu \parallel {}_{0,\Omega }=(\mu ,\mu ){}^{\frac{1}{2}},\parallel v\parallel {}_{1,\Omega }=(\parallel v\parallel {}_{0,\Omega }^2+\parallel \nabla v\parallel {}_{0,\Omega }^2){}^{\frac{1}{2}},\parallel q\parallel {}_{{\rm H}(\text{d}\text{i}\text{v};\Omega )}=(\parallel q\parallel {}_{0,\Omega }^2+\parallel \nabla \cdot q\parallel {}_{0,\Omega }^2){}^{\frac{1}{2}}$,其中(·,·)表示空间(L2(Ω))n中的内积。作时间剖分τ=T/N,N为正的常数,在时间tn=nτ,n=1,2,…,N处逼近方程的解。$\text{Ω}={\displaystyle \underset{{\rm K}\in \tau {}_h}{\cup }{\rm K}}$是Ω的一个剖分,h=max{diam(K):K∈Th}。设Vh,Λh和Wh分别是V,Λ和W的有限维子空间,具有通常逼近性质。Vh,Λh和Wh的一个标准的选法是k,r和s次的分片多项式空间。

记λ=-∇u,σ=aλ,并用一阶差分近似ut,则式(1)可变成离散形式I

$\{\begin{array}{l}\overline{\partial }{}_{t}u{}^n=\nabla \cdot \sigma {}^n=f{}^n+R{}_1^n,x\in \Omega ,\\ \lambda {}^n+\nabla u{}^n=0,x\in \Omega ,\\ a\lambda {}^n-\sigma {}^n=0,x\in \Omega ,\\ u{}^n(x)=0,x\in \Omega ,\\ u{}^0(x)=u{}_0(x)x\in \Omega \end{array}(3)$

(3)式中$\overline{\partial }{}_{t}u{}^n=(u{}^n-u{}^{n-1})/\tau ,R{}_1^n=\overline{\partial }{}_{t}u{}^n-u{}_t^n$。

对于(u,λ,σ)∈V×Λ×W和(v,u,q)∈V×Λ×W,定义双线性形式

$\begin{array}{l}a{}_{{\rm N}}(u,\lambda ,\sigma ;\upsilon ,\mu ,q)=(u+\tau \nabla \cdot \sigma ,\upsilon +\tau \nabla \cdot q)+\\ \tau (\lambda +\nabla u,\mu +\nabla \upsilon )+\tau (a{}^{-1}(a\lambda -\sigma ),a\mu -q)(4)\end{array}$

根据离散形式式(3),省略时间截断误差定义最小二乘扩展混合元格式

格式Ⅰ 给定初始值u${}_h^0$∈Vh,对于n=1,2,…,N,求(u,λ,σ)∈Vh×Λh×Wh,使得

$\begin{array}{l}a{}_{{\rm N}}(u{}_n^n,\lambda {}_n^n,\sigma {}_n^n;\upsilon {}_h,\mu {}_h,q{}_h)=(u{}_h^{n-1}+\tau f{}^n,\upsilon {}_h+\\ \tau \nabla \cdot q{}_h),\forall (\upsilon {}_h,\mu {}_h,q{}_h)\in V{}_h\times \Lambda {}_h\times W{}_h。(5)\end{array}$

2 解的存在唯一性

为证明格式Ⅰ解的唯一性,首先说明双线性形式aN(·;·)是强制的。

引理1 对任意的(v,μq)∈V×A×W存在c>0,使得aN(v,μ,q;v,μq)≥c(‖v‖${}_{1,n}^2$+‖μ‖${}_{0,n}^2$+‖q‖H(div n))。 (6)

证明 引入适当选取的常数β>0,式(4)可整理为

$\begin{array}{l}a{}_{{\rm N}}(\upsilon ,\mu ,q;\upsilon ,\mu ,q)=(\tau \nabla \cdot q+(1-\beta )\upsilon ,\\ \tau \nabla \cdot q+(1-\beta )\upsilon )+\tau (q-a(\mu +\beta \upsilon ),\\ a{}^{-1}(q-a(\mu +\beta \upsilon )))+\\ 2\tau (\mu +(1-a\beta )\nabla \upsilon ,\mu +(1-a\beta )\nabla \upsilon )+\\ (2\beta -\beta {}^2)(\upsilon ,\upsilon )。(7)\end{array}$

由于υ(x)=0,x∈Γ,所以由Poincaré-Friedrichs不等式知,存在正常数CF,使得

‖v‖${}_{0,n}^2$≤CF‖∇v‖${}_{0,n}^2$ (8)

结合式(2)有

$\begin{array}{l}a{}_{{\rm N}}(\upsilon ,\mu ,q;\upsilon ,\mu ,q)\ge \beta (2(\tau \alpha {}_{*}+C{}_F)-\\ \beta (\tau \alpha {}^{*2}+\tau \alpha {}^{*}+C{}_F))\parallel \nabla \upsilon \parallel {}_{0,\Omega }^2。(9)\end{array}$

取β=τα*/(τα*2+τα*+CF),则

$a{}_{{\rm N}}(\upsilon ,\mu ,q;\upsilon ,\mu ,q)\ge C\parallel \nabla \upsilon \parallel {}_{0,\Omega }^2\ge C\parallel \upsilon \parallel {}_{1,\Omega }^2。(10)$

从而由式(10)和式(4),知

$\begin{array}{l}\parallel \mu \parallel {}_{0,\Omega }^2\le 2\parallel \mu +\nabla \upsilon \parallel {}_{0,\Omega }^2+2\parallel \nabla \upsilon \parallel {}_{0,\Omega }^2\le Ca{}_{{\rm N}}\\ (\upsilon ,\mu ,q;\upsilon ,\mu ,q),(11)\\ \parallel q\parallel {}_{0,\Omega }^2\le 2\parallel a\mu -q\parallel {}_{0,\Omega }^2+2\parallel a\mu \parallel {}_{0,\Omega }^2\le Ca{}_{{\rm N}}\\ (\upsilon ,\mu ,q;\upsilon ,\mu ,q),(12)\\ \parallel \nabla \cdot q\parallel {}_{0,\Omega }^2\le \frac{2}{\tau }\parallel \upsilon +\tau \nabla \cdot q\parallel {}_{0,\Omega }^2+\frac{2}{\tau }\parallel \upsilon \parallel {}_{0,\Omega }^2\le Ca{}_{{\rm N}}\\ (\upsilon ,\mu ,q;\upsilon ,\mu ,q),(13)\end{array}$

结合式(10)—式(13)便证得式(6)成立。

定理1 设f∈L2(Ω),则格式Ⅰ存在唯一解(μ${}_h^n$,λ${}_h^n$,σ${}_h^n$)∈Vh×Λh×Wh。

证明 考虑空间V×Λ×W,其范数定义为‖υ‖1,Ω+‖μ‖0,Ω+‖Q‖H(div;Ω)。由引理1知双线民生形式αN(·;·)在空间Vh×Λh×Wh中是强制的;而且由Cauchy不等式易证αN(·;·)在Vh×Λh×Wh中是有界的。因此由Lax-Milgram引理知格式Ⅰ解存在唯一。

3 数值算列

为验证所提格式的有效性,在区域[0,1]×[0,1]上考虑方程

$\{\begin{array}{l}u{}_t-(au{}_x){}_x=f,(x,t)\in (0,1)\times (0,1],\\ u(0,t)=u(1,t)=0,t\in [0,1],\\ u(x,0)=\sin \pi x,x\in [0,1]\end{array}(14)$

(14)式中$a=\frac{1}{100},f=\left(1+\frac{\pi {}^2}{100}\right)\text{e}{}^t\sin \pi x$。

容易验证方程的真解为$u=\text{e}{}^t\sin \pi x,\lambda =-\pi \text{e}{}^t\cos \pi x,\sigma =-\frac{1}{100}\pi \text{e}{}^t\cos \pi x$。取Λh为分片常数多项式空间,Vh和Wh为分片线性多项式空间,用MATLAB编程计算近似解uh,λh和σh,绘出真解和最小二乘扩展混合元解的比较图,见图1—图3。

摘要：对二阶抛物问题提出了最小二乘扩展混合元数值模拟格式。数值算例验证了格式的有效性。

关键词：二阶抛物问题,最小二乘,扩展混合有限元方法

参考文献

[1]Thome V.Galerkin finite element methods for parabolic problems.Lecture Notes in Mathematics,Springer,1984

[2]Johnson C V.Thome,Error estimates for some mixed finite element methods for parabolic type problem.RAIRO Numer Anal,1981;15:41—78

[3] Guo Ling,Chen Huanzhen.An expanded characteristic mixed finite element method for a convection-dominated transport problem.Journal of Coputational Mthematics,2005;23:479—490

最小二乘准则 篇7

近年来,加权总体最小二乘是测绘学科数据处理领域的热点研究问题。针对观测向量L和系数矩阵A同时存在误差的情况,在1980 年,Golub和Van Loan[1]从数值分析的观点首次对这种情况进行了全面的分析,并正式命名为总体最小二乘。通过对总体最小二乘问题的进一步研究,发现观测向量L和系数矩阵A并不可能都是独立等精度的,而在总体最小二乘算法的求解过程中是将两者当做独立等精度的情况进行求解的, 面对这一问题,Schaffrin和Wieser[2]提出了加权总体最小二乘方法。目前,比较新的加权总体最小二乘算法除了Schaffrin和Wieser[2]提出的Schaffrin-Wieser算法,还有就是2012 年Mahboub[3]提出的Mahboub算法,以及2012 年Amiri-Simkooei[4]提出的Amiri-Simkooei算法。本文主要研究了这三种算法,通过近似垂直、近似水平和普通直线拟合的模拟算例以及在六参数坐标转换模型中的应用对这三种算法进行了比较和分析。对今后在实际工程中的应用,具有一定的理论价值和现实意义。

1 加权总体最小二乘算法

1.1 Schaffrin-Wieser算法

误差方程式为:

式中,y ∈ Rn ×1为观测向量,A ∈ Rn × m为系数矩阵,∈ Rm ×1为未知参数,ey∈ Rn ×1、EA∈ Rn × m分别为观测向量和系数矩阵的误差向量,Qy∈ Rn × n、QA∈ Rmn × mn分别为观测向量和系数矩阵的协方差矩阵,其中QA被分解为Q0∈ Rm × m、Qx∈ Rn × n,Py∈ Rn × n、PA∈Rmn × mn分别为观测向量和系数矩阵的权阵,vec(·) 表示矩阵的拉直运算,表示kronecker积。

加权总体最小二乘的平差准则或目标函数为:

误差模型在该准则下的拉格朗日极值函数为:

通过推导,得Schaffrin-Wieser算法如下[2]:

(6) 当时,计算结束。

1. 2 Mahboub算法

2012 年Mahboub[3]提出了五条原则来构造系数矩阵的协方差矩阵[5]:

(1) 如果系数阵的某个元素重复出现,认为这两个元素100% 相关,因此这两个元素之间的协方差等于重复元素的自方差;

(2) 假如系数矩阵的某个元素以其相反数的形式重复,认定这两个元素100% 负相关,因此这两个元素之间的协方差等于重复元素的自方差的相反数;

(3) 如果系数阵的某个元素是常数,认为其方差为零;

(4) 系数阵中两个不同元素,若两者明显相关,他们的协方差即为其实际值,否则为零;

(5) 上述规则在同方差情况中同样适用,若元素是随机数只需用数字1 作为其方差,若是常数其自方差为零。

利用上述5 条规则,Mahboub在重新构造系数矩阵的协方差阵的基础上,提出了一种新的算法。具体解算步骤如下[3]:

(1) 根据上述五条规则生成QA;

(5) 当时,计算结束。

1. 3 Amiri-Simkooei算法

2012 年,Amiri-Simkooei根据经典最小二乘理论,并结合文献[3] 的迭代方法,提出了一种类似于最小二乘解形式的加权总体最小二乘算法。具体的解算步骤如下[4]:

(3) 当时,计算结束。

2 算例及分析

2. 1 直线拟合

为了比较各种加权总体最小二乘方法计算得到的待定系数估值的差别,本文利用MATLAB软件编程进行数值模拟实验,分别用上述三种算法,以及总体最小二乘方法( TLS)、最小二乘方法( LS)和混合总体最小二乘方法( LS-TLS) 进行求解系数,对直线的不同情况进行拟合,最后比较分析各种方法在不同情况下的差异。 本文将迭代的Frobenius范数差的阈值 δ0设为10- 8。

(1) 近似垂直的直线

模拟一条近似垂直于x轴的直线y=ax+b,其直线的设计模型为y=2000x-2000,其自变量x从0.1到10之间每隔0.1模拟一个观测点,因变量y的范围为(-1800,18000),共100组观测数据,然后在每组观测点的x和y的坐标上均加上均值0,方差为σ2i=1/i×diag(0.01,0.01)的随机误差,其中i∈(1,…,100)表示点号。作1000次随机试验,利用MATLAB软件计算,各种方法的参数估计结果见表1。图1为各方法求出的参数估计值与真值差值的范数分布图,图2为利用Schaffrin-Wieser算法、Mahboub算法、Amiri-Simkooei算法和TLS算法对近似垂直直线拟合做100次随机试验,得到的参数估计值与真值的差值分布图。

从表1 的数据和图1 中可以看出,在近似垂直直线拟合中使用加权总体最小二乘方法( WTLS)进行求解,拟合的结果比TLS、LS以及LS-TLS方法拟合的结果都要好。而LS-TLS方法拟合的结果比TLS拟合的结果好。这是由于在LS-TLS方法中考虑到了在系数矩阵中存在常数项,不需要进行误差改正,从而使得得到的EA比较可靠。而在加权总体最小二乘方法中,不仅考虑到系数矩阵中某些元素是常数,不存在误差,而且考虑到各元素之间是不等精度的,以及各元素之间存在相互联系,所以采用WTLS解算的结果更加接近真值。从图2 中可以看出,在近似垂直直线拟合中Schaffrin-Wieser算法、Mahboub算法、Amiri-Simkooei算法这3 种算法所求得的参数估值波动较小,而TLS算法求得的参数估值与真值的差值波动较大,这也直观反映出,在系数矩阵和观测向量中各元素精度不同时,采用加权总体最小二乘算法更好。

(2) 近似水平的直线

模拟一条近似水平的直线y = ax + b,其直线设计模型为y = x /500 + 1,自变量x从0. 1 到10 之间每隔0. 1 模拟一个观测点,因变量y的范围为(1. 0002,1. 02),共100 组观测数据,然后在每组观测点的x和y的坐标上均加上均值为0,方差为σi2= 1 / i × diag(0. 01,0. 01) 的随机误差,其中i ∈(1,…,100) 表示点号,作1000 次随机试验。表2为各种方法对近似水平直线模拟的估计值。图3 为各方法求出的参数估计值与真值差值的范数分布图,图4 为利用Schaffrin-Wieser算法、Mahboub算法、Amiri-Simkooei算法和TLS算法对近似水平直线拟合做100 次随机试验,得到的参数估计值与真值的差值分布图。

从表2的数据和图3中可以看出,在近似水平直线拟合中使用加权总体最小二乘方法(WTLS)进行求解,拟合的结果比TLS、LS以及LS-TLS方法拟合的结果都要好。从图4中也能较直观地看出Schaffrin-Wieser算法、Mahboub算法、Amiri-Simkooei算法这3种算法所求得的参数估值波动较小,比较接近真值,且这3种算法所求得的参数估值的结果都较为接近。

(3) 普通直线

模拟普通直线y = ax + b,直线的设计斜率范围为( - 5 ≤ a ≤5,且a ≠0),b = 1 自变量x从0. 1 到10 之间每隔0. 1 模拟一个观测点,共100 组观测数据,然后在每组观测点的x和y的坐标上均加上均值为0,方差为 σi2= 1 / i × diag(0. 01,0. 01) 的随机误差,其中i ∈ (1,…,100) 表示点号。表3 代表不同斜率直线的截距的模拟结果,b的真值恒为1。表4 代表不同斜率直线的斜率参数模拟结果。

从表3 和表4 的数据可以看出,在普通直线拟合中使用加权总体最小二乘方法( WTLS) 进行求解,拟合的结果更为接近实际真值。而三种加权算法得到的参数估值也基本相同。从以上直线拟合的算例中可以看出,上述三种加权算法虽然表现形式不同,但实际上是等价的,它们都是通过构造Lagrange目标函数,在这一基础下,通过一定目的的公式变形得到参数估值的迭代求解公式。

2. 2 坐标转换算例及分析

六参数坐标转换模型可写成:

当控制点个数为n(n ≥ 3) 时,六参数转换模型可以表示为:

其中(Xs,Ys)、(Xt,Yt) 分别为原始坐标系和目标坐标系中控制点的坐标值;sX、sY为沿X轴和Y轴的缩放比例,β 为旋转参数,c1、c2 为平移参数,δ 为关于坐标系两个坐标轴旋转的尺度因子。

通过MATLAB软件进行数值模拟实验。选取13 个点的二维坐标,设转换参数a1 = 0. 9,b1 = -0. 8,c1 = 1,a2 = 0. 6,b2 = 0. 7,c2 = 5,按六参数转换模型得到新坐标系下的坐标值(x,y),实验数据如表5 所示。在两套坐标上均加上均值为0,标准差为 σi= i × diag(0. 001,0. 001) 的随机误差,其中i表示点号,以使各坐标点具有不同的精度。

分别用Mahboub算法、Amiri-Simkooei算法、TLS、LS、LS-TLS方法求解坐标转换参数, 由于Schaffrin-Wieser算法在系数矩阵协方差结构上的限制,无法应用于二维坐标转换模型,因此在本算例中未使用该方法。最后将各方法求解的参数估值列于表6。

从表6 的数据和图5 可以看出使用加权总体最小二乘方法(WTLS) 进行求解,拟合的结果更为接近实际真值。而WTLS方法得到的结果比LS-TLS得到的结果更好,则可以看出,由于在WTLS方法中引入了PA、Py权阵,达到部分修改系数矩阵A的作用,即只修改数据元素,而固定常数元素,以及使用了Mahboub算法中构造系数矩阵协因数阵的五条原则[3],使得到的系数矩阵的残差阵EA更为准确,从而较LS-TLS方法得到更合理的模型。

3 总结

上述三种算法有相同点,同时也存在着明显的不同点以及适用范围:

(1) Schaffrin-Wieser算法的特点在于系数矩阵的定权方面,提出了QA构造法QA= Q0Qx,其中QA被表示成两个矩阵的直积形式[2],能够使迭代计算更为简便。但是该方法的缺点在于拆分后的矩阵Q0可能是奇异的,从而导致QA也是奇异的,以及当QA中行向量和列向量没有固定形式时,无法使用,例如在二维坐标转换模型应用中;

(2) Mahboub算法是在Schaffrin-Wieser算法的基础上提出的,在算法中利用五条原则来重新构造系数矩阵的协方差矩阵,这使得在变量误差模型中系数矩阵的协方差阵结构得到更好的保持。这一算法也使得实施过程中协方差阵不必限制于其他条件,如Schaffrin-Wieser算法中的直积形式[2]; 但是该方法仍存在不足之处,由于在构造系数矩阵时需要人工计算或通过编写程序代码自动获得系数矩阵中各个元素的协因数值而使得该方法操作过程复杂,计算所需时间长;

(3) Amiri-Simkooei算法是在Schaffrin-Wieser算法和Mahboub算法的基础上,依据最小二乘理论提出的,因此它更容易被理解,应用起来更为方便、灵活,在高斯- 马尔科夫模型中,参数估值是根据线性迭代算法求得的, 因此应用AmiriSimkooei算法进行参数求解时能够避开需要将公式线性化这一麻烦[4],这样使得计算过程更为简便。而且利用最小二乘理论可以直接得到参数估值的协方差矩阵,而不必重新推导; 但是该算法在表示方面为了与最小二乘方法形式一致,在算法中引入了较多代换量,例如这也使得算法显得冗长。

在实际应用中,需要根据实际情况判断选用哪一种加权总体最小二乘方法更为合适。

摘要：介绍了求解加权总体最小二乘问题的Schaffrin-Wieser算法和Mahboub算法以及Amiri-Simkooei算法,通过MATLAB软件编程进行数值模拟实验比较了三种方法在近似垂直、近似水平和普通直线拟合中的应用以及在二维坐标转换参数求解中的应用,并分析了三种算法的区别与联系。