最小二乘映射

关键词： 参数 优化 设计 性能

最小二乘映射（精选八篇）

最小二乘映射 篇1

工程设计和产品开发越来越注重优化设计, 如何为众多的参数选择最合适的数值, 这将在很大程度上影响产品的质量、性能及成本。优化设计的主要工作是性能的重复分析, 每次系统几何参数或物理参数的改变, 都需要对结构重新进行计算。对于复杂结构多参数系统而言, 用传统的计算方法循环求解不同输入下的输出, 计算费用非常高, 因此, 受制于计算成本, 很难对结构系统进行自适应设计和优化, 或对参数进行稳健的估计。目前解决这类问题的一条途径是对结构进行降阶建模和分析, 如Guyan降阶法[1]、Ritz向量降阶法[2]、正常正交分解法[3]等, 但这些方法在保持原模型的物理特性和提高效率上仍存在一定的局限性, 即使是降阶后, 求解仍比较耗时。

减基法作为一种新的快速计算方法, 在20世纪70年代被提出, 其基本思想是, 当系统由多个参数来描述时, 这些参数的不同组合会使系统方程有不同的解, 而系统在新参数下的解可以用事先设计的样本参数组所对应解的线性组合来得到。近年来, 不少学者对该方法进行了研究, Maday等[4,5]提出了减基法计算的预收敛理论;Rovas[6]把减基法应用于不同种类偏微分方程的求解中;Liu等[7]将减基法扩展到反问题中, 并进行了相关研究。国内对该方法的研究很少, 刘玉秋等[8]将修正减基法用于船舶设计, 雷飞等[9]将减基法用于车身复杂结构大规模问题的快速分析。

上述文献中采用的减基法都是先构造减基空间, 再通过Galerkin映射把原方程投影到低阶方程进行求解。本文在此基础上提出一种新的减基法:在获取减基空间后, 利用最小二乘映射把系统特征矩阵及载荷向量向减基空间投影得到减缩系统, 求解减缩系统并将解映射到原空间得到问题的近似解。本文以某赛车车架为例, 分别通过计算耗时和误差对比, 验证了该减基方法的可靠性和有效性。

1 基于最小二乘映射的减基法理论及算法

1.1基于最小二乘映射减基法理论

对于大型结构静力学问题, 有限元方法通常用偏微分方程弱形式的矩阵表示:

K (μ) U (μ) =F (1)

式中, Κ (μ) 为刚度矩阵;U (μ) 为场变量;F为载荷向量;μ为输入参数, 对于多参数问题, μ为参数向量。

通过拉丁超立方采样, 在参数域Ω中采样N个参数点, 得到集合

SN={μ1, μ2, …, μN} (2)

式中, N为样本个数, 且N≪n;n为原系统总自由度数。

系统在这些参数样本下的响应可表示为{ζ1, ζ2, …, ζN}, ζi为系统平衡方程式 (1) 在μ=μi处的解向量, 记Φ= (ζ1, ζ2, …, ζN) , 同时这些解向量构成系统解空间或减基空间:

WN=span{ζ1, ζ2, …, ζN} (3)

当采样方法科学, 样本点数目选取合理时, 所构造的减基空间具有良好的完备性, 从而保障了降阶投影后原模型的物理特性不发生改变或缺失。因此, 当参数域发生变化时, 新参数下的解能用减基空间中基向量的线性组合表示, 即

$\tilde{\mathit{U}}(\mathit{\mu })=\alpha {}_1(\mathit{\mu })\mathit{\zeta }{}_1+\alpha {}_2(\mathit{\mu })\mathit{\zeta }{}_2+\cdots +\alpha {}_{{\rm N}}(\mathit{\mu })\mathit{\zeta }{}_{{\rm N}}=\mathit{\Phi }\mathit{\alpha }(\mathit{\mu })$ (4)

α (μ) = (α1 (μ) , α2 (μ) , …, αN (μ) ) T (5)

式中, ζi为减基空间的基;αi (μ) 为对应基的权系数。

将近似解方程式 (4) 代入平衡方程式 (1) , 得

K (μ) (Φ α (μ) ) =F (6)

对于静力学问题, 当参数发生变化时, 结构对应的数学模型会发生相应的变化, 即刚度矩阵改变。为了避免参数变化后重新进行前处理以及保证计算的高效性, 将系统刚度矩阵显式分解成与参数相关的系数部分和与参数无关的矩阵部分, 将刚度矩阵表示为

$\mathit{{\rm K}}(\mathit{\mu })={\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i$ (7)

式中, σi (μ) 为与参数相关的函数项;Ki (i=1, 2, …, m) 为刚度矩阵中的不变项;m为刚度矩阵分离项数。

当参数发生变化时, 无需对结构重新进行前处理, 只需要对与参数相关的系数部分进行修正即可快速得到新的系统矩阵。

将式 (7) 代入式 (6) , 得

${\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i(\mathit{\Phi }\mathit{\alpha }(\mathit{\mu }))=\mathit{F}$ (8)

根据矩阵乘法性质, 将式 (8) 写成

$({\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i\mathit{\Phi })\mathit{\alpha }(\mathit{\mu })=\mathit{F}$ (9)

在式 (9) 中, $({\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i\mathit{\Phi })$为n×N维矩阵, α (μ) 为N×1维向量, 且N≪n, 用最小二乘法求解该方程。方程两边同乘以$({\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i\mathit{\Phi }){}^{\text{Τ}}$, 得

$\begin{array}{l}({\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i\mathit{\Phi }){}^{\text{Τ}}({\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i\mathit{\Phi })\mathit{\alpha }(\mathit{\mu })=\\ ({\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i\mathit{\Phi }){}^{\text{Τ}}\mathit{F}(10)\end{array}$

注意到Ki为对称矩阵, 分别展开方程左右可得

$\begin{array}{l}({\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i\mathit{\Phi }){}^{\text{Τ}}({\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i\mathit{\Phi })=\\ {\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m\sigma }}{}_i(\mathit{\mu })\sigma {}_j(\mathit{\mu })\mathit{\Phi }{}^{\text{Τ}}\mathit{{\rm K}}{}_i\mathit{{\rm K}}{}_j\mathit{\Phi }(11)\end{array}$

$({\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{{\rm K}}{}_i\mathit{\Phi }){}^{\text{Τ}}\mathit{F}={\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{\Phi }{}^{\text{Τ}}\mathit{{\rm K}}{}_i\mathit{F}$ (12)

令

Ai j=ΦTKiKjΦ (13)

Bi=ΦTKiF (14)

$\mathit{{\rm K}}{}^{\text{Ν}}(\mathit{\mu })={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m\sigma }}{}_i(\mathit{\mu })\sigma {}_j(\mathit{\mu })\mathit{A}{}_{ij}$ (15)

$\mathit{F}{}^{\text{Ν}}={\displaystyle \sum _{i=1}^m\sigma }{}_i(\mathit{\mu })\mathit{B}{}_i$ (16)

则式 (10) 可写成

KN (μ) α (μ) =FN (17)

式 (13) 和式 (14) 中, Ki为刚度矩阵中不变项, F为载荷常向量, 故Ai j、Bi的值恒定而不受参数变化的影响, 只需在离线阶段计算一次, 其结果可以在在线阶段反复调用。KN (μ) 为减缩系统方程的刚度矩阵, FN为缩减系统方程的力向量。在经过上述最小二乘变换后, 原线性系统由n×n阶降为了N×N阶, 从而使问题求解计算量大大缩减。求解式 (17) 得到减缩系统的响应α (μ) , 将α (μ) 回代式 (4) 即可得到原系统的近似响应。

1.2基于最小二乘映射减基法计算流程

根据减基法计算特点, 将计算分为离线计算与在线计算两部分。

(1) 离线阶段。

①对参数域进行采样并求解系统在样本点的解;②构造减基空间;③分离刚度矩阵, 投影与参数无关的矩阵。

(2) 在线阶段。

①将参数无关矩阵与设计变量集成, 构造新参数下的减缩系统;②求解减缩系统;③以减缩解为系数对解空间基向量进行线性组合, 得到新参数下的近似解。

2 基于最小二乘映射减基法误差定义

由于减基法是通过求解高维问题降阶后的低维系统而快速得到原系统的近似解, 所以在进行降阶的过程中必然存在数值误差。为了验证减基法的计算精度, 在相同参数下比较有限元法计算结果和减缩计算结果, 利用欧几里德范数定义相对误差:

$\mathit{r}{}_i=\parallel \mathit{U}{}_i-\tilde{\mathit{U}}{}_i\parallel /\parallel \mathit{U}{}_i\parallel i=1,2,\cdots ,n$ (18)

式中, Ui为有限元求解的解向量。

计算所有参数下的相对误差, 得到误差向量r= (r1, r2, …, rn) , 将该向量的无穷范数定义为减缩误差

$\epsilon =\parallel \mathit{r}\parallel {}_{\infty }=\underset{1\le i\le n}{{\displaystyle \max }}|r{}_i|$ (19)

3 赛车车架刚度计算算例

以某大学生方程式赛车的车架刚度计算为例, 车架刚度是评价赛车可靠性及安全性的一个重要指标, 它直接影响赛车能承受的最大载荷以及冲击韧性等, 而车架材料及几何尺寸都是影响刚度的重要因素。当进行车架设计时, 需要在各参数范围内选取多组参数组合, 分别进行计算分析以找出最优设计方案。基于传统有限元的计算方法在参数改变时需要重新进行前处理并集成新的特征矩阵, 不仅繁琐而且耗时, 以减基法为基础的快速算法对处理此类多参数变化问题非常有优势。

根据车架各部位承受载荷状况的不同, 该赛车车架的主要受力构件采用直径为25mm的空心圆钢管焊接而成, 次要受力构件及支撑件采用直径为20mm的空心圆钢管焊接而成。取其壁厚作为设计参数, 直径20mm的钢管壁厚t1变化范围为1.2～2.8mm, 直径25mm的钢管壁厚t2变化范围为1.0～2.5mm。根据钢管加工过程包含的杂质元素及热处理工艺的不同[10], 取其弹性模量E变化范围为192～320GPa, 泊松比υ=0.33。计算边界条件如图1所示, 约束A点六个由度, B、C、D三点约束三个平动自由度, 保留三个转动自由度, 在车架顶部四个对称位置分别施加垂直向下的集中载荷, 通过测量车架变形大小来计算车架刚度。用梁单元将车架离散成409个单元, 371个节点, 共2226个自由度。

用拉丁超立方采样法在参数域采取10个样本点, 用减基法把原来2226×2226的系数矩阵缩减为10×10的矩阵, 大幅度提高了求解效率。计算过程用Fortran程序实现, 采用减基法计算, 离线计算时间为27.217s, 在线计算时间为0.00016s, 而用有限元法计算一次的时间为2.018s。当问题规模不大或计算次数不多时, 减基法并不能充分体现出其优势。随着计算次数的增加, 减缩系统的计算时间与有限元计算时间的比较见图2。由图2知, 当计算次数大于某一临界值时, 减缩计算的时间总是小于有限元计算时间, 且随着计算次数的增加, 减缩计算的优势越明显。在需要修改参数进行反复迭代的大规模计算过程中, 采用减缩法能有效提高求解效率。

为验证该减缩算法的稳定性, 在参数域中随机选取50组参数组合:

Ψ={μtest-1, μtest-2, …, μtest-50}

在10个基下分别计算其响应。选取图1中E点作为响应观测点, 考虑位移变化较大的y方向和z方向, 减缩计算与有限元计算结果相对误差分别如图3和图4所示。由图3和图4可知, 所有参数点两个方向的误差均小于3.5×10-4, 说明该方法对参数空间具有良好的适应性, 同时表明通过拉丁超立方采样构造的减缩空间能较好地保持原系统的物理特性。

减缩计算误差与基空间维数之间的关系如图5和图6所示。同样选取y方向和z方向进行观测, 图5和图6表明, 最小二乘映射减基法与基于Galerkin映射的传统减基法相比, 两者的最大减缩误差具有一致收敛特性, 当基空间维数较低时, 最小二乘映射减基法的精度略优于基于Galerkin映射的减基法, 随着基空间维数的增加, 两者的减缩误差都迅速减小并趋于一个稳定值, 且之后再增加基向量的个数不会提高减缩精度。

4 结论

本文提出了一种基于最小二乘映射的减基法。该方法通过样本点求解构造减基空间, 并将系统特征矩阵及载荷向量向减基空间进行最小二乘映射, 从而将求解大规模线性方程组变为解小型线性方程组, 与常规计算方法相比, 节省了求解资源, 有效地提高了计算效率, 同时具有良好的计算精度。赛车车架的算例表明该方法是有效且可靠的。此外, 该方法对大型复杂结构的反复迭代计算优势明显, 故可应用到工程优化设计和反问题求解领域。

摘要：针对机械工程中复杂结构多参数问题, 提出一种新的基于最小二乘映射的减基法。该方法通过在参数域采集样本点, 计算系统在有限个样本点下的响应以构造减基空间, 利用最小二乘映射把原方程向减基空间进行投影得到减缩方程, 在减基空间快速求解该减缩系统, 获得原问题的减缩解, 并把减缩解还原到原空间, 得到问题的近似解。当系统参数发生变化时, 能通过减缩系统快速得到新参数下的响应, 极大地提高了计算效率。最后将该方法用于赛车车架刚度计算, 结果表明方法是有效且可靠的。

关键词：减基法,最小二乘映射,多参数问题,拉丁超立方采样

参考文献

[1] Guyan R J.Reduction of Stiffness and Mass Matrices[J].AIAA Journal, 1965, 3 (2) :380-381.

[2] Wilson E L, Bayo E P.Use of Special Ritz Vectors in Dynamic Substructure Analysis[J].AIAA Journal, 1967, 5 (7) :1944-1954.

[3] Machiels L, Maday Y, OliveiraI B, et al.Output Bounds for Reduced-basis Approximations of Symmetric Positive Definite Eigenvalue Problems[J].Comptes Rendus de l’Académie des Sciences-series I-Mathematics, 2000, 331 (2) :153-158.

[4] Maday Y, Patera A T, Turinici G.Global a Priori Convergence Theory for Reduced-basis Approximation of Single-parameter Symmetric Coercive Elliptic Partial Differential Equations[J].C R Acad. Sci. Paris, S′erie I, 2002, 335 (3) :289-294.

[5] Maday Y, Patera A T, Turinici G.A Priori Convergence Theory for Reduced-basis Approximations of Single-parameter Elliptic Partial Differential Equations[J].Journal of Scientific Computing, 2002, 17 (1/4) :437-446.

[6] Rovas D V.Reduced-basis Output Bound Methods for Parametrized Partial Differential Equations[D].Cambridge, MA:Massachusetts Institute of Technology, 2002.

[7] Liu G R, Zaw Khin, Wang Y Y.Rapid Inverse Parameter Estimation Using Reduced-basis Approximation with Asymptotic Error Estimation[J].Comput. Methods Appl. Mech. Engrg., 2008, 197:3898-3910.

[8] 刘玉秋, 聂武.修正缩减基法在船舶直接设计法中的应用[J].哈尔滨工程大学学报, 1997, 18 (4) :13-18.

[9] 雷飞, 韩旭.车身复杂结构大规模问题的缩减计算[J].中国机械工程, 2009, 20 (17) :2127-2131, 2141.

最小二乘估计教学方案 篇2

教学目标：1、掌握最小二乘法的思想

2、能根据给出的线性回归方程系数公式建立线性回归方程

教学重点：最小二乘法的思想

教学难点：线性回归方程系数公式的应用

教学过程

回顾：上节我们讨论了人的身高与右手一?长之间的线性关系，用了很多种方法刻画这种线性关系，但是这些方法都缺少数学思想依据。

问题1、用什么样的线性关系刻画会更好一些?

想法：保证这条直线与所有点都近(也就是距离最小)。

最小二乘法就是基于这种想法。

问题2、用什么样的方法刻画点与直线的距离会方便有效?

设直线方程为y=a+bx，样本点A(xi，yi)

方法一、点到直线的距离公式

方法二、

显然方法二能有效地表示点A与直线y=a+bx的距离，而且比方法一更方便计算，所以我们用它表示二者之间的接近程度。

问题3、怎样刻画多个点与直线的接近程度?

例如有5个样本点，其坐标分别为(x1，y1)，(x2，y2)，(x3，y3)，(x4，y4)，(x5，y5)与直线y=a+bx的接近程度：

从而我们可以推广到n个样本点：(x1，y1)，(x2，y2)，…(xn，yn)与直线y=a+bx的接近程度：

使得上式达到最小值的直线y=a+bx就是我们所要求的直线，这种方法称为最小二乘法

问题4、怎样使 达到最小值?

先讨论3个样本点的情况

设有3个点(x1，y1)，(x2，y2)，(x3，y3)，则由最小二乘法可知直线y=a+bx与这3个点的接近程度由下面表达式刻画：

整理成为关于a的一元二次函数 ，如下所示：

利用配方法可得

从而当 时，使得函数 达到最小值。

将 代入①式，整理成为关于b的一元二次函数 ，

同样使用配方法可以得到，当

时，使得函数 达到最小值。

从而得到直线y=a+bx的系数a，b，且称直线y=a+bx为这3个样本点的线性回归方程。

用同样的方法我们可以推导出n个点的线性回归方程的系数：

其中

由 我们知道线性回归直线y=a+bx一定过 。

例题与练习

例1 在上一节练习中，从散点图可以看出，某小卖部6天卖出热茶的杯数(y)与当天气温(x)之间是线性相关的`。数据如下表

气温(xi)/oC261813104-1

杯数(yi)/杯202434385064

(1)试用最小二乘法求出线性回归方程。

(2)如果某天的气温是-3 oC，请预测可能会卖出热茶多少杯。

解：(1)先画出其散点图

ixiyixi2xiyi

12620676520

21824324432

31334169442

41038100380

545016200

6-1641-64

合计7023012861910

可以求得

则线性回归方程为

y =57.557-1.648x

(2)当某天的气温是-3 oC时，卖出热茶的杯数估计为：

练习1 已知x，y之间的一组数据如下表，则y与x的线性回归方程y=a+bx必经过点 ( D )

x0123

y1357

(A)(2，2) (B)(1.5，0) (C)(1，2) (D)(1.5，4)

练习2 某连锁经营公司所属5个零售店某月的销售额和利润额资料如下表：

商店名称ABCDE

销售额(x)/千万元35679

利润额(y)/百万元23345

(1)画出销售额和利润额的散点图;

(2)若销售额和利润额具有相关关系，计算利润额y对销售额x的回归直线方程。

解：(1)

(2)数据如下表：

ixiyixi2xiyi

13296

2532515

3633618

4744928

5958145

合计301712

可以求得b=0.5，a=0.4

线性回归方程为：

小结

1、最小二乘法的思想

2、线性回归方程的系数：

最小移动二乘抠图 篇3

抠图算法可分为3类:基于采样的方法、基于传播的方法、采样和传播结合的方法。其中，基于传播的方法有着举足轻重的作用，Levin等根据局部线性模型[1]推导出的拉氏矩阵给出邻域像素的alpha值间的线性关系，被广泛应用在抠图算法中。拉氏矩阵有其局限性，拉氏矩阵表示空间邻域内像素间的关系，但不能体现非邻域间像素间的关系。拉氏矩阵建立在空间连续的假设基础上，在某些前景和背景分量突变的区域，拉氏矩阵难以得到理想的效果。

本文提出一种改进的拉氏矩阵计算方法，使用最小移动二乘法替代最小二乘法推导出移动拉氏矩阵。相对于最小二乘法，移动最小二乘法求解的线性条件更为准确。笔者使用KNN邻域替代空间邻域，使得拉氏矩阵可以反映非邻域间像素的alpha值的关系。试验结果表明本文方法是有效的。

1 相关工作

早期的抠图算法例如贝叶斯抠图基于颜色空间上的采样，对任意未知像素，基于采样的方法在三分图的已知前景和背景区域中寻找最优的前景和背景对匹配当前像素，像素对应的alpha值与像素的空间位置无关，该类方法不能保持alpha图的空间连续性。在前景和背景颜色混杂时，未知像素可能和多个的前景背景颜色匹配，导致生成错误的alpha图。近年来发展的基于采样的方法如共享式采样[2]假定在小邻域内前背景分量一致。共享式采样为了保证结果的光滑性，并提高正确率，改进了信任度计算方法，并对结果进行多次平滑操作，

Rajan等人[3]提出了颜色及纹理采样方法，不仅考虑到了像素的信息，还考虑纹理的信息，从而使得采样更为可靠鲁棒。但基于采样的方法不能很好的保证alpha图在空间上的连续性，也不能有效的利用空间连续性计算alpha图，不同于基于采样的方法基于传播的方法充分考虑相邻像素间的空间关系，基于传播的方法构造像素在alpha图上的关系，alpha值的求解类似于alpha值从已知区域向未知区域传播的过程，例如如泊松抠图、随机游走抠图[4]、闭形式抠图[1]等。

由于抠图问题是典型的病态问题，因此需要alpha图满足一定的假定条件才可求出唯一解，例如泊松抠图假定像素对应的前景分量和背景分量是光滑的，在此基础上闭形式算法使用，并由此给出相邻像素间alpha图的变化信息。闭形式抠图假定像素对应的前景分量和背景分量在局部小邻域内是不变，从而根据局部线性假设得到alpha值在局部邻域空间上的线性关系，从而以拉氏矩阵的形式给出alpha图在空间上的关系。旧的传播方法如泊松抠图随机游走抠图需要调整局部或全局参数，相比之下，闭形式方法不需要调整参数且效果更佳，因此被广泛的应用的现有的抠图算法中。

闭形式抠图方法需要求解大型稀疏线性方程，因此效率不高。闭形式方法由于基于小邻域下的传播模型，在某些前背景分量无法传播到的区域无法取得良好的效果。何凯明等提出了大核抠图方法[5]，何等发现拉氏矩阵在大核情况下运用共轭梯度法求解收敛更快，而且大邻域可以覆盖更多的前景和背景区域，从而可以传播更多的前景和背景分量，从而在某些有洞的区域得到较好的效果。大核方法假设内像素的alpha值在相对较大的邻域内保持线性关系，其假定条件对比于小核更不易成立，因此在复杂纹理处效果不佳。Zheng等人[6]提出了基于学习的抠图方法，将局部线性模型扩展为局部非线性模型，从而在非线性关系下也能取得较好的效果。

由于闭形式方法获取邻域上像素的alpha值关系，本质上属于邻域方法，不能获取像素在非邻域上的关系。近几年的基于传播的方法侧重于从非邻域上获取alpha值的关系，Lee等人[7]提出了非邻域抠图方法，从非邻域方法构造拉氏矩阵。KNN(K Nearest Neighbor)方法[8,9]使用KNN邻域替代空间邻域。Chen等则提出了邻域和非邻域结合[10]的方法，该方法假定alpha图能保持图像空间上的流形结构，从而得出alpha图在图像流形上的非邻域结构关系，再和邻域上的拉氏矩阵结合，从邻域和非邻域二个方面得出新的传播模型。

拉氏矩阵在图像编辑中有大量的应用[11,12]，从现有的抠图方法来看，闭形式对应的拉氏矩阵仍然在现有抠图方法中起到重要的作用，而抠图方法在一系列应用中有着基础性的应用价值[13,14]。本文使用KNN邻域替代空间邻域，从而可以获取非邻域像素在alpha图上的线性关系，并用移动最小二乘替代最小二乘，从而计算出移动拉氏矩阵，并得到alpha图，实验结果表明移动拉氏矩阵更为有效。

2 最小移动二乘抠图

闭形式方法基于局部线性假设，表示如下

当局部邻域内假设条件不成立时，特别是邻域比较大且纹理复杂的情况下效果不佳。假设在邻域内alpha值满足线性条件，不同于闭形式方法[1]使用最小二乘法求解局部线性关系，在窗口wi内使用移动最小二乘法求解局部线性关系，表示如下

与闭形式抠图不同处在于:在最小化式(2)和式(3)中增加了权值ω，移动最小二乘在距离当前像素越远的地方权值ω越小，因此移动最小二乘法能求解出更准确的局部线性关系，比最小二乘法求解的线性关系更为有效。ωi是邻域wk中的权值。式(2)和式(3)可以表示为以下矩阵的形式

对于每个邻域wk,Gk定义为矩阵。Gk每行包括向量(Ii,1)。Wk是每行向量对应的权值向量。G'k为Gk的Wk加权，对应的每行向量表示为(Wk·Ii,Wk)。αk是邻域内所有像素对应的alpha值组成的向量。

系数ak,bk求解如下

令可表示为

式中:δi,j是Kronecker delta函数;μk和σ2分别是小窗口wk内的基于Wk的加权均值和方差;是窗口内像素的个数。

2.1 彩色模型下的移动最小二乘抠图

彩色模型下类似于闭形式算法，用式(9)表示彩色图像各通道间的线性关系

式中:c为彩色图像的通道数，在考虑各个通道信息后，式(2)和式(3)转化为

对式(10)进行化简后，解得彩色模型下移动拉氏矩阵如下

式中:I为小邻域内所有像素对应3×1颜色向量组成的矩阵;μk为I的Wk加权平均;Σk是I在Wk加权下的协方差矩阵。

2.2 KNN邻域

由于拉氏矩阵不能反映像素的非邻域关系，在此借鉴KNN抠图引入KNN邻域，将拉氏矩阵中的空间邻域扩展到KNN邻域，KNN空间的点由(R,G,B,X,Y)五维共同决定。使用KD-TREE实现KNN邻域的高效查找。由于KNN邻域反映空间非邻域上像素间的关系。因此本文方法结合了非邻域抠图的优点。

3 移动二乘抠图中的大核求解

由于移动抠图算法中，设核的大小为r，图像像素个数为imagesize，存储拉氏矩阵L所需要的空间复杂度为imagesize×r2，计算空间复杂度随着核的增大而急剧增大。借鉴大核方法中[5]的计算技巧，使用改进的共轭梯度法求解alpha值。

对于方程Lx=b，共轭梯度法的关键在于构造共轭向量p，并求其对应的残差。共轭梯度法可以用迭代方法求解。在每次迭代过程中，新共轭向量由求解如下

共轭方向的系数求解如下

新的x值与残差用求解如下

在共轭梯度求解过程中关键的步骤在于求解向量Lp，直接求解L的空间复杂度过大，但Lp的维数为imagesize，因此需要避免直接求解L，而直接用下式求解Lp向量中点i对应的元素qi

式中:Wk是像素k对应的邻域;是邻域的大小;i是包围像素k邻域Wk中的一个像素;qi为q向量的第i个元素;I为像素i对应的3维向量，表示R,G,B三个通道;pi为共轭向量中像素i对应的元素;μk是3维向量，为邻域Wk中Ii向量的均值;为邻域Wk中元素i对应的共轭向量pi的均值;ak*是像素k的对应的3维向量;bk*为像素k对应的标量。

(Lp)i计算公式的正确性由以下定理保证:

定理1:式(17)计算得出的(Lp)i与利用式(12)计算出的(Lp)i是等价的。

证明:

令q=Lp，由于q与p为线性关系，因此只需要证明式(21)

将式(20)代入式(17)消除bk*，可得到

此外，有

根据式(18)，并对pj做偏导，可得

将式(23)和式(24)代入式(22)，可得到

式(25)即式(12)中对应的拉氏矩阵L。

4 实验结果

当前的抠图方法中，大多使用了基于传播和基于采样结合的方式提高算法的准确性，在此将最常见的几种传播方法闭形式方法、KNN方法、基于学习的方法、大核方法进行了比较。

从图1可知，闭形式算法对细节处理较好，由于其传播模型只考虑到了空间小邻域间的信息，而Doll毛发附近的字符间有间隔，是空间上的非邻域关系，因此alpha值不能在字符间顺利传播，也无法干净的抠除这些字符。在plant图中，由于植物树叶间的空洞与三分图中的背景分量并没有空间上的邻域关系，因为闭形式方法同样难以将背景分量传播至植物中的空洞中。

在大核方法中，由于传播模型在相对较大的邻域内进行传播。在大核情况下，毛发附近的字符是邻域关系，因此算法成功的抠除了毛发附近的字符，但由于植物树叶间的空洞与背景空间距离较远，所以大核仍然无法解决植物中的空洞问题。此外，局部线性假设不易在大核时成立，因此算法对复杂纹理处理的效果不好。

由于KNN方法建立在非邻域的基础上，因此KNN方法可以在非邻域上进行alpha值的传播。KNN方法因此可以在植物的空洞以及Doll毛发附近的字符处取得良好的效果。

相对于闭形式方法，KNN方法对应的拉氏矩阵建立在全局统一的参数基础上，因此细节处理不好。在Doll图片中，KNN方法使得毛发周围较为模糊，不能较好地提取毛发，毛发左边的英文字母也没有清除干净。在plant图中算法在叶子周边没有抠除干净。在Plastic bag图片中，KNN方法同样在绳子附近留下大量噪声。而本文方法在Doll周边的毛发处理较为干净，特别是英文字符全部被干净的清除，在Plastic bag中，绳子附近处理也较为干净。Tree图中在有洞的区域同样得到较好的效果。从表1的参数可以看到，net图中的方法效果明显好于其他方法，由于net图中的未知区域非常大，而本文方法提供了更为准确的拉氏矩阵，因此得出了较好的结果。

摘要：交互式抠图技术在有限的用户交互下抠取图像的前景,被广泛地应用在图像及视频编辑、三维重建等领域中,有极高的应用价值。近年来的抠图技术中,拉氏矩阵给出alpha图上像素间的线性关系,对alpha图的估计起到了重要作用。提出了一种改进拉普拉斯矩阵的方法,使用移动最小二乘法替代最小二乘法,结合最近邻(KNN)方法给出移动拉氏矩阵,并使用移动拉氏矩阵计算alpha图。实验结果证明了移动拉氏矩阵的有效性。

多维正交整体最小二乘应用研究 篇4

在利用最小平差对线性模型的参数数据求解过程中, 通常假定数据矩阵与结构矩阵不存在误差, 如果实际测量的观测数据中, 自变量和因变量都不可避免的含有误差, 那么最小二乘估计从统计观点上看就不是最优解。此时就得利用Golub等人提出的整体最小二乘来对此问题进行求解[1], 并导出整体最小二乘求解的奇异值分解法[2]。但是, 由于SVD法的复杂性限制了TLS的广泛应用, Schaffrin等人提出了整体最小二乘法的递推算法[3,4], 同时鲁铁定对整体最小二乘的基本算法也进行了变换[5]。以往对整体最小二乘模型进行研究是基于一元线性回归为准则, 但是当由一元延伸到多元的时候, 基于SVD分解的整体最小二乘并不是完全在VTV=min的情况下进行解算, 同时解算的方程只考虑前n-1行, 而没有完全顾及整个方程, 从而造成了解算结果误差的偏执。而基于鲁铁定提出的ax+by=1模型, 为了强调解算的方便, 忽略了过多的误差因子δV, 解算过程仍然依据VTV=min的原则进行的, 因此解算的结果明显出现了扭曲的现象, 并不是真正意义上按照VTV=min的原则进行解算的。本文基于此, 对于一元线性回归, 避开引入误差变量, 利用通过二维数据点的中心直接寻求二维数据点到直线的垂直距离最短的原则来构造正交整体最小二乘模型。对于此模型进行延伸, 使其通过多维数据点的中心, 寻找多维数据点到曲面的垂直距离最短, 从而形成了多维正交整体最小二乘模型。此模型不仅解决了整体最小二乘应用广泛性的问题, 而且也遵循了VTV=min的原则。

1 多维整体最小二乘法

目前, 平差处理用的最广泛的是最小二乘。但是最小二乘只考虑了因变量误差, 而在测绘数据处理实践中, 其自变量存在误差的情况较多, 此时就有必要研究系数矩阵和观测矩阵向量都存在误差的数据处理理论, 整体最小二乘将有助于解决此类问题。假设整体最小二乘的误差方程式为:

对式 (1) 进行变形得:

式 (2) 中, B=[1 x y z];D=[0 VxVyVz];Z=[a b c-1]T。

利用Lagrange乘数法求解, 目标函数为:

根据求极值定理确定BTBZ=KZ

Lagrang乘数K应该选择为矩阵BTB的最小的特征值λm+1, 再由BTBZ=KZ的前3个式中解得参数a, b, c的值。

对于模型:a (x+Vx) +b (y+Vy) +c (z+Vz) =1, 令a=a0+δa, b=b0+δb。

整理成:FV+BδB-W=0, 利用间接平差求得参数:

式 (4) 中,

上式模型解算的过程中忽略了δa Vx, δb Vy, δc Vz, 所以极小化原则中的Vx, Vy, Vz并不是实际系数矩阵与观测矩阵的误差。但解算依然依据VxTVx+VyTVyVzTVz=min的原则, 所以此时解算的结果势必会造成误差的偏执。

2 多维正交整体最小二乘模型的构造

2.1 模型构造的基础

假设有n个数据观测点 (x1, y1) , (x2, y2) , …, (xn, yn) , 对这些数据点进行拟合一条直线, 为了能够拟合一条最佳的直线, 则需要考虑这n个数据点到拟合直线的距离平方和最小, 即在顾及x、y方向误差的过程中, 寻找点到直线距离平方和最短的直线方程, 如图一所示。

首先, 令直线方程为ax+by-c=0。假设此直线方程通过点 (x0, y0) , 则c=ax0+by0。因此将直线方程写作:

由几何及微积分可知, 点 (xi, yi) 到直线a (xx0) +b (y-y0) =0的距离d为:

根据文献[2], 拟合直线必须通过n个数据点的中心, 才能使偏差d2最小, 则基本直线方程变为:

式 (7) 中,

2.2 正交整体最小二乘模型的构造

在所有观测数据点到式 (8) 直线距离之和最小的情况下, 求出a、b的值, 便是整体最小二乘的最佳解。为了求拟合直线的法向量[a, b]T, 考虑如何使偏差d2最小, 可将d2写成2×1单位向量t= (a2+b2) -1/2[a, b]T与n×2矩阵M的乘积, 即:

式 (9) 中,

根据n×2矩阵M的形式可知2×2矩阵MTM是一个实的对称半正定矩阵。对其进行特征值分解:

式 (11) 中, σ1>σ2。于是:

由于向量范数相对于正交矩阵是不变的, 即|UTt|2=||t||2, 同时t= (a2+b2) -1/2[a, b]T为Euclidean范数等于1的向量, 则式 (12) 简化为:

另一方面,

由于t, u1, u2的Euclidean范数均等于1, 且特征向量u1和u2相互正交, 若t=u2, 则:

由式 (14) 、 (15) 可知, 当t=u2时, 距离平方和取最小值σ22。

向量[a, b, c]T解算结果为MTM最小特征值所对应的特征向量, 其中, 假设观测数据点[x, y, z]为独立同精度观测, 则可令

根据观测数据的标准差公式及协因数传播率可得单位权中误差:

3 实例分析

以文献[6]中数据作为样本, 如表一所示。选取坝体温度和水位压力作为自变量x、y, 大坝水平位移值为观测量z。分别利用本文方法与最小二乘和其他模型的整体最小二乘进行参数解算, 将模型解算结果统一规划为z=a+bx+cy形式, 并对残差及单位权中误差进行对比分析, 结果如表二所示。

从表二可以看出, 本文基于正交多维整体最小二乘拟合残差及单位权中误差较优于最小二乘和另外两种模型的整体最小二乘。

4 结束语

本文针对奇异值分解解算整体最小二乘的复杂性及同时顾及自变量、因变量解算过程中产生误差变异的情形, 将正交距离引入到整体最小二乘的解算过程中。此方法没有直接考虑自变量、因变量误差, 解算模型过程中不存在Lagrange乘数K的确定、方程组解算的复杂性及影响解算过程而忽略δV的问题。因此不仅对二维还是多维都具有较高的拟合精度及拟合稳定性, 具有较广泛的应用。

本文只是对所有观测数据点拟合精度进行考虑, 而没有顾及到数据点本身粗差的影响, 因此, 如果观测数据点中存在粗差, 在没有对其剔除的情况下, 依然用此方法对其进行拟合及误差分析, 将会产生变异的误差曲线或误差曲面。因此, 怎样将观测数据点粗差剔除的方法与本文正交整体最小二乘方法结合起来、实现无缝链接, 将有待进一步探讨。

参考文献

[1]Golub G H, van Loan C F.An Analysis of theTotal Least Squares Problem[J].SIAM Journal on Numerical Analysis, 1980, 17 (06) :8832-893.

[2]Beltrami E.Sulle funzioni bilineari, Giomaledi Mathematiche ad Uso Studenti Delle Uninersita.1873, 11:98-106.An English translation by D Boley isavailable as University of Minnesota Department ofComputer Science, Technical Report 90-37, 1990.

[3]Schaffrin B.A Note on Constrained Total LeastSquares Estimation[J].Linear Algebra and Its Applications, 2006, (417) :245-258.

[4]Schaffrin B, Wieser A.On Weighted TotalLeast Squares Adjustment for Linear Regression[J].Geod, 2008, (82) :415-421.

[5]鲁铁定, 陶本藻, 周世健.基于整体最小二乘法的线性回归建模和解法[J].武汉大学学报 (信息科学版) , 2008, 33 (05) :504-507.

基于最小二乘原理的影子定位分析 篇5

日照形成的影子会随着太阳光线的变化而发生变化, 而太阳对地球的照射会因照射时间、照射地点的经纬度的不同而不同。通过分析太阳照射下直杆影子的长度和位置变化, 可以得出太阳的运动轨迹, 进而确定出直杆所处的地理位置。

2太阳角度相关量说明和计算公式[1]

年积日N:某日期在年内的顺序号;

太阳赤纬δ[2]:地球赤道平面与太阳和地球中心连线之间的夹角:

太阳时角ω:某一时刻, 日心到地心连线所在子午圈与观察者所在子午圈之间的夹角

其中ST为当地时间, t为北京时间, λ为时差, φ为当地经度。

太阳高度角β:太阳光入射方向和地平面之间的夹角;

太阳方位角α:太阳光入射方向和当地经线之间的夹角;

3影子长度变化模型

3.1影子长度计算

影子的推移随着时间发生着变化, 根据文献[3]我们得到高度为h的直杆的投影端点坐标 (以正北为y轴, 正东为x轴建立直角坐标系) 的计算公式:

式中上午时刻取负号, 下午时刻取正号, β>0即影子轨迹的计算范围为某地某天太阳的日出和日落时间点之间的曲线。

假设直杆所处的坐标为 (0, 0) , 根据公式 (5) 求出影子端点坐标即可得到影子的

3.2影子长度变化规律

从上述所建立的模型中可以看出, 影子的长度直接与直杆的高度h和太阳高度角有关, 可得到以下结论:

(1) 对于同一地点, 同一时刻, 杆子越长, 影子的长度越大; (2) 对于同一直杆, 同一纬度地区, 影子长度变化轨迹相同, 且越往东, 影子出现最短的时刻越早; (3) 对于同一直杆, 同一地点, 同一日期, 时刻越靠近正午, 影子长度越小。 (4) 对于同一直杆, 同一时刻, 在太阳直射点处的杆子影子长度最小, 影子长度随着杆子位置向南北两侧递减。

这都符合日常客观规律。

4影子定位模型[4]

很多情况下, 我们想知道某次影子测量的测量地点, 但是我们只知道测量时间和影子的顶点坐标, 且不明确影子顶点坐标系的建立规则。因此, 我们建立一个影子定位模型来确定影子的测量地点。

4.1基于最小二乘的影子定位模型

在这个模型中, 已知:某直角坐标系下影子n个顶点坐标 (x[n], y[n]) ;影子顶点坐标测量时间间隔为t0, 影子开始测量时间为北京时间t1, 影子结束测量时间为北京时间t2。

由于每隔t0时段, 影子转过的角度是确定的, 即

利用公式 (2) , (3) , (4) 可得出理论上该夹角的值:

我们利用最小二乘法求得一个最优解:

4.2实例说明

4.2.1实例介绍

有一组数据[5]坐标系以直杆底端为原点, 水平地面为XY平面。直杆垂直于地面。测量日期:2015年4月18日。

4.2.2约束范围确定

在北京时间14:42到15:42测量地点处于白天。由于时区之间的差异得到经度范围在西经0.5度到东经169度之间, 且在北京时间14:42到15:42测量地点的影子长度是在逐渐增加的, 可以得出, 此时当地时间已经过了正午12点, 则当当地时间为正午12点时, 北京时间t<14.7, 根据公式 (2) 可得经度大于东经79.5度, 则测量地点:E79.5°<φ<E169°

4.2.3测量地点确定

利用基于最小二乘的影子定位模型, 得:

根据地图查找, 该点处于南海区域。考虑到测量的实际操作性, 以及模型存在的误差, 我们将范围扩大到了距离该点最近的沿岸区域。将直杆所在地点大致定在中国海南乐东黎族自治县的莺歌海镇, 立国镇等地区 (N18.52°, E108.69°) 或越南顺化峰牙-已傍国家公园 (N17.34°, E106.76°) 。

5结束语

在影子长度变化模型中, 建立以正北方向为X轴, 正东方向为Y轴的平面直角坐标系, 利用太阳高度角β和太阳方位角α, 确立了长为h的直杆的影子的顶点坐标, 并进一步确定了直杆长度变化计算公式。

在直杆影子定位模型中, 由于如何建立XY直角坐标平面未知, 不能直接以arctany/x的值来当作实际测得的α值, 因此, 通过利用两个测量时间之间的影子夹角建立了基于最小二乘的影子定位模型。

参考文献

[1]刘海波, 王建芳, 于海芹, 等.太阳能工程中几种相关角度的计算及应用[J].创新与技术.

[2]于贺军.气象用太阳赤纬和时差计算方法研究[J].气象水文海洋仪器, 2006 (3) :50-53.

[3]夏利江.青藏铁路旱桥桥面遮阳对桥下及周边冻土太阳辐射影响[J].地球科学进展, 2014, 29 (3) :380-387.

[3]汽车减振器的发展与现状[D].吉林:吉林大学珠海学院.

[4]武琳.基于太阳阴影轨迹的经纬度估计技术研究[J].Computer Vision and Image Understanding, 2010, 114 (8) :915-927.

最小二乘映射 篇6

Hammerstein模型是由静态非线性模块和动态线性模块相互连接构成。最常用的Hammerstein系统辨识方法有3种:过参数化辨识方法[1]、迭代方法[2]和关键变量分离方法[3]。过参数化辨识方法的参数向量包含了原系统静态非线性模块参数与线性动态模块参数的乘积项,这使得过参数模型的参数向量维数(即参数数目)大大增加,导致辨识算法的计算量也相应增加。

针对Hammerstein OEMA 模型,将关键变量分离原理[3]与辅助模型思想[4,5]相结合,提出基于关键变量分离的辅助模型递推增广最小二乘辨识方法。关键变量分离原理是分离出线性环节的关键变量,将非线性环节表达式代入到其中去,从而使整个系统的输出表示为所有待估参数的回归方程形式,不存在冗余参数,算法计算量小。

1 Hammerstein OEMA 模型描述

如图1 所示Hammerstein OEMA 模型表示如下

$y(t)=\frac{B(z)}{A(z)}x(t)+D(z)v(t)(1)$

x(t)=f[u(t)] (2)

其中u(t)和y(t) 分别为系统的输入和输出; x(t) 、r(t)和w(t) 为未知中间变量; f[·]是非线性函数;v(t) 为零均值、不相关随机白噪声(不可测)。A(z),B(z)和D(z)均为单位后移算子z-1 的多项式[z-1y(t)=y(t-1)] ,且

A(z)=1+a1z-1+a2z-2+…+anaz-na,

B(z)=b0+b1z-1+b2z-2+…bnbz-nb,

D(z)=1+d1z-1+d2z-2+…+dndz-nd。

假设t≤0 时,u(t)=0, y(t)=0, v(t)=0, 且阶次na,nb,nd已知。非线性环节x(t) =f[u(t)]是一个分段线性函数。当u(t) ≥ 0时,x(t)=m1u(t);当u(t)≤0时, x(t)=m2u(t)。引入开关函数

$h(t)=h[u(t)]=\{\begin{array}{l}1,u(t)＞0‚\\ 0‚u(t)＜0。\end{array}$

则函数x(t) 可记为

x(t)=m2u(t)+(m1-m2)h(t)u(t) (3)

任意非零常数α,对于方程(1)而言,$\left[\alpha x(t),\frac{B(z)}{\alpha }\right]$与[x(t),B(z)]产生的效果完全相同。为了得到唯一解,令b0=1。

本文的目标是: 使用所采集的数据{u(t),y(t)},借助于关键变量分离原理和辅助模型思想,提出新的算法,估计系统的参数mi,ai,bi,di。

2 使用关键变量分离的辅助模型辨识方法

先引入符号:“A=:X”或“X:=A”表示“A 记作(定义为) X”之意(因为符号$\begin{array}{l}\Delta \\ =\end{array}$没有左右之分,含义模糊);上标T表示矩阵转置。

参照图1, 定义中间变量:

r(t):$=\frac{B(z)}{A(z)}x(t)(4)$

w(t):=D(z)v(t) (5)

由式(1)可得

y(t)=r(t)+w(t) (6)

由式(4)可得

r(t)=x(t)+b1x(t-1)+…+bnbx(t-nb)-a1r(t-1)-…-anar(t-na)

参见文献[3],只将式(3)代入上式中的被分离出来的关键变量x(t) 中,其余的(没有被分离出来的)关键变量x(t-i),(i=1,2,…,nb)保持不变。于是有

r(t)=m2u(t)+(m1-m2)h(t)u(t)-a1r(t-1)-…-anar(t-na)+b1x(t-1)+…+bnbx(t-nb) (7)

定义如下信息向量

φs(t):=[u(t),h(t)u(t),-r(t-1),…,-r(t-na),

x(t-1),…,x(t-nb)]T∈Rna+nb+2,

φn(t):=[v(t-1),v(t-2),…,v(t-nd)]T∈Rnd。以及参数向量

θs:=[m2,m1-m2,a1,…,ana,b,…,bnb]T∈Rna+nb+2。

θn:=[d1,d2,…,dny]T∈Rni

那么式(7)和式(5)写成向量形式分别为

r(t)=φsT(t)θs (8)

w(t)=φnT(t)θn+v(t) (9)

把式(8)和式(9)代入式(6)可得

y(t)=φsT(t)θs+φnT(t)θn+v(t)=:φT(t)θ+v(t) (10)

式(10)中

注意到信息向量φ(t) 中包含未知中间变量r(t-i), x(t-i)和不可测噪声项v(t-i)。 为此在辅助模型辨识思想的基础上[4,5],构造一个辅助模型ra(t)=φTa(t)θa,rn(t)是辅助模型的输出。φ(t):=φs(t)φn(t)∈Rna+nb+nd+2,θ:=θsθn∈Rna+nb+nd+2

信息向量φ(t)中φs(t)里的未知项r(t-i)用辅助模型的输出ra(t-i)代替,不可测中间变量(未分离出来的关键变量)x(t-i)用其估计值$\widehat{x}(t-i)$代替,作了代替后的φs(t)记为$\widehat{\phi }$s(t),即

$\begin{array}{l}\widehat{\phi }{}_{\text{s}}(t)=[u(t),h(t)u(t),-r{}_a(t-1),-r{}_a(t-2),\\ \cdots ,-r{}_a(t-n{}_a),\widehat{x}(t-1),\widehat{x}(t-2),\cdots ,\\ \widehat{x}(t-n{}_b]{}^{\text{Τ}}(11)\end{array}$

信息向量φ(t)中φn(t)里的未知项v(t-i)用估计残差$\widehat{v}(t-i)$代替,替代后的φn(t)记为$\widehat{\phi }$n(t),即

$\widehat{\phi }{}_n(t)=[\widehat{v}(t-1),\widehat{v}(t-2),\cdots ‚\widehat{v}(t-n{}_d)]{}^{\text{Τ}}(12)$

选择$\widehat{\phi }$s(t)作为辅助模型的信息向量φa(t),估计$\widehat{\theta }$s(t)作为辅助模型的参数向量θa,则

$r{}_a(t)=\widehat{\phi }{}_s^{\text{Τ}}(t)\widehat{\theta }{}_s(t)(13)$

用$\widehat{\phi }$s(t)和$\widehat{\phi }$n(t)代替φ(t)中未知的φs(t)和φn(t),代替后记作为

$\widehat{\phi }$(t):$=[u(t),h(t)u(t),-r{}_a(t-1),\cdots ,-r{}_a(t-n{}_a),\widehat{x}(t-1),\cdots ,\widehat{x}(t-n{}_b),\widehat{v}(t-1),\cdots ‚\widehat{v}(t-n{}_d)]{}^{\text{Τ}}(14)$

用估计值$\widehat{\phi }$(t)和$\widehat{\theta }$(t)分别代替式(10)中φ(t)和θ(t),经变换得到估计残差的计算式,

$\widehat{v}(t)=y(t)-\widehat{\phi }{}^{\text{Τ}}(t)\widehat{\theta }(t)(15)$

作了上述准备, 极小化准则函数

$J(\theta )={\displaystyle \sum _{i=1}^t[}y(i)-\widehat{\phi }{}^{\text{Τ}}(i)\theta ]{}^2$。仿照标准递推最小二乘法的推导过程,可得到Hammerstein OEMA 模型基于关键变量分离的辅助模型递推增广最小二乘辨识方法,

$\widehat{\theta }(t)=\widehat{\theta }(t-1)+L(t)[y(t)-\widehat{\phi }{}^{\text{Τ}}(t)\widehat{\theta }(t-1)](16)$

$L(t)={\rm P}(t-1)\widehat{\phi }(t)[1+\widehat{\phi }{}^{\text{Τ}}(t){\rm P}(t-1)\widehat{\phi }(t)]{}^{-1}(17)$

${\rm P}(t)=[{\rm I}-L(t)\widehat{\phi }{}^{\text{Τ}}(t)]{\rm P}(t-1),{\rm P}(0)=p{}_0{\rm I}(18)$

$\widehat{x}(t)=\widehat{m}{}_2(t)u(t)+[\widehat{m}{}_1(t)-\widehat{m}{}_2(t)]h(t)u(t)(19)$

$\begin{array}{l}\widehat{\theta }{}_{\text{s}}(t)=[\widehat{m}{}_2(t),\widehat{m}{}_1(t)-\widehat{m}{}_2(t),\widehat{a}{}_1(t),\cdots ,\widehat{a}{}_{n{}_a}(t),\\ \widehat{b}{}_1(t),\cdots ,\widehat{b}{}_{n{}_b}(t)]{}^{\text{Τ}}(20)\end{array}$

$\widehat{\theta }{}_{\text{n}}(t)=[\widehat{d}{}_1(t),\widehat{d}{}_2(t),\cdots ,\widehat{d}{}_{\text{n}{}_d}(t)]{}^{\text{Τ}}(21)$

3 仿真例子

考虑如下Hammerstein OEMA 模型,

$\begin{array}{l}y(t)=\frac{B(z)}{A(z)}x(t)+D(z)v(t),x(t)=\{\begin{array}{l}m{}_{1}u(t),u(t)＞0;\\ m{}_{2}u(t),u(t)＜0。\end{array}\\ m{}_1=1.50,m{}_2=3,\\ A(z)=1+a{}_{1}z{}^{-1}+a{}_{2}z{}^{-2}=1-1.60z{}^{-1}+0.80z{}^{-2},\\ B(z)=b{}_0+b{}_{1}z{}^{-1}+b{}_{2}z{}^{-2}=1+0.20z{}^{-1}-0.30z{}^{-2},\\ D(z)=1+d{}_{1}z{}^{-1}=1+0.90z{}^{-1}。\end{array}$

θ=[m2,m1-m2,a1,a2,b1,b2,d1]T=[3.00,-1.50,-1.60,0.80,0.20,-0.30,0.90]T。

仿真时,系统输入{u(t)}采用零均值、单位方差、不相关可测随机信号,{v(t)} 采用零均值方差为σ2的随机白噪声信号。当σ2=0.102和σ2=0.502时,采用本文提出的算法估计这个系统的参数,参数估计误差$\delta =\parallel \widehat{\theta }(t)-\theta \parallel /\parallel \theta \parallel$随数据长度t变化的曲线如图2 所示。

由图2 可以看出:参数估计误差随着数据长度的增加不断减小;噪信比越小,参数估计收敛于真值的速度越快,且估计误差也越小。

4 结论

研究了Hammerstein模型关键变量分离的辅助模型递推增广最小二乘辨识方法。算法原理具有特色,容易实现。仿真结果证明了算法的有效性。

摘要：针对Hammerstein输出误差自回归(OEMA)模型,将关键变量分离原理与辅助模型辨识思想相结合,提出了基于关键变量分离的辅助模型递推增广最小二乘辨识方法。该方法能获得系统参数估计和噪声参数估计,且能实现在线辨识。

关键词：Hammerstein模型,关键变量分离原理,辅助模型,递推辨识

参考文献

[1] Ding F,Chen T.Identification of Hammerstein nonlinear ARMAXsystems.Automatica,2005;41(9):1479—1489

[2] Liu Y,Bai E W.Iterative identification of Hammerstein systems.Automatica,2007;43(2):346—354

[3] V r s J.Recursive identification of Hammerstein systems with discon-tinuous nonlinearities containing dead-zones.IEEE Transactions onAutomatic Control,2003;48(12):2203—2206

[4]王冬青,丁锋.基于辅助模型的多新息广义增广随机梯度算法.控制与决策,2008;23(9):999—1003,1010

最小二乘影像匹配算法的实现与研究 篇7

现阶段常用的影像匹配的方法有:最小二乘影像匹配、基于特征影像匹配、多点匹配、松弛法影像匹配、附有几何约束的多影像匹配等[1]。其中, 1986年由德国Ackermann教授提出的最小二乘影像匹配算法具有灵活、可靠、高精度等特点[2], 在数字摄影测量中得到了广泛的应用, 有很高的研究价值。

1 算法原理

最小二乘影像匹配算法的基本原理是将辐射畸变和几何畸变两大影像灰度系统变形参数引入到影像匹配的模型中, 按最小二乘原则求解变形参数, 利用变形参数在初匹配的匹配结果基础上进一步提高了匹配精度[3]。

该算法的具体步骤为[2]:

(1) 几何变形改正。根据几何变形改正参数a0, a1, a2, b0, b1, b2将左方影像窗口的影像坐标变换至右方影像阵列:

(2) 重采样。由于换算所得之坐标x2, y2一般不可能是右方影像阵列中的整数行列号, 因此需要进行重采样。

(3) 辐射畸变改正。利用由最小二乘影像匹配所求得辐射畸变改正参数h0, h1, 对上述重采样的结果作辐射改正。

(4) 计算左方影像窗口与经过几何、辐射改正后的右方影像窗口的灰度阵列之间的相关系数, 判断是否需要继续迭代。若相关系数小于前一次迭代后所求得的相关系数, 则可认为迭代结束, 也可以根据几何变形参数是否小于某个预定的阈值判断是否结束迭代。

(5) 采用最小二乘影像匹配, 利用误差方程解求变形参数的改正值dh0, dh1, da0, …。

(6) 计算变形参数。根据下列公式计算变形参数:

(7) 计算最佳匹配的点位。可用梯度的平方为权, 在左方影像窗口内对坐标作加权平均, 以它作为目标点坐标。

2 程序实现关键问题

最小二乘影像匹配算法程序有两个重要的部分, 一是利用误差方程求解变形参数改正数, 一个是几何、辐射变形的改正, 在这两个部分的程序实现过程中应注意两个关键问题:一是在每次迭代过程中均需进行几何变形改正, 即根据几何变形改正参数将左方影像窗口的影像坐标变换至右方影像阵列。其中计算出的坐标有可能超出影像范围, 导致程序错误, 应利用条件语句加以限制避免;二是变形参数改正数的解算需要利用上次迭代的改正结果, 结果应由原始影像坐标值、像素值与上次迭代计算出的改正参数求得, 避免误差累积。

3 实验分析与结论

影像匹配算法的评价标准有多种。数字摄影测量学中, 由于待匹配影像数据量较大, 匹配出的同名点坐标将用于后续的相对定向与绝对定向中, 因而对影像匹配速度与精度要求较高。以初匹配算法为相关系数法的最小二乘影像匹配算法为例, 下面将主要从匹配精度和速度两方面对最小二乘影像匹配性能进行分析, 并将其与相关系数法进行比较。

最小二乘影像匹配算法的匹配速度受到很多因素的制约, 包括初匹配算法的计算量、待匹配影像的数据量、匹配同名点对数等。其中, 可由使用者自行设置的计算窗口的大小对匹配速度影响很大。在选取相关系数法为初匹配算法的情况下, 对824*1265大小的影像进行单点匹配, 计算窗口为7*7、11*11、15*15时, 匹配时间分别为1.105s、1.928s、2.548s。由此可以说明, 计算窗口增大, 计算量随之增大, 匹配耗时也随之增长。因此应根据影像的特点选取合适的计算窗口, 避免由于计算窗口过大导致匹配速度过慢的情况。

最小二乘影像匹配算法充分利用了影像窗口内的信息进行平差计算, 考虑了辐射畸变与几何畸变两大影像灰度变形参数, 使影像匹配可以达到1/10甚至1/100像素的高精度, 即影像匹配精度可达到子像素等级。而相关系数法则采用计算窗口相对于待匹配影像不断移动一个像素并计算相关系数, 选取相关系数最大的窗口中心位置为匹配点的方法进行影像匹配[2], 匹配点坐标结果以整像素为单位, 无法达到子像素级, 故而当选取相关系数法为初匹配方法时, 最小二乘影像匹配算法的计算量明显大于单独的相关系数法影像匹配, 匹配速度相对较慢, 但匹配精度有所提高。

4 最小二乘影像匹配算法改进建议

实验结果分析可知, 最小二乘影像匹配算法虽然具有灵活、可靠、高精度等优点, 但仍存在一些待改善的问题。其中, 针对该算法计算量较大、匹配速度受限的情况, 在不降低匹配精度的基础上, 可从待匹配影像等方面入手, 减少待匹配位置, 进而减少计算量, 达到提高算法的匹配速率的目的。例如, 对于航空影像的影像匹配, 可将核线影像的制作与最小二乘影像匹配结合起来, 在条件充足的情况下, 先将普通的航空影像转换成核线影像, 再根据核线影像的特征进行影像匹配, 可在一定程度上提高匹配速度;也可将点特征提取算法与最小二乘影像匹配算法相结合, 在提取出两张影像特征点的基础上, 对左右影像上的特征点进行影像匹配, 忽略非点特征的待匹配位置, 进而提高匹配速度。实验证明, 使用核线影像进行最小二乘影像匹配和基于点特征提取的最小二乘影像匹配均可提高算法的匹配速度, 但仍存在待解决的问题——核线影像的制作需要已知影像的外方位元素等条件, 条件获取困难;两张影像上提取出的点特征可能不一致或出现错检漏检, 导致影像匹配错误, 因此还需根据具体的情况进行进一步的研究。

摘要：最小二乘影像匹配算法是数字摄影测量中常用的影像匹配方法之一。文章探讨了最小二乘影像匹配算法的基本原理, 介绍了该算法程序实现的关键问题, 分析了其实验结果并与相关系数法进行了比较, 最后为最小二乘影像匹配算法的改进提供了建议。

关键词：影像匹配,最小二乘,点特征

参考文献

[1]刘荣, 张祖勋, 张剑清.附有约束条件的最小二乘影像匹配及其在交通事故中的应用[J].华东地质学院学报, 2003 (3)

[2]张剑清, 潘励等.摄影测量学[M].武汉大学出版社, 2003

最小二乘映射 篇8

经典最小二乘平差通常考虑偶然误差只存在观测值中。然而,在数据的获取过程中,由于观测条件的限制,测量误差是不可避免的,所以在实际问题参数估计所建立的函数模型中观测向量和系数矩阵都应该包含误差。

近年来,总体最小二乘是测绘学科数据处理领域的热点研究问题。文献[1]针对总体最小二乘解算问题,应用测量平差中的间接平差原理推导了总体最小二乘的迭代逼近解算公式;文献[2~9]研究了基于总体最小二乘的大地测量反演理论及其应用问题,主要包括总体最小二乘解的性质[6,7]、算法[8,9]和在地壳应变参数[3,4]、火山形变[2]、断层位错模型参数[5]等反演问题中的应用;文献[10~13]研究了总体最小二乘方法在二维或三维坐标转换中的应用。

文献[14]系统阐述了整体最小二乘基本思想、原理,基于奇异值分解原理详细地阐述了总体最小二乘的解算方法,并且在对普通最小二乘和总体最小二乘深入比较的基础上,提出要注意总体最小二乘应用条件的问题;文献[15]研究了总体最小二乘法在工程测量上的应用。文献[16]和文献[17]综述了目前总体最小二乘的研究进展,剖析了总体最小二乘算法,其中主要有两大类,即基于奇异值分解(SVD)的方法和Schaffrin等发展的基于拉格朗日求极值的迭代解法。本文主要研究总体最小二乘平差问题的奇异值分解法和迭代解法,通过直线拟合的模拟算例对两类求解方法进行了比较和分析。

1 总体最小二乘平差问题求解方法

设线性函数模型为:

式(1)中:L为n×1观测值向量,A为n×m系数矩阵,X为m×1未知参数。

总体最小二乘的基本思想可以归纳为[1]:不仅观测向量L含有观测误差ΔL,系数矩阵A也存在误差ΔA,但是本文主要研究观测向量L含有观测误差和系数矩阵A误差服从同一分布,即不加权的总体最小二乘。

其线性方程为:

该方程可以改写为:

1.1 基于矩阵奇异值分解的总体最小二乘方法(SVD1)

将方程(2)表示为误差方程的形式[1]:

即:

式中其限制条件约束条件为:

或tr(E·ET)=tr(EAEAT+eeT)=vec(EA)Tvec(EA)+eeT=min,式中tr表示矩阵的迹。

对于误差模型的解算,数学上通常采用矩阵的奇异值分解来求解参数的TLS解。其解算过程为[1]:首先对增广矩阵[A L]进行奇异值分解:

式(7)中

其中,U=[U1U2]∈Rn×n为矩阵[A L]·[A L]T=AAT+LLT的n个特征值向量组成的正交阵,

V为矩阵[A L]·[A L]T=AAT+LLT的m+1个特征向量组成的正交矩阵。

为矩阵C的奇异值,其奇异值按顺序δ1≥δ2≥…≥δm+1排列。

未知参数X的估值为[1]:

根据上述解算原理,可以得到SVD解算的计算步骤[1]:

(1)建立函数模型,

(2)对增广矩阵进行奇异值分解,

(3)如果V22非奇异,则可得到

1.2 基于最小奇异值解法的总体最小二乘方法(SVD2)

对系数矩阵A进行SVD分解

对增广矩阵[A L]进行奇异值分解如上方法公式,由于奇异值向量vi∈Rm+1为矩阵的特征向量,所以X^满足如下特征向量方程[1]:

由式(13)第一行可以得到[1]:

1.3 迭代解法

迭代解法是基于拉格朗日极值的总体最小二乘解法,其平差准则或目标函数为:

式(15)中,vec(·)是将误差矩阵按列拉直得到的列向量,排列的顺序为从左到右。

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

式中,

,λ为维数n×1的Lagrange因子。

通过推导,得总体最小二乘的迭代逼近法如下[1]:

2 算例与分析

为了比较各种总体最小二乘方法计算得到的待定系数估值的差别,本文利用MATLAB软件给数据点位加入服从正态分布的偶然误差,然后分别用最小二乘方法(LS)、基于矩阵奇异值分解的总体最小二乘方法(SVD1)、基于最小奇异值解法的总体最小二乘方法(SVD2)以及基于拉格朗日极值的总体最小二乘迭代解法进行求解系数,对直线的不同情况进行了拟合,最后比较分析各种方法在不同情况下的差异。本文将两次迭代的Frobenius范数差设为0.001。

2.1 近似垂直的直线

模拟一条近似垂直于x轴的直线y=ax+b,其直线的设计模型为y=2000x-2000,自变量x从0.1到10之间每隔0.1模拟一个观测点,因变量y的范围为(-1800,18000),共100组观测数据,然后在每组观测点的x和y坐标分别添加服从正态分布N(0,0.01)的随机误差,以使得系数矩阵A和观测向量b都有误差。利用MATLAB软件计算了LS、SVD1、SVD2和迭代解法的参数估计结果,见表1。由于三种方法的估计值较为接近,三条拟合直线几乎重合,故拟合图省略。

2.2 近似水平的直线

模拟一条近似水平的直线y=ax+b,其直线的设计模型为y=x/500+1,自变量x从0.1到10之间每隔0.1模拟一个观测点,因变量y的范围为(1.0002,1.02),共100组观测数据,然后在每组观测点的x和y坐标分别添加服从正态分布N(0,0.01)的随机误差,以使得系数矩阵A和观测向量b都有误差。

表2为四种方法对近似水平直线模拟的估计值;同样三种方法的拟合图几乎重合,在此省略。

2.3 普通直线

模拟普通直线y=ax+b,直线的设计斜率的范围(10≤a≤10,且a≠0),b=1自变量x从0.1到10之间每隔0.1模拟一个观测点,共100组观测数据,然后在每组观测点的x和y坐标分别添加服从正态分布N(0,0.01)的随机误差,以使得系数矩阵A和观测向量b都有误差。三种方法的拟合图同样非常接近,在此省略。

2.4 算例分析

在表1中模拟近似垂直的直线拟合,设计的参数值a=2000,b=-2000;由表1中数据可知,LS解法参数估计的结果不能令人满意,这是由于直线斜率越大,x值波动误差对参数值的影响越大,因此无法得到精确的拟合结果。而用TLS方法拟合,无论直线的斜率多大,只要满足拟合条件,即所有点位到直线的垂直距离的平方和最小,故拟合的结果更为接近实际真值。由表1的数据可以看出,SVD1和SVD2拟合的结果几乎相同,两种方法实际上是等价的。本次迭代只进行了两次就结束了,奇异值解法与迭代解法的参数估计值也是很接近的,结果相差不大,这种情况下两种方法得到的估计值是一致的。

在表2中模拟近似水平的直线拟合,设计的参数值a=0.002,b=1;三种方法所得参数估计值也是很接近的,利用TLS方法的所得参数a和b的值要比LS方法更接近真值。分析其主要原因是在接近水平的直线中,y值的波动误差对参数的估计值影响较大,利用LS方法拟合的直线参数比起TLS方法与真值的偏差更大。

由于近似垂直的直线和近似水平的直线具有特殊性,故模拟普通直线的拟合,普通直线的斜率范围为(10≤a≤10),直线的截距均为1。其模拟结果显示,与真值相比,大多数直线TLS方法的参数估计值比LS方法更为接近真实值。奇异值分解的两种方法在普通直线的情况下所得到的结果是一样的,且经过迭代解算得到参数的估计值与奇异值分解法的结果几乎没有什么差异,模拟数据验证了两种算法的一致性。

3 结束语

奇异值分解(SVD)算法没有顾及EIV的随机模型,仅仅是基于数值逼近理论的总体最小二乘解,并非真正统计意义上的总体最小二乘解;因此,尽管SVD算法简单,但是显然无法获得大地测量领域普遍存在的观测数据不等精度以及相关情况下的平差模型的统计意义上的最佳估计值,这也是自SVD算法提出至今在大地测量领域的应用受到限制的主要原因[17]。

Schaffrin等发展的基于拉格朗日求极值的迭代解法是结合测量数据处理的特点,从测量平差的角度建立拉格朗日条件极值的目标函数通过推导构建起来的顾及系数矩阵含有误差情况下的一种扩展最小二乘方法;理论上,在等精度独立情况下,迭代解法与SVD算法是等价的。

本文在已有文献的基础上研究了总体最小二乘平差的奇异值分解法和迭代解法,运用MATLAB软件对直线的多种情况进行了数据模拟;通过对比分析,发现奇异值分解法与迭代解法的结果是一致的。迭代解法是从测量平差角度推导得到的,更适用于测量数据处理。总体最小二乘属于非线性估计,模型和算法的复杂性要远高于最小二乘估计[17]。本文只是讨论了等精度独立情况下常用的总体最小二乘解法,对于目前已有的加权总体最小二乘方法的比较和分析,并在此基础上提出算法简单、计算高效的加权总体最小二乘方法是今后研究的重点。

摘要：本文研究了总体最小二乘平差的两类方法,包括奇异值分解法和迭代解法,并且对两类方法 进行比较。运用MATLAB软件对直线的多种情况进行了数据模拟,比较分析了总体最小二乘平差的不同方法以及最小二乘方法的结果,发现奇异值分解法与迭代解法的结果是一致的;迭代解法是从测量平差角度推导得到的,更适用于测量数据处理。