精华内容
下载资源
问答
  • 下列四个数值最小
    千次阅读 多人点赞
    2018-08-12 01:12:52

    [深入浅出多旋翼飞控开发][姿态篇][四][非线性最小二乘与飞控传感器校准]

    作者:王伟韵
    QQ : 352707983
    Github

    前言

    搞好了传感器,那意味着飞控已经完成了一半。

    不用猜了,这句话正是鄙人说的:)。飞控的软硬件相关工作,绝大部分是围绕传感器展开的,它们的重要性,可能比你想象中的更大。

    硬件上,为了寻找性能优秀又符合产品实际需求的传感器如陀螺仪、加速度计、磁力计等,我们经常需要对来自不同厂商的样品进行实际性能测试与评估,即使实际上常用的型号就那么几个。为了提升飞控的性能,许多传感器问题我们还需要从硬件上进行解决,如降低电源噪声、隔离机身震动甚至还要为传感器提供恒定温度的工作环境。

    软件上,传感器校准、滤波与融合更是飞控程序中的重点,作为一名渣渣飞控工程师,浑浑噩噩工作数年,发现到头来大部分时间都是花在了传感器数据的处理上面,而且未来需要解决的问题,依然和传感器有关。

    本篇文章将结合实际飞控代码,以循序渐进的方式向读者讲述,如何使用最优化算法,解决飞控中最基本的传感器校准问题。

    一.传感器校准简述

    飞控上会使用各种各样的传感器,但对于多旋翼而言,最核心的是陀螺仪无疑。离开了陀螺仪,你甚至无法让一架多旋翼飞行器正常离地。其次则是加速度计,提供姿态观测,让飞行器可以得到一个较为客观的自身姿态观测信息,使得自稳飞行成为可能(能够自动保持自身姿态水平或一定角度)。除此之外,还有许许多多的的其它传感器,共同为飞控提供稳定飞行所需要的数据。

    其中,陀螺仪、加速度计和磁力计在飞行前(出厂前),均需要进行校准,才能达到正常工作所需要的性能。这便引出了一个问题:什么是传感器校准(标定)。

    答案很简单,所有事物都不是非黑即白,而是时刻处于一个不可预知的混沌状态。工厂按照特定的技术及工艺标准制造生产出一批传感器,这些传感器中所有个体的实际状态处于一个随机分布且均值为期望误差的序列上。简单来说,我们买来10个传感器,这些传感器均存在着特定误差且这10个传感器的误差各不一致。这就要求我们在实际使用这些传感器时,需要对每一个进行单独的校准,得到特定的校准参数,从而正常发挥传感器的性能。

    二.传感器误差类型

    接下来我们会面临下面两个问题:

    • 传感器有哪些误差
    • 如何测量传感器的误差

    如果从我们目前所使用的MEMS传感器的自身原理及制造工艺出发,分析传感器的误差来源,那可能需要长篇大论一番了,但那是超出作者能力的事情,也不是本篇文章所讨论的重点。

    根据实际经验,我们可以把传感器误差分为两大类:

    • 零偏误差(offset)
    • 刻度误差(scale)

    这两个是我们最容易理解,也是飞控中三大传感器均需要解决的误差类型,其中

    零偏误差: 指传感器在零激励下的输出误差。如陀螺仪在静止状态下,其角速度测量值理论为0,而实际上我们将飞控静止放置,陀螺仪依然会输出诸如1°/s之类的数据,这便是该陀螺仪传感器的零偏误差。

    刻度误差: 指传感器在一定激励输入下,其输出值与输入值的比值。比如我们将飞控放在一个以500°/s恒定角速度运动的转台上,而陀螺仪实际测量得到的角速度输出为495°/s,500 / 495 ≈ 1.01 便是该陀螺仪的刻度误差。

    从上述说明中我们可以看出,要正确测量传感器的误差,必须向传感器施加特定的激励信号,通过比较传感器的实际输出与激励输入量,计算得到传感器的误差。

    不同类型的传感器,所需要的激励信号源也不一样。如陀螺仪,要提供精确的激励信号,必须依靠专业的转台设备,这是个人及小公司均不具备的条件。好在对于大多数应用而言,陀螺仪的刻度误差在一个可以忍受的范围内,因此可被我们忽略。但对于精益求精的产品而言,要提升飞控性能到极致,刻度误差校准还是必须的,如某行业龙头,其飞控产品及整机产品,出厂前均需使用转台进行校准。缺少这些条件的个人爱好者,就只能暂时忽略这一项,只校准陀螺仪的零偏误差了,而大多数情况下,陀螺仪的零偏校准方式很简单,使用均值校准即可。因此,陀螺仪的校准不在本篇文章的讨论范围内。

    而本文中讨论的重点:加速度计磁力计,其信号激励源均来自于自然界,前者为地球的重力场,后者为地球的地磁场。

    三.传感器校准的简单方式

    1.加速度计

    假设有下面的情形,我们将飞控水平放置于水平面,测得加速度计z轴输出值 g 1 = 0.97 g g_1= 0.97g g1=0.97g,然后将飞控翻转并同样水平置于水平面,再测得加速度计z轴输出值 g 2 = − 0.99 g g_2 = -0.99g g2=0.99g

    设z轴零偏误差为 s c a l e z scale_z scalez,刻度误差为 o f f s e t z offset_z offsetz,已知实际重力加速度为g,可以得到下列两条等式:

    • ( g 1 − o f f s e t z ) ∗ s c a l e z = g (g_1 - offset_z) * scale_z = g (g1offsetz)scalez=g
    • ( g 2 − o f f s e t z ) ∗ s c a l e z = − g (g_2 - offset_z) * scale_z = -g (g2offsetz)scalez=g

    从而可以计算得到:

    • s c a l e z ≈ 1.02 scale_z ≈ 1.02 scalez1.02
    • o f f s e t z = − 0.01 g offset_z = -0.01g offsetz=0.01g

    将这个方法推广到其余两轴,便能得到其校准参数 s c a l e x , s c a l e y , o f f s e t x , o f f s e t y scale_x, scale_y,offset_x,offset_y scalex,scaley,offsetx,offsety

    缺陷:

    这个校准方法较为简单,但是其可靠性建立于一个理想条件下:加速度计在采集每一个方向的数据时,必须严格贴合水平面。否则只能得到一个不太准确的校准数值。实际应用中,我们很难保证校准时让飞控的每一个面保持水平,因为飞行器本身就是一个不规则物体,更多时候我们采集到的加速度数据为 [ 0.05 g , 0.08 g , 0 , 95 g ] [0.05g,0.08g,0,95g] [0.05g,0.08g,0,95g],而非 [ 0 , 0 , 0.98 g ] [0,0,0.98g] [0,0,0.98g]。或者说,这种方式的实现成本过高,难以实际应用。

    2.磁力计

    磁力计的简单校准方法原理等同于加速度计,即采集得到磁力计的三个轴的测量极值(最大与最小值),然后单独计算出每个轴的零偏误差与刻度误差。实现区别在于由于传感器激励信号源不一样,我们需要将磁力计在三维空间中进行各个方向的360°旋转,才能尽可能准确采集到极值。

    缺陷:

    实际应用中,磁力计校准对于飞行器使用者来说可能颇为频繁,而飞行器大小又各不一致,在很多情况下,我们是很难将飞行器翻来覆去地在各个方向上进行旋转,其可操作性非常低,且不方便普通用户的使用。如果旋转不全面,便无法采集到磁力计的真实极值,校准效果会大打折扣,甚至校准值严重偏离真实值。

    四.建立传感器误差模型

    对于程序员而言,解决实际问题的标准姿势如下:

    提出实际问题—>建立数学模型—>转换为程序—>求解

    而上一节中,我们描述了一个较为简单的用于加速度计与磁力计的校准方法,但是其前提条件较为苛刻,且校准效果不佳。于是,在这里我们提出了一个待解决的问题:找到一种使用便捷并且效果更好的校准算法。使用便捷,意味着约束条件越少越好,即在校准算法对传感器的采样点的数量和位置要求不能过于严格的情况下依然可以取得良好的校准效果。

    为了解决这个问题,得到传感器的误差参数,首先我们需要建立传感器的误差模型,并将其转换为数学方程,从而进行求解。

    以加速度计为例建立传感器误差模型

    设三轴加速度计的零偏误差与刻度误差分别为
    o x , o y , o z , s x , s y , s z o_x ,o_y,o_z,s_x,s_y,s_z ox,oy,oz,sx,sy,sz
    某一时刻该三轴加速度计的测量值分别为
    x , y , z x,y,z x,y,z
    设当前时刻真实加速度向量为
    a = ( a x , a y , a z ) a=(a_x,a_y,a_z) a=(ax,ay,az)
    有以下关系:
    { a x = ( x − o x ) s x a y = ( y − o y ) s y a z = ( z − o z ) s z \begin{cases} a_x=(x-o_x)s_x\\ a_y=(y-o_y)s_y\\ a_z=(z-o_z)s_z\\ \end{cases} ax=(xox)sxay=(yoy)syaz=(zoz)sz

    静止状态下,不管加速度计方向如何,其测量得到的加速度为重力加速度在各轴上的投影。而在地球上,重力加速度向量模值总是等于 1 g 1g 1g。假设 a x , a y , a z a_x,a_y,a_z ax,ay,az单位为 g g g,便有以下等式成立:
    a x 2 + a y 2 + a z 2 = 1 a_x^2+a_y^2+a_z^2=1 ax2+ay2+az2=1

    ( x − o x ) 2 s x 2 + ( y − o y ) 2 s y 2 + ( z − o z ) 2 s z 2 = 1 (x-o_x)^2s_x^2+(y-o_y)^2s_y^2+(z-o_z)^2s_z^2=1 (xox)2sx2+(yoy)2sy2+(zoz)2sz2=1

    至此,我们将加速度计的误差求解转化成为了一个数学问题: 已知 x , y , z x,y,z x,y,z ,由上一条等式,求解出 o x , o y , o z , s x , s y , s z o_x ,o_y,o_z,s_x,s_y,s_z ox,oy,oz,sx,sy,sz

    五.求解误差方程

    上一节中,我们得到了传感器的误差方程: ( x − o x ) 2 s x 2 + ( y − o y ) 2 s y 2 + ( z − o z ) 2 s z 2 = 1 (x-o_x)^2s_x^2+(y-o_y)^2s_y^2+(z-o_z)^2s_z^2=1 (xox)2sx2+(yoy)2sy2+(zoz)2sz2=1
    求解方程,对于经历过九年义务教育的我们来说是一件非常熟悉的事情,但这里存在着几个棘手的问题:

    • 该方程为六元二次方程,存在着6个未知数,如果只存在 x , y , z x,y,z x,y,z一组观测值,那将会有无穷多解。因此我们需要获得多组观测值 x i , y i , z i x_i,y_i,z_i xi,yi,zi,以构建方程组,从而求出唯一解,满足方程组里的每一个方程

    • 对于六元二次方程而言,使用六组独立的观测值,构建出包含6个方程的方程组,理论上可以求出方程组的唯一解。可是事实上,我们从传感器上获取到的观测值,都是带有一定随机噪声的,即观测值并非完全准确,从而导致求得的方程解未必精确

    • 使用计算机,我们可以轻松求出线性方程组的唯一精确解(如高斯消元法),而非线性方程问题,无论是从理论上还是计算方法上,都要复杂得多,大部分时候,对于无特定形式的非线性方程 f ( x ) = 0 f(x)=0 f(x)=0,是没有办法直接求得精确解的。

    为了解决上述问题,这里引进了一种数学方法:最小二乘法

    最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。

    对于现实中的绝大部分问题,我们都是无法求得精确解的,于是退而求其次,使用最小二乘法,我们可以得到其近似解。

    以传感器误差方程为例,设 r r r为某特定解下方程的误差,有
    r = 1 − ( x − o x ) 2 s x 2 + ( y − o y ) 2 s y 2 + ( z − o z ) 2 s z 2 r=1 - (x-o_x)^2s_x^2+(y-o_y)^2s_y^2+(z-o_z)^2s_z^2 r=1(xox)2sx2+(yoy)2sy2+(zoz)2sz2
    对于方程组而言,其误差的平方和函数为
    S = ∑ i = 0 n r i 2 S=\sum_{i=0}^{n} r_i^2 S=i=0nri2
    所谓最小二乘法,便是求得一组特定解,使得 S S S最小(求函数极值),该条件下的解,便是方程组的最优近似解。

    对于线性方程组,可将误差函数 S = 0 S=0 S=0转换为矩阵方程并进行求解即可,可惜现实问题大部分均为非线性问题,如本文中的传感器误差方程便是一个非线性方程,是无法变换为矩阵方程形式的,因此我们必须换一个思路来解决问题。

    假定该函数存在着局部最小值,那给定一个初始解,通过不断地迭代计算,从而找到误差平方和函数 S S S的局部最优解,这便是非线性最小二乘的基本思想。

    接下来,我们将讨论求解非线性最小二乘的具体算法。

    六.高斯牛顿法

    非线性最小二乘的求解算法有许多种,比较常见的有高斯牛顿法Levenberg-Marquardt法(LM),DogLeg法等。其中高斯牛顿法是最为经典的一种算法,APM飞控中便是使用了该算法,而LM法可以看成是高斯牛顿法的改良版,鲁棒性更好,在实际工程中应用更广泛,天穹飞控中则使用了LM法。而LM法的基础部分和高斯牛顿法完全一致,因此我们从高斯牛顿法开始讲起。

    由于算法原理细节部分较为复杂,需要用到大量公式,为了方便理解,将结合天穹飞控的实际代码进行讲解,高斯牛顿法和LM法的具体实现大部分是一样的,所以可以结合着来看。对于初学者来说彻底弄懂这一部分的原理对自己能力的提升帮助较大,因为此类最优化算法在飞控之外的非常多工程领域都是有着大量应用的,如计算机视觉和机器学习等。

    天穹飞控中的LM法完整代码

    要搞明白高斯牛顿法,我们要先来了解一下牛顿法

    1.牛顿法

    为求得函数最优解(极值点),牛顿法将问题转换为了求解 f ′ ( x ) = 0 f^{'} (x)=0 f(x)=0

    首先给定一个初始解 x 0 x_0 x0,然后在 x 0 x_0 x0处进行一阶泰勒展开,得到
    f ′ ( x ) ≈ f ′ ( x 0 ) + ( x − x 0 ) f ′ ′ ( x 0 ) f^{'} (x) \approx f^{'} (x_0) + (x-x_0)f^{''} (x_0) f(x)f(x0)+(xx0)f(x0)

    f ′ ( x 0 ) + ( x − x 0 ) f ′ ′ ( x 0 ) = 0 f^{'} (x_0) + (x-x_0)f^{''} (x_0) = 0 f(x0)+(xx0)f(x0)=0
    求解得到 x x x,相比于 x 0 x_0 x0,有 f ′ ( x ) < f ′ ( x 0 ) f^{'} (x) < f^{'} (x_0) f(x)<f(x0),即新解 x x x比初始解更接近最优解。当然该解和最优解之间还有一定未知的距离,于是我们将新解作为初始解,重复这一步骤,经过不断地迭代计算,最终会收敛于 f ′ ( x ) = 0 f^{'} (x)=0 f(x)=0

    牛顿法寻找极值点的动画示意图:
    牛顿法寻找极值点的动画示意图
    牛顿法的优点:

    • 收敛速度快,且可用于求解任意连续函数的最优化问题

    牛顿法的缺点:

    • 可以看出,牛顿法实际上等同于寻找函数的局部极值点,因此必须选取一个相对接近的初始点,否则会导致无法正确收敛

    • 需要求解目标函数的Hessian矩阵(二阶导)的逆矩阵,计算比较复杂

    2.高斯牛顿法

    高斯牛顿法是一种非线性最小二乘最优化方法。 其利用了目标函数的泰勒展开式把非线性函数的最小二乘化问题化为每次迭代的线性函数的最小二乘化问题

    高斯牛顿法实际上是牛顿法的简化形式,或者说改进版。多维情况下,当维数较大时,牛顿法中的Hessian矩阵求逆是往往行不通的,而高斯牛顿法对牛顿法中的Hessian矩阵进行了化简,使得问题可解。

    下面讲解高斯牛顿法的具体实现

    首先定义目标函数
    S ( β ) = ∑ i = 0 5 r i 2 ( β ) S(\beta)= \sum_{i=0}^{5}r _i^2( \beta) S(β)=i=05ri2(β)
    其中
    β = [ β 0 , β 1 , β 2 , β 3 , β 4 , β 5 ] \beta=[\beta_0,\beta_1,\beta_2,\beta_3,\beta_4,\beta_5] β=[β0,β1,β2,β3,β4,β5]
    对应了加速度计误差方程中的六个变量,即
    [ o x , o y , o z , s x , s y , s z ] [o_x ,o_y,o_z,s_x,s_y,s_z] [ox,oy,oz,sx,sy,sz]
    r i 2 ( β ) r _i^2( \beta) ri2(β)为传感器误差方程的误差平方和函数,即
    r i 2 ( β ) = ( 1 − ( x i − o x ) 2 s x 2 + ( y i − o y ) 2 s y 2 + ( z i − o z ) 2 s z 2 ) 2 r _i^2( \beta)=(1 - (x_i-o_x)^2s_x^2+(y_i-o_y)^2s_y^2+(z_i-o_z)^2s_z^2)^2 ri2(β)=(1(xiox)2sx2+(yioy)2sy2+(zioz)2sz2)2

    设定一个初解 β α \beta_\alpha βα,将目标函数 S ( β ) S(\beta) S(β) β α \beta_\alpha βα附近泰勒展开,并舍弃高阶无穷项,近似有:
    S ( β ) = S ( β α ) + [ ▽ S ( β α ) ] T ( β − β α ) + 1 2 ( β − β α ) T H S ( β α ) ( β − β α ) S(\beta)= S(\beta_\alpha)+[\bigtriangledown S(\beta_\alpha)]^T(\beta-\beta_\alpha)+\frac{1}{2}(\beta-\beta_\alpha)^TH_S(\beta_\alpha)(\beta-\beta_\alpha) S(β)=S(βα)+[S(βα)]T(ββα)+21(ββα)THS(βα)(ββα)
    我们的目标是求出特定解,使 S ( β ) S(\beta) S(β)有最小值。当函数有最小值(属于极小值)时,其一阶导为0,即
    ▽ S ( β ) = ▽ S ( β α ) + H S ( β α ) ( β − β α ) = 0 \bigtriangledown S(\beta)=\bigtriangledown S(\beta_\alpha)+H_S(\beta_\alpha)(\beta-\beta_\alpha)=0 S(β)=S(βα)+HS(βα)(ββα)=0
    从而得到了高斯牛顿法的迭代公式
    β = β α − [ H S ( β α ) ] − 1 ▽ S ( β α ) \beta=\beta_\alpha - [H_S(\beta_\alpha)]^{-1}\bigtriangledown S(\beta_\alpha) β=βα[HS(βα)]1S(βα)

    △ = [ H S ( β α ) ] − 1 ▽ S ( β α ) \bigtriangleup = [H_S(\beta_\alpha)]^{-1}\bigtriangledown S(\beta_\alpha) =[HS(βα)]1S(βα)
    则有
    β = β α − △ \beta=\beta_\alpha - \bigtriangleup β=βα
    这里便是迭代的关键了,给定初解后,每一次计算出 △ \bigtriangleup ,便得到了一个新解,然后使用新解不断进行迭代计算,使其逐渐逼近最优解。

    可能看到这里,你依然是一脸懵逼,但没有关系,你只需要知道,为了求解出函数的最优解,也就是传感器的最优校准参数,我们只需要不断去计算 △ \bigtriangleup 就行了。

    接下来便讲解如何计算 △ \bigtriangleup (这里是整个算法实现的重点)

    △ \bigtriangleup 中的 H S H_S HS S S S函数的Hessian矩阵(实际上是S函数的二阶导)
    H S = [ ∂ 2 S ∂ β 0 2 ⋯ ∂ 2 S ∂ β 0 β 5 ⋮ ⋱ ⋮ ∂ 2 S ∂ β 5 β 0 ⋯ ∂ 2 S ∂ β 5 2 ] H_S=\left[ \begin{matrix} \frac{\partial^2 S}{\partial \beta_0^2}&\cdots&\frac{\partial^2 S}{\partial \beta_0 \beta_5}&\\ \vdots&\ddots&\vdots&\\ \frac{\partial^2 S}{\partial \beta_5 \beta_0}&\cdots&\frac{\partial^2 S}{\partial \beta_5^2}&\\ \end{matrix} \right] HS=β022Sβ5β02Sβ0β52Sβ522S
    前文提到,高斯牛顿法对牛顿法所做的改进,便是对Hessian矩阵进行了简化,其中有
    ∂ 2 S ∂ β i β j = 2 ∑ k ( ∂ r k ∂ β i ∂ r k ∂ β j + r k ∂ 2 r k ∂ β i β j ) \frac{\partial^2 S}{\partial \beta_i \beta_j}=2\sum_{k}(\frac{\partial r_k}{\partial \beta_i} \frac{\partial r_k}{\partial \beta_j} + r_k\frac{\partial^2 r_k}{\partial \beta_i \beta_j} ) βiβj2S=2k(βirkβjrk+rkβiβj2rk)
    忽略二阶微分项,近似有
    ∂ 2 S ∂ β i β j = 2 ∑ k ∂ r k ∂ β i ∂ r k ∂ β j \frac{\partial^2 S}{\partial \beta_i \beta_j}=2\sum_{k}\frac{\partial r_k}{\partial \beta_i} \frac{\partial r_k}{\partial \beta_j} βiβj2S=2kβirkβjrk
    于是有
    H S = 2 [ ∂ r 0 ∂ β 0 ⋯ ∂ r 0 ∂ β 5 ⋮ ⋱ ⋮ ∂ r 5 ∂ β 0 ⋯ ∂ r 5 ∂ β 5 ] T [ ∂ r 0 ∂ β 0 ⋯ ∂ r 0 ∂ β 5 ⋮ ⋱ ⋮ ∂ r 5 ∂ β 0 ⋯ ∂ r 5 ∂ β 5 ] = 2 J r T J r H_S=2 \left[ \begin{matrix} \frac{\partial r_0}{\partial \beta_0}&\cdots&\frac{\partial r_0}{\partial \beta_5}&\\ \vdots&\ddots&\vdots&\\ \frac{\partial r_5}{\partial \beta_0}&\cdots&\frac{\partial r_5}{\partial \beta_5}&\\ \end{matrix} \right]^T \left[ \begin{matrix} \frac{\partial r_0}{\partial \beta_0}&\cdots&\frac{\partial r_0}{\partial \beta_5}&\\ \vdots&\ddots&\vdots&\\ \frac{\partial r_5}{\partial \beta_0}&\cdots&\frac{\partial r_5}{\partial \beta_5}&\\ \end{matrix} \right] =2J_r^TJ_r HS=2β0r0β0r5β5r0β5r5Tβ0r0β0r5β5r0β5r5=2JrTJr
    这里的 J r J_r Jr实际上就是函数 r r r(传感器误差方程的误差函数)的雅可比矩阵

    △ \bigtriangleup 中的另外一个成员 ▽ S \bigtriangledown S S S S S函数的一阶导的列向量,即
    ▽ S = [ ∂ S ∂ β 0 ⋮ ∂ S ∂ β 5 ] \bigtriangledown S = \left[ \begin{matrix} \frac{\partial S}{\partial \beta_0}\\ \vdots\\ \frac{\partial S}{\partial \beta_5}\\ \end{matrix} \right] S=β0Sβ5S
    又有
    ∂ S ∂ β i = ∑ k ∂ ( r k 2 ) ∂ β i = 2 ∑ k ∂ r k ∂ β i r k \frac{\partial S}{\partial \beta_i}=\sum_{k}\frac{\partial (r_k^2)}{\partial \beta_i} =2 \sum_{k}\frac{\partial r_k}{\partial \beta_i} r_k βiS=kβi(rk2)=2kβirkrk
    故$\bigtriangledown S $可改写为
    ▽ S = [ ∂ S ∂ β 0 ⋮ ∂ S ∂ β 5 ] = 2 [ ∂ r 0 ∂ β 0 ⋯ ∂ r 0 ∂ β 5 ⋮ ⋱ ⋮ ∂ r 5 ∂ β 0 ⋯ ∂ r 5 ∂ β 5 ] T [ r 0 ⋮ r 5 ] = 2 J r T r \bigtriangledown S =\left[ \begin{matrix} \frac{\partial S}{\partial \beta_0}\\ \vdots\\ \frac{\partial S}{\partial \beta_5}\\ \end{matrix} \right] =2 \left[ \begin{matrix} \frac{\partial r_0}{\partial \beta_0}&\cdots&\frac{\partial r_0}{\partial \beta_5}&\\ \vdots&\ddots&\vdots&\\ \frac{\partial r_5}{\partial \beta_0}&\cdots&\frac{\partial r_5}{\partial \beta_5}&\\ \end{matrix} \right]^T \left[ \begin{matrix} r_0\\ \vdots\\ r_5\\ \end{matrix} \right] =2J_r^Tr S=β0Sβ5S=2β0r0β0r5β5r0β5r5Tr0r5=2JrTr

    于是我们终于得到了 △ \bigtriangleup 的最终形式
    △ = [ J r ( β α ) T J r ( β α ) ] − 1 J r ( β α ) T r ( β α ) \bigtriangleup = [J_r(\beta_\alpha)^T J_r(\beta_\alpha)]^{-1} J_r(\beta_\alpha)^T r(\beta_\alpha) =[Jr(βα)TJr(βα)]1Jr(βα)Tr(βα)
    或者换一个形式
    J r ( β α ) T J r ( β α ) △ = J r ( β α ) T r ( β α ) J_r(\beta_\alpha)^T J_r(\beta_\alpha) \bigtriangleup =J_r(\beta_\alpha)^T r(\beta_\alpha) Jr(βα)TJr(βα)=Jr(βα)Tr(βα)
    这就成了一个线性方程组,未知量是 △ \bigtriangleup ,使用高斯消元法,便能求解出 △ \bigtriangleup 的唯一解。

    至此,高斯牛顿法的所有计算步骤均结束,我们只需要不断进行迭代计算:求解出 △ \bigtriangleup ,从而得到新解 β \beta β,然后进行下一次迭代,不断逼近函数最优解。

    随着迭代次数的增加,最终会得到一个无限接近最优的解,但是考虑到时间成本和实际使用需求,在校准程序中我们设置一个阈值 e p s eps eps,当方程解达到一定精度时(定义为迭代步长的平方 △ 2 < e p s \bigtriangleup ^2 < eps 2<eps),停止迭代,结束此次校准。

    看完了,还是有点懵,怎么办?

    下面贴上高斯牛顿法的完整实现代码,结合代码进行理解和分析,事半功倍,如果还没看懂,那就多看几遍。

    主程序

    /**********************************************************************************************************
    *函 数 名: GaussNewton
    *功能说明: 高斯牛顿法求解传感器误差方程,得到校准参数
    *形    参: 传感器采样数据(6组) 零偏误差指针 比例误差指针 数据向量长度
    *返 回 值: 无
    **********************************************************************************************************/
    void GaussNewton(Vector3f_t inputData[6], Vector3f_t* offset, Vector3f_t* scale, float length)
    {
        uint32_t cnt    = 0;
        double   eps    = 0.000000001;
        double   change = 100.0;
        float    data[3];  
        float    beta[6];      //方程解
        float    delta[6];     //迭代步长
        float    JtR[6];       //梯度矩阵
        float    JtJ[6][6];    //Hessian矩阵
       
        //设定方程解初值
        beta[0] = beta[1] = beta[2] = 0;
        beta[3] = beta[4] = beta[5] = 1 / length;
    
        //开始迭代,当迭代步长小于eps时结束计算,得到方程近似最优解
        while(change > eps) 
        {
            //矩阵初始化
            ResetMatrices(JtR, JtJ);
    
            //计算误差方程函数的梯度JtR和Hessian矩阵JtJ
            for(uint8_t i=0; i<6; i++) 
            {
                data[0] = inputData[i].x;
                data[1] = inputData[i].y;
                data[2] = inputData[i].z;
                UpdateMatrices(JtR, JtJ, beta, data);
            }
    
            //高斯消元法求解方程:JtJ * delta = JtR,得到delta
            GaussEliminateSolveDelta(JtR, JtJ, delta);
    
            //计算迭代步长
            change = delta[0]*delta[0] +
                     delta[0]*delta[0] +
                     delta[1]*delta[1] +
                     delta[2]*delta[2] +
                     delta[3]*delta[3] / (beta[3]*beta[3]) +
                     delta[4]*delta[4] / (beta[4]*beta[4]) +
                     delta[5]*delta[5] / (beta[5]*beta[5]);
    
            //更新方程解
            for(uint8_t i=0; i<6; i++) 
            {
                beta[i] -= delta[i];
            }
                
            //限制迭代次数
            if(cnt++ > 100)
                break;
        }
    
        //更新校准参数
        scale->x  = beta[3] * length;
        scale->y  = beta[4] * length;
        scale->z  = beta[5] * length;
        offset->x = beta[0];
        offset->y = beta[1];
        offset->z = beta[2];
    }
    

    主程序中调用了三个子函数,其中:

    • ResetMatrices,用于初始化矩阵变量
    static void ResetMatrices(float JtR[6], float JtJ[6][6])
    {
        int16_t j,k;
        for(j=0; j<6; j++) 
        {
            JtR[j] = 0.0f;
            for(k=0; k<6; k++) 
            {
                JtJ[j][k] = 0.0f;
            }
        }
    }
    
    • UpdateMatrices,用于计算求解 △ \bigtriangleup 所用到的 J r ( β α ) T J r ( β α ) J_r(\beta_\alpha)^T J_r(\beta_\alpha) Jr(βα)TJr(βα)(Hessian矩阵)和 J r ( β α ) T r ( β α ) J_r(\beta_\alpha)^T r(\beta_\alpha) Jr(βα)Tr(βα)这两个矩阵
    static void UpdateMatrices(float JtR[6], float JtJ[6][6], float beta[6], float data[3])
    {
        int16_t j, k;
        float dx, b;
        float residual = 1.0;
        float jacobian[6];
        
        for(j=0; j<3; j++) 
        { 
            b = beta[3+j];
            dx = (float)data[j] - beta[j];
            //计算残差 (传感器误差方程的误差)
            residual -= b*b*dx*dx;
            
            //计算雅可比矩阵
            jacobian[j] = 2.0f*b*b*dx;
            jacobian[3+j] = -2.0f*b*dx*dx;
        }
        
        for(j=0; j<6; j++) 
        {
            //计算函数梯度
            JtR[j] += jacobian[j]*residual;
            
            for(k=0; k<6; k++) 
            {
                //计算Hessian矩阵(简化形式,省略二阶偏导),即雅可比矩阵与其转置的乘积
                JtJ[j][k] += jacobian[j]*jacobian[k];
            }
        }
    }
    
    • GaussEliminateSolveDelta,使用高斯消元法求解 △ \bigtriangleup
    static void GaussEliminateSolveDelta(float JtR[6], float JtJ[6][6], float delta[6])
    {
        int16_t i,j,k;
        float mu;
       
        //逐次消元,将线性方程组转换为上三角方程组
        for(i=0; i<6; i++)
        {
            //若JtJ[i][i]不为0,将该列在JtJ[i][i]以下的元素消为0
            for(j=i+1; j<6; j++)
            {
                mu = JtJ[i][j] / JtJ[i][i];
                if(mu != 0.0f) 
                {
                    JtR[j] -= mu * JtR[i];
                    for(k=j; k<6; k++)
                    {
                        JtJ[k][j] -= mu * JtJ[k][i];
                    }
                }
            }
        }
    
        //回代得到方程组的解
        for(i=5; i>=0; i--)
        {
            JtR[i] /= JtJ[i][i];
            JtJ[i][i] = 1.0f;
            
            for(j=0; j<i; j++)
            {
                mu = JtJ[i][j];
                JtR[j] -= mu * JtR[i];
                JtJ[i][j] = 0.0f;
            }
        }
    
        for(i=0; i<6; i++)
        {
            delta[i] = JtR[i];
        }
    }
    

    事实上,在实际应用中使用高斯牛顿法,会出现许许多多的问题,因此下面介绍一种基于高斯牛顿改良优化而来的算法。

    七.Levenberg-Marquardt法

    实际应用中,高斯牛顿法有着如下缺点:

    • 初始点选取较为严苛,若初始点距离极小值点过远,迭代步长过大会导致迭代下一代的函数值不一定小于上一代的函数值
    • 高斯牛顿法中的Hessian矩阵在某些情况下可能是非正定的(不可逆)
    • 残差 r r r过大时可能导致高斯牛顿法无法收敛

    而Levenberg-Marquardt法(LM法,又称阻尼最小二乘法)改善了高斯牛顿法的不足之处,并结合了高斯牛顿法和梯度下降法的优点,而实际上所做的改动非常小。

    在高斯牛顿法中,迭代步长 △ \bigtriangleup 被定义为
    △ = [ H S ( β α ) ] − 1 ▽ S ( β α ) \bigtriangleup = [H_S(\beta_\alpha)]^{-1}\bigtriangledown S(\beta_\alpha) =[HS(βα)]1S(βα)
    而在LM法中,加入了阻尼因子 μ \mu μ,有
    △ = [ H S ( β α ) + μ I ] − 1 ▽ S ( β α ) \bigtriangleup = [H_S(\beta_\alpha) + \mu I]^{-1}\bigtriangledown S(\beta_\alpha) =[HS(βα)+μI]1S(βα)
    从而消除了Hessian矩阵非正定的情况,更重要的是迭代过程中阻尼因子 μ \mu μ是可变的,用于调整下降方向,当:

    • μ \mu μ趋近于0时,算法退化为高斯牛顿法
    • μ \mu μ很大时,效果相当于梯度下降法

    使用LM法时, 需要加入一些 μ \mu μ的调整策略,使得距离解较远的时候,使用梯度下降法,距离极小值较近的时候,使 μ \mu μ减小,这样就近似高斯牛顿法,能得到比梯度下降法更快的收敛速度。

    算法的具体实现这里就不贴出来了,可以到这里查看:天穹飞控中的LM法

    八.传感器数据采集

    前几节介绍了如何使用非线性最小二乘法校准传感器,下面简单说明校准时采集所需传感器数据的具体步骤。

    1.加速度计

    在天穹飞控中,加速度校准使用标准的六面采集法,需要将飞控按照上、下、前、后、左、右六个方向静止放置一段时间(2s左右),顺序随意且飞控方向大致准确即可,飞控自动检测当前加速度计方向,并采集多组数据取平均值,最后得到6组方向各异的数据,导入LM校准函数,便能求出加速度计的6个校准参数。

    加速度校准采集代码

    天穹飞控兼容mavlink,可以使用QGroundControl地面站完成加速度计的校准步骤,如下图:
    使用QGC校准加速度计

    点击左侧的Accelerometer按钮,确认OK,飞控便会进入加速度校准流程,同时QGC会显示出当前校准进度以及图片提示,当对应方向采集数据完毕时,对应图片外框会变绿,红色的则是还没有采集数据的方向。

    2.磁力计

    磁力计的数据采集方式采用比较通用的两圈旋转方式,即:将飞控水平旋转一圈,然后再竖直旋转一圈。在这个过程中,记录下各轴的最大和最小值数据。

    磁力计校准采集代码

    同样可以使用QGC进行磁力计校准,点击Compass并确认OK。开始按照上述步骤旋转飞控,完成校准。

    然而相对加速度计校准而言,磁力计校准对算法的要求更高。由于加速度计校准环境通常比较单一,且静止状态下加速度计数据基本不受外界干扰,所以采集到的数据都是比较准确的,使得算法很容易收敛,实测大概只需4-5次的迭代运算,就可以达到所需精度。

    而磁力计数据采集在实际应用过程中,可能存在以下问题:

    • 校准环境比较复杂,而且校准频率非常频繁。用户可能是在室内校准,也可能是在城市的建筑间校准,或者是在充满了未知磁场干扰的环境下校准。这便导致了在校准期间,可能由于各种未知的外部干扰,使得磁力计数据采集出现“野值”

    • 为了简化校准步骤,我们采用了较为简单的两圈旋转方式,但由于地磁场的分布方向原因,实际上这样我们采集到的磁力计数据并不是在整个球面上分布均匀的,如下图所示(图片来自正在开发中的天穹地面站):
      磁力计数据分布示意图
      为了解决问题1,我们对采集到的数据做了一定的滤波平滑处理,同时实时计算数据的一个平均模值,当新值的大小超过平均模值一定比例时,抛弃这个值,避免在构建误差方程组时,引入包含较大干扰量的观测值,导致算法无法收敛或最终计算出来的传感器校准参数存在误差。

    问题2考验了算法的性能,即在样本分布不均匀的情况下是否依然可以收敛,经过测试和验证,天穹飞控中使用的LM法是能够经受得住这个考验的。

    九.总结

    本篇文章介绍了飞控中最基本的传感器误差问题,并提出了有效的解决方法。其中所使用到的非线性最小二乘算法,在现代计算机工程领域有着非常广泛的应用,所以即使对飞控技术兴趣不大的童鞋,本文也是值得一看的。

    而对于我们的飞控而言,完成了传感器基本误差的校准,只是相当于迈出了第一步,后面还有无数难题待解决,再接下来的教程中,会陆续针对不同的问题寻求最优的解决方案,当然由于本人水平有限,只能起到抛砖引玉的作用,希望有更多的有识之士,加入这个开源飞控项目,一起完善,共同进步。

    最后再留一下传送门:

    天穹飞控
    飞控交流论坛
    飞控交流Q群:472648354

    更多相关内容
  • 数值分析:数据的最小二乘拟合

    千次阅读 2021-10-10 15:34:34
    在已知某天在不同时间的前温度高低,借用最小二乘法确定这一天的气温变化规律。通过MATLAB编程,选取适当函数进行求解绘图。 2实验内容 假定某天的气温变化记录如下表所示: 时间(t) 0 ...

     1 实验目的

    在已知某天在不同时间的前温度高低,借用最小二乘法确定这一天的气温变化规律。通过MATLAB编程,选取适当函数进行求解绘图。

    2 实验内容

    假定某天的气温变化记录如下表所示:

    时间(t)

    0

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    温度(x(t))

    15

    14

    14

    14

    14

    15

    16

    18

    20

    22

    23

    25

    28

    时间(t)

    13

    14

    15

    16

    17

    18

    19

    20

    21

    22

    23

    24

    温度(x(t))

    31

    32

    31

    29

    27

    25

    24

    22

    20

    18

    17

    16

    试用最小二乘法确定这一天的气温变化规律,考虑用下列类型的函数,计算误差平方和,并作图比较效果。

    (1)二次函数

    (2)三次函数

    (3)四次函数

    (4)函数式中a,b,c为常数(拟合参数)

    3 实验知识点

    最小二乘法、非线性拟合

    4 算法思想

    所谓拟合是指已知某函数的若干离散函数值{f1,f2,…,fn},通过调整该函数中若干待定系数f(λ1, λ2,…,λn),使得该函数与已知点集的差别最小。拟合的方法除了最小二乘法外,还有拉格朗日插值法、牛顿插值法、牛顿迭代法、区间二分法、弦截法、雅克比迭代法和牛顿科特斯数值积分发等方法。

    最小二乘拟合步骤:

    2为最小,按ni=1这样的标准定义的bai拟合函数称为du最小二乘拟合,是离散情形下的最佳平方逼近.对给定数据点{(Xi,Yi)}(i=0,1,…,m),在取定的函数类Φ 中,求p(x)∈Φ ,使误差的平方和E2最小,E2=∑[p(Xi)-Yi]2。从几何意义上讲,就是寻求与给定点 {(Xi,Yi)}(i=0,1,…,m)的距离平方和为最小的曲线y=p(x)。函数p(x)称为拟合函数或最小二乘解,求拟合函数p(x)的方法称为曲线拟合的最小二乘法。

    5 实验代码

    5.1 二次函数

    t=0:24;

    y=[15 14 14 14 14 15 16 18 20 22 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16];

    [p2,e2]=polyfit(t,y,2);

    p=polyfit(t,log(y),2)

    b=-p(1)

    c=p(2)/(2*b)

    a=exp(p(3)+b*c^2)

    x1=0:0.1:24;

    y2=polyval(p2,x1);

    plot(t,y,'*',x1,y2)

    xlabel('时间(t)');

    ylabel('气温(℃)');

    legend('气温散点图','拟合图');

    title('二次函数');

    e2

    5.2 三次函数

    t=0:24;

    y=[15 14 14 14 14 15 16 18 20 22 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16];

    [p3,e3]=polyfit(t,y,3);

    p=polyfit(t,log(y),3)

    b=-p(1)

    c=p(2)/(2*b)

    a=exp(p(3)+b*c^2)

    x1=0:0.1:24;

    y3=polyval(p3,x1);

    plot(t,y,'r--',x1,y3)

    xlabel('时间(t)');

    ylabel('气温(℃)');

    legend('气温散点图','拟合图');

    title('三次函数');

    e2

    5.3 四次函数

    t=0:24;

    y=[15 14 14 14 14 15 16 18 20 22 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16];

    [p4,e4]=polyfit(t,y,4);

    p=polyfit(t,log(y),4)

    b=-p(1)

    c=p(2)/(2*b)

    a=exp(p(4)+b*c^2)

    x1=0:0.1:24;

    y4=polyval(p3,x1);

    plot(t,y,'r--',x1,y3)

    xlabel('时间(t)');

    ylabel('气温(℃)');

    legend('气温散点图','拟合图');

    title('四次函数');

    e2

    6.4 函数,式中a,b,c为常数(拟合参数)

    t=0:24;

    y=[15 14 14 14 14 15 16 18 20 22 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16];

    [p2,e2]=polyfit(t,y,2);[p3,e3]=polyfit(t,y,3);[p4,e4]=polyfit(t,y,4);

    p=polyfit(t,log(y),2)

    b=-p(1)

    c=p(2)/(2*b)

    a=exp(p(3)+b*c^2)

    x1=0:0.1:24;

    ye=a*exp(-b*(x1-c).^2);y5=a*exp(-b*(t-c).^2);e5=sqrt(sum((y-y5).^2));

    plot(t,ye,'g-')

    xlabel('时间(t)');ylabel('气温(℃)');legend('气温散点图','拟合图');

    e2

     

    6.5 综合对比

    t=0:24;

    y=[15 14 14 14 14 15 16 18 20 22 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16];

    [p2,e2]=polyfit(t,y,2);[p3,e3]=polyfit(t,y,3);[p4,e4]=polyfit(t,y,4);

    p=polyfit(t,log(y),2)

    b=-p(1)

    c=p(2)/(2*b)

    a=exp(p(3)+b*c^2)

    x1=0:0.1:24;

    y2=polyval(p2,x1);

    y3=polyval(p3,x1);

    y4=polyval(p4,x1);

    ye=a*exp(-b*(x1-c).^2);y5=a*exp(-b*(t-c).^2);e5=sqrt(sum((y-y5).^2));

    plot(t,y,'*',x1,y2,'r--',x1,y3,'b-',x1,y4,'k-',x1,ye,'g-')

    xlabel('时间(t)');ylabel('气温(℃)');legend('气温散点图','拟合图');

    title('二次函数');

    e2

    展开全文
  • 最小二乘拟合结语 前言 Hello!小伙伴! 非常感谢您阅读海轰的文章,倘若文中有错误的地方,欢迎您指出~   自我介绍 ଘ(੭ˊᵕˋ)੭ 昵称:海轰 标签:程序猿|C++选手|学生 简介:因C语言结识编程,随后...

    前言

    Hello!小伙伴!
    非常感谢您阅读海轰的文章,倘若文中有错误的地方,欢迎您指出~
     
    自我介绍 ଘ(੭ˊᵕˋ)੭
    昵称:海轰
    标签:程序猿|C++选手|学生
    简介:因C语言结识编程,随后转入计算机专业,有幸拿过一些国奖、省奖…已保研。目前正在学习C++/Linux/Python
    学习经验:扎实基础 + 多做笔记 + 多敲代码 + 多思考 + 学好英语!
     
    初学Python 小白阶段
    文章仅作为自己的学习笔记 用于知识体系建立以及复习
    题不在多 学一题 懂一题
    知其然 知其所以然!

    往期文章

    Python数学建模系列(一):规划问题之线性规划

    Python数学建模系列(二):规划问题之整数规划

    Python数学建模系列(三):规划问题之非线性规划

    1. 一维插值

    插值:求过已知有限个数据点的近似函数。

    插值函数经过样本点,拟合函数一般基于最小二乘法尽量靠近所有样本并穿过。常见差值方法有拉格朗日插值法、分段插值法、样条插值法。
    image.png

    interp1d(x, y) 计算一维插值

    ​1.1 线性插值与样条插值(B-spline)

    例1:某电学元件的电压数据记录在0~2.25πA范围与电流关系满足正弦函数,分别用线性插值和样条插值方法给出经过数据点的数值逼近函数曲线。

    Demo代码

    import matplotlib
    import numpy as np
    from scipy import interpolate
    import matplotlib.pyplot as plt
    
    # 引入中文字体 
    font = {
        "family": "Microsoft YaHei"
    }
    matplotlib.rc("font", **font)
    
    # 初始数据量 0 - 2.25pi 分为10份 均匀分
    x = np.linspace(0, 2.25 * np.pi, 10)
    y = np.sin(x)
    
    # 得到差值函数 (使用线性插值)
    f_linear = interpolate.interp1d(x, y)
    # 新数据  0 - 2.25pi 分为100份 均匀分 (线性插值)
    x_new = np.linspace(0, 2.25 * np.pi, 100)
    y_new = f_linear(x_new)
    
    # 使用B-spline插值
    tck = interpolate.splrep(x, y)
    y_bspline = interpolate.splev(x_new, tck)
    
    # 可视化
    plt.xlabel(u'安培/A')
    plt.ylabel(u'伏特/V')
    plt.plot(x, y, "o", label=u"原始数据")
    plt.plot(x_new, f_linear(x_new), label=u"线性插值")
    plt.plot(x_new, y_bspline, label=u"B-spline插值")
    plt.legend()
    plt.show()
    

    输出:
    image.png
    涉及知识点:

    • numpy.linspace
    • scipy.interpolate.interp1d
    • scipy.interpolate.splrep

    1.2 高阶样条插值

    随着插值节点增多,多项式次数也变高,插值曲线在一些区域出现跳跃,并且越来越偏离原始曲线,称为龙格现象。

    例2:某电学元件的电压数据记录在010A范围与电流关系满足正弦函数,分别用05阶样条插值方法给出经过数据点的数值逼近函数曲线。

    Demo代码

    import matplotlib
    import numpy as np
    from matplotlib import pyplot as plt
    from scipy import interpolate
     
    font = {
        "family": "Microsoft YaHei"
    }
    matplotlib.rc("font", **font)
    # 创建数据点集
    x = np.linspace(0, 10, 11)
    y = np.sin(x)
    # 绘制数据点集
    plt.figure(figsize=(12, 9))
    plt.plot(x, y, 'ro')
    # 根据kind创建interp1d对象f、计算插值结果
    xnew = np.linspace(0, 10, 101)
    # 邻接 0阶 线性 二阶
    for kind in ['nearest', 'zero', 'linear', 'quadratic']:
        f = interpolate.interp1d(x, y, kind=kind)
        ynew = f(xnew)
        plt.plot(xnew, ynew, label=str(kind))
    plt.xticks(fontsize=20)
    plt.yticks(fontsize=20)
    plt.legend(loc="lower right")
    plt.show()
    

    输出:
    image.png

    分别对每一种插值方式进行查看

    1.当kind = nearest时

    import matplotlib
    import numpy as np
    from matplotlib import pyplot as plt
    from scipy import interpolate
     
    font = {
        "family": "Microsoft YaHei"
    }
    matplotlib.rc("font", **font)
    
    # 创建数据点集
    x = np.linspace(0, 10, 11)
    y = np.sin(x)
    
    # 得插值函数
    f = interpolate.interp1d(x, y, kind='nearest')
    
    # 新数据
    x_new = np.linspace(0,10,101)
    y_new = f(x_new)
    
    # 可视化
    plt.plot(x, y, 'o', x_new, y_new, '-')
    plt.show()
    

    image.png
    2.当kind = zero时

    import matplotlib
    import numpy as np
    from matplotlib import pyplot as plt
    from scipy import interpolate
     
    font = {
        "family": "Microsoft YaHei"
    }
    matplotlib.rc("font", **font)
    
    # 创建数据点集
    x = np.linspace(0, 10, 11)
    y = np.sin(x)
    
    # 得插值函数
    f = interpolate.interp1d(x, y, kind='zero')
    
    # 新数据
    x_new = np.linspace(0,10,101)
    y_new = f(x_new)
    
    # 可视化
    plt.plot(x, y, 'o', x_new, y_new, '-')
    plt.show()
    

    image.png
    3.当kind = linear时

    import matplotlib
    import numpy as np
    from matplotlib import pyplot as plt
    from scipy import interpolate
     
    font = {
        "family": "Microsoft YaHei"
    }
    matplotlib.rc("font", **font)
    
    # 创建数据点集
    x = np.linspace(0, 10, 11)
    y = np.sin(x)
    
    # 得插值函数
    f = interpolate.interp1d(x, y, kind='linear')
    
    # 新数据
    x_new = np.linspace(0,10,101)
    y_new = f(x_new)
    
    # 可视化
    plt.plot(x, y, 'o', x_new, y_new, '-')
    plt.show()
    

    image.png
    4.当kind = quadratic时

    import matplotlib
    import numpy as np
    from matplotlib import pyplot as plt
    from scipy import interpolate
     
    font = {
        "family": "Microsoft YaHei"
    }
    matplotlib.rc("font", **font)
    
    # 创建数据点集
    x = np.linspace(0, 10, 11)
    y = np.sin(x)
    
    # 得插值函数
    f = interpolate.interp1d(x, y, kind='quadratic')
    
    # 新数据
    x_new = np.linspace(0,10,101)
    y_new = f(x_new)
    
    # 可视化
    plt.plot(x, y, 'o', x_new, y_new, '-')
    plt.show()
    

    image.png
    5.当kind = cubic时

    import matplotlib
    import numpy as np
    from matplotlib import pyplot as plt
    from scipy import interpolate
     
    font = {
        "family": "Microsoft YaHei"
    }
    matplotlib.rc("font", **font)
    
    # 创建数据点集
    x = np.linspace(0, 10, 11)
    y = np.sin(x)
    
    # 得插值函数
    f = interpolate.interp1d(x, y, kind='cubic')
    
    # 新数据
    x_new = np.linspace(0,10,101)
    y_new = f(x_new)
    
    # 可视化
    plt.plot(x, y, 'o', x_new, y_new, '-')
    plt.show()
    

    image.png

    2. 二维插值

    interp2d(x, y, z, kind=“’’) 计算二维插值

    2.1 图像模糊处理——样条插值

    例3:某图像表达式为 z = ( x + y ) e − 5 ∗ ( x 2 + y 2 ) z = (x+y)e^{-5*(x^2 + y^2)} z=(x+y)e5(x2+y2),完成图像的二维插值使其变清晰。

    Demo代码

    import numpy as np
    from scipy import interpolate
    import matplotlib.pyplot as plt
     
     
    def func(x, y):
        return (x + y) * np.exp(-5.0 * (x ** 2 + y ** 2))
    
    # X-Y轴分为15*15的网格
    # x, y = np.mgrid[-1:1:15j, -1:1:15j]
    x = np.linspace(-1, 1, 15)
    y = np.linspace(-1, 1, 15)
    x, y = np.meshgrid(x, y)
    fvals = func(x, y)
    
    # 二维插值
    newfunc = interpolate.interp2d(x, y, fvals, kind='cubic')
    
    # 计算100*100网格上插值
    xnew = np.linspace(-1, 1, 100)
    ynew = np.linspace(-1, 1, 100)
    fnew = newfunc(xnew, ynew)
    xnew, ynew = np.meshgrid(xnew, ynew)
    
    plt.subplot(121)
    # extent x轴和y轴范围
    im1 = plt.imshow(fvals, extent=[-1, 1, -1, 1], interpolation="nearest", origin="lower",cmap="Reds")
    plt.colorbar(im1)
    
    plt.subplot(122)
    im2 = plt.imshow(fnew, extent=[-1, 1, -1, 1], interpolation="nearest", origin="lower",cmap="Reds")
    plt.colorbar(im2)
    plt.show()
    

    输出:
    在这里插入图片描述

    2.2 二维插值的三维图

    例4:某图像表达式为 z = ( x + y ) e − 5 ∗ ( x 2 + y 2 ) z = (x+y)e^{-5*(x^2 + y^2)} z=(x+y)e5(x2+y2),完成三维图像的二维插值可视化。

    其实就是在3二维插值基础上 实现了三维图像的绘制

    Demo代码

    import numpy as np
    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib as mpl
    from scipy import interpolate
    import matplotlib.cm as cm
    import matplotlib.pyplot as plt
    
    def func(x, y):
        return (x + y) * np.exp(-5.0 * (x ** 2 + y ** 2))
    
    # X-Y轴分为20*20的网格
    x = np.linspace(-1, 1, 20)
    y = np.linspace(-1, 1, 20)
    x, y = np.meshgrid(x, y)
    fvals = func(x, y)
    
    # 绘制分图1
    fig = plt.figure(figsize=(9, 6))
    ax = plt.subplot(1, 2, 1, projection='3d')
    surf = ax.plot_surface(x, y, fvals, rstride=2, cstride=2, cmap=cm.coolwarm, linewidth=0.5, antialiased=True)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('f(x,y)')
    plt.colorbar(surf, shrink=0.5, aspect=5)  # 添加颜色条标注
    
    # 二维插值
    newfunc = interpolate.interp2d(x, y, fvals, kind='cubic')
    
    # 计算100*100网格上插值
    xnew = np.linspace(-1, 1, 100)
    ynew = np.linspace(-1, 1, 100)
    fnew = newfunc(xnew, ynew)
    xnew, ynew = np.meshgrid(xnew, ynew)
    
    ax2 = plt.subplot(1, 2, 2, projection='3d')
    surf2 = ax2.plot_surface(xnew, ynew, fnew, rstride=2, cstride=2, cmap=cm.coolwarm, linewidth=0.5, antialiased=True)
    ax2.set_xlabel('xnew')
    ax2.set_ylabel('ynew')
    ax2.set_zlabel('fnew(x,y)')
    plt.colorbar(surf2, shrink=0.5, aspect=5)
    
    # 标注
    plt.show()
    

    输出:
    image.png

    3. 最小二乘拟合

    拟合指的是已知某函数的若干离散函数值{f1,f2,…,fn},通过调整该函 数中若干待定系数f(λ1, λ2,…,λn),使得该函数与已知点集的差别(最小二乘意义)最小。

    如果待定函数是线性,就叫线性拟合或者线性回归(主要在统计中),否则叫作非线性拟合或者非线性回归。表达式也可以是分段函数,这种情况下叫作样条拟合。

    从几何意义上讲,拟合是给定了空间中的一些点,找到一个已知形式、未知参数的连续曲面来最大限度地逼近这些点;而插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。

    选择参数c使得拟合模型与实际观测值在曲线拟合各点的残差(或离差)ek=yk-f(xk,c)的加权平方和达到最小,此时所求曲线称作在加权最小二乘意义下对数据的拟合曲线,这种方法叫做最小二乘法。

    涉及知识点

    from scipy.optimize import leastsq
    

    例5:对下列电学元件的电压电流记录结果进行最小二乘拟合,绘制相应曲线。
    电流(A)8.19 2.72 6.39 8.71 4.7 2.66 3.78
    电压(V)7.01 2.78 6.47 6.71 4.1 4.23 4.05

    在这里插入图片描述

    Demo代码

    import matplotlib
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.optimize import leastsq
     
    # 引入中文字体
    font = {
        "family": "Microsoft YaHei"
    }
    matplotlib.rc("font", **font)
    
    # 设置图字号
    plt.figure(figsize=(9, 9))
    
    # 初始数据值
    X = np.array([8.19, 2.72, 6.39, 8.71, 4.7, 2.66, 3.78])
    Y = np.array([7.01, 2.78, 6.47, 6.71, 4.1, 4.23, 4.05])
     
    # 计算以p为参数的直线与原始数据之间误差
    def f(p):
        k, b = p
        return (Y - (k * X + b))
     
    # leastsq使得f的输出数组的平方和最小,参数初始值k、b设为[1,0]
    r = leastsq(f, [1, 0])
    
    # 得到计算出的最优k、b
    k, b = r[0]
    
    # 可视化
    plt.scatter(X, Y, s=100, alpha=1.0, marker='o', label=u'数据点')
    x = np.linspace(0, 10, 1000)
    y = k * x + b
    ax = plt.gca()
    plt.plot(x, y, color='r', linewidth=5, linestyle=":", markersize=20, label=u'拟合曲线')
    plt.legend(loc=0, numpoints=1)
    leg = plt.gca().get_legend()
    ltext = leg.get_texts()
    plt.setp(ltext, fontsize='xx-large')
    plt.xlabel(u'安培/A')
    plt.ylabel(u'伏特/V')
    plt.xlim(0, x.max() * 1.1)
    plt.ylim(0, y.max() * 1.1)
    plt.xticks(fontsize=20)
    plt.yticks(fontsize=20)
    plt.legend(loc='upper left')
    plt.show()
    

    输出:
    image.png

    结语

    学习来源:B站及其课堂PPT,对其中代码进行了复现

    链接:https://www.bilibili.com/video/BV12h411d7Dm?from=search&seid=5685064698782810720

    文章仅作为学习笔记,记录从0到1的一个过程

    希望对您有所帮助,如有错误欢迎小伙伴指正~

    我是 海轰ଘ(੭ˊᵕˋ)੭

    如果您觉得写得可以的话,请点个赞吧

    谢谢支持 ❤️

    在这里插入图片描述

    展开全文
  • 博士研究高等数值分析试题.doc
  • 有效数字、数值修约及运算法则.ppt
  • 最小二乘法—多项式拟合非线性函数

    万次阅读 多人点赞 2019-07-10 15:02:57
    二、误差函数—最小二乘法 三、引出案例函数曲线 、目标函数 五、优化目标函数 六、优化目标函数—梯度下降法 七、优化目标函数—求解线性方程组 八、python编程实现拟合曲线函数 九、结果分析 一、函数...

    原文地址:https://www.jianshu.com/p/af0a4f71c05a

    目录

    一、函数的近似表示—高次多项式

    二、误差函数—最小二乘法

    三、引出案例函数曲线

    四、目标函数

    五、优化目标函数

    六、优化目标函数—梯度下降法

    七、优化目标函数—求解线性方程组

    八、python编程实现拟合曲线函数

    九、结果分析


     

    一、函数的近似表示—高次多项式

    为了研究一些复杂的函数,我们希望对函数自变量进行有限的加、减、乘法三种算数运算,便可以求出原函数值,因此我们通常使用高次多项式来近似表示函数

                                                   高次多项式近似表示f(x)

    可以看到,我们可以用n阶多项式的系数线性组合来近似表示f(x)

    特别说明,由泰勒展开式可知,当pn(x)的各阶导数和f(x)的各阶导数都相等,则f(x)和pn(x)的误差只是(x-x0)^n的高阶无穷小

                                                              泰勒公式

    二、误差函数—最小二乘法

    我们用f(x)的真实函数值减去多项式函数的结果的平方,来表示f(x)和多项式函数的误差关系,即

                                                      最小二乘法表示误差

    我们令x0=0,则最小二乘法表示的误差为

                                                   最小二乘法表示误差

    三、引出案例函数曲线

    有了上述数学知识,下面我们用一个案例来介绍最小二乘法拟合非线性函数的算法步骤

    案例:求一个多项式来拟合下列函数的函数值

                                                                    案例函数

    其中我们加入了噪点数据noise,使得函数在定义域内随机的震荡

    我们画出f(x)在[-1, 1]之间的图像,即

                                                                                                  案例函数图像

    案例问题即为:已知上述N个数据点,求其函数f(x)的表达式?

    显然,f(x)也就是机器学习要学习的目标

    四、目标函数

    下面我们开始推导机器学习的过程

    机器学习的目标为:

    (1) 学习一个f(x)多项式,可以拟合真实数据的变化趋势

    (2)f(x)的目标:使每一个真实数值到f(x)的拟合数值的距离之和最小

    翻译为数学语言,即

                             学习到的f(x)

                                                   目标函数

    五、优化目标函数

    为了使目标函数最优,我们对每个系数求其偏导数为0,即

                                                     优化目标函数

    至此,我们需要根据上述k的等式,求出所有的系数ak,有两种方法可以求解

    (1)梯度下降法

    (2)求解线性方程

    六、优化目标函数—梯度下降法

    采用梯度下降法,可以避免我们用数学方法直接一步来求解上述k个方程的最优解,而是采用迭代逼近最优解的思路,其具体步骤为:

    (1)初始化k个系数值,开始迭代

    (2)每次迭代,分别求出各个系数ak对应的梯度值

    (3)用梯度值和学习率来更新每一个系数ak

    (4)保证每次更新完所有系数,对应的损失函数值都在减少(说明梯度一直在下降)

    翻译为数学语言为:

    计算ak对应的梯度值

                                                     计算ak的梯度值

    更新ak

                                  更新ak

    七、优化目标函数—求解线性方程组

    求解线性方程组让我们直接一步求出偏导数最优解,其具体步骤为:

    (1)将最优化偏导数的方程组写为矩阵乘法形式:XA=Y

    (2)利用数学消元算法来求解XA=Y

    我们对k个偏导数=0的等式进行化简处理,即

                                                        k个偏导数等式

    下面我们将其写为矩阵相乘的形式,将所有只含有x的系数写为第一个矩阵X,即

                                                 含有x的系数矩阵

    将只含有待求解的多项式系数ak写为第二个矩阵A,即

          含有a的系数矩阵

    最后将含有y的系数写为第三个矩阵Y,即

                  含有y的系数矩阵

    此时,我们就可以用矩阵的乘法来表示最初的k个偏导数为0的等式,即

             k个偏导数等式的矩阵形式

    注意:其中X是k*k的矩阵,A是k*1的矩阵,Y是k*1的矩阵,符合矩阵的乘法规则

    从矩阵表达式可以看出,X和Y矩阵是已知的,A矩阵只包含待求解的f(x)的所有多项式系数,则问题即转化为线性方程组:XA=Y

    求解线性方程组,大概有如下算法

    (1)高斯消元法(全主元、列主消元)

    (2)LU三角分解法

    (3)雅克比迭代法

    (4)逐次超松弛(SOR)迭代法

    (5)高斯-赛德尔迭代法

    这里我们采用高斯消元法来求解A,具体算法步骤较为复杂,我们在单独的章节介绍高斯消元算法

    八、python编程实现拟合曲线函数

    生成带有噪点的待拟合的数据集合

                                                                           生成带有噪点的待拟合的数据集合

    计算最小二乘法当前的误差

                                                                                计算最小二乘法当前的误差

    迭代解法:最小二乘法+梯度下降法

                                                                           迭代解法:最小二乘法+梯度下降法

    数学解法:最小二乘法+求解线性方程组

                                                                           数学解法:最小二乘法+求解线性方程组

    九、结果分析

    我们来拟合N=10次幂的多项式,分别用梯度下降法和求解线性方程组的算法来比较拟合结果

    (1)使用梯度下降法来优化目标函数,其拟合结果为:

                                                                                         梯度下降法拟合结果

    其误差变化为

                                                                                               梯度下降法优化的误差

    (2)使用求解线性方程组来优化目标函数,其拟合结果为:

                                                                                        求解线性方程组拟合结果

    其误差变化为

                               求解线性方程组优化的误差

    比较上述两种解法,以及多项式拟合的思想,我们可以总结出:

    • 利用函数的近似多项式表示,和最小二乘法来表示目标函数,我们可以高效的拟合非线性函数
    • 梯度下降法可以一步步迭代逼近目标函数的极值,而不用直接求出偏导数最优解
    • 求解方程组,使用数学消元的算法,直接一步优化完成了目标函数的偏导数最优解
    • 从结果上比较,求解方程组的效率和拟合结果要比梯度下降法好

    案例代码见:最小二乘法—多项式拟合非线性函数

    示例论文预测:基于非线性最小二乘曲线拟合法的电子商务预测模型

    展开全文
  • 提供了数值分析实验的八实验,可供参考。实验一 函数插值方法;实验二 函数逼近与曲线拟合;实验三 数值积分与数值微分;实验 线方程组的直接解法等等
  • 数据库最小的存储单位是

    千次阅读 2021-02-03 01:16:51
    【简答题】读48元辅音 (10.0分) 【简答题】笔记 (5.0分) 【单选题】段是表空间中一种逻辑存储结构,以下( )是 ORACLE 数据库使用的段类型 【简答题】道路交通安全法对乘车人作了哪些禁止性的规定? 【简答题】道路...
  • python最小二乘法拟合实例

    千次阅读 多人点赞 2019-11-05 14:50:04
    最小二乘法主要解决的是插值中出现的三问题: ① 高多项式会龙格现象; ② 定函数关系不合理(数据误差); ⑶评估逼近效果: 假设“逼近”规律的近似函数为y=f(x), 即有yi﹡= f(xi)(i=1,2,…,m).它与观测值yi之差...
  • 本文代码实现基本按照《数据结构》课本目录顺序,外加大量的复杂算法实现,一篇文章足够。能换你一收藏了吧?
  • 例题:下列给定程序中,函数fun的功能是:从str所指字符串中,找出s所指子串的个数作为函数值返回。 例如,当str 所指字符串中的内容为asdfghasdfgh,s所指字符串的内容为as,则函数返回整数2。 注意:不要改动main...
  • 例题:给定程序中,函数fun 的功能是:将形参student所指结构体数组中年龄最小者的数据作为函数值返回,并在main函数中输出。 注意:请勿改动主函数main与其它函数中的任何内容,仅在fun函数的横线上填写所需的若干...
  • 例题:下列给定程序中,函数fun的功能是:统计一无符合整数中各位数字值为0的个数,通过形参传回主函数,并把该整数中各位上最大的数字值作为函数值返回。 例如。若输入无符号整数10080,则数字值为0的个数为3,...
  • 最小费用流及其求法

    万次阅读 多人点赞 2019-05-04 08:33:48
    【3】树:基本概念与最小生成树 【4】匹配问题: 匈牙利算法 、最优指派、相等子图 【5】Euler 图和 Hamilton 图 【6】计划评审方法和关键路线法【统筹方法】:广泛地用于系统分析和项 目管理 【7】最小费用流...
  • 1、一个四位二进制补码的表示范围是( B ) A、0~15 B、-8~7 C、-7~7 D、-7~8 过程:二进制补码取值范围为. 2、十进制数-48 用补码表示为( B ) A、10110000 B、11010000 C、11110000 D、11001111 过程:-48 的...
  • 超硬核!数据结构学霸笔记,考试面试吹牛就靠它

    万次阅读 多人点赞 2021-03-26 11:11:21
    2)四个并列程序段:分别为O(N),O(N^2),O(N^3),O(N*logN),整个程序时间复杂度为O(N^3),因为随着N的增长,其它项都会忽略 3)算法分析的目的是分析算法的效率以求改进 4)对算法分析的前提是算法必须正确 5)原地...
  • 简言之,回归最简单的定义就是: 给出一点集,构造一函数来拟合这点集,并且尽可能的让该点集与拟合函数间的误差最小,如果这函数曲线是一条直线,那就被称为线性回归,如果曲线是一条三次曲线,就被称为三...
  • A、既不能移动,也不能改变大小 B、仅可以移动,不能改变大小 C、仅可以改变大小,不能移动 D、既能移动,也能改变大小 8、下列四个不同数制表示的数中,数值最大的是( )。 A、二进制数11011101 B、八进制数334 C...
  • 19计算机应用基础.doc

    2022-07-03 09:34:13
    《计算机应用基础》统考练习题三 1.计算机中存储数据的最小...这组数中,三个数值相同是 。 A.277,10111111,BF B.203,10000011,83 C.247,10100111,A8 D.213,10010110,96 9.二进制数0.101转换成十进制数是 。 A. 0.
  • 例题:下列给定程序中,函数fun的功能是:实现两整数的交换。 例如,给x和y分别输入60和65,输出为:x=65 y=60。 注意:不要改动main函数,不能增行或删行,也不能更改程序的结构。 代码如下: #include<stdio....
  • 从上面的公式可以看到,因为有求和,参数的更新是基于所有的样本,其能得到全局最优解,求解的参数满足风险函数最小化。但是这种参数更新存在的问题是如果样本数量很大,会导致性能很差,这种更新方式叫做批量梯度...
  • A、运算器 B、控制器 C、存储器 D、中央处理器 4、键盘打字键区的8基准键位是指( )。 A、SDFGHJKL B、ASDFGHJK C、ASDFHJKL D、ASDFJKL; 5、下列设备属于外部设备的是( )。 A、CPU B、主板 C、...
  • 【判断题】查看变量类型的Python内置函数是type____【判断题】Python内置函数sum____用来返回数值型序列中所有元素之和。【单选题】Python中,用于获取用户输入的命令为:【多选题】以下关于机器学...
  • 青少年软件编程(Scratch)等级考试试卷(级)2020.9 分数:100.00 题数:30 一、单选题(共15题,每题2分,共30分) 1.执行下面程序,输入4和7后,角色说出的内容是?( ) A、4,7 B、7,7 C、7,4 D...
  • 完整的计算机系统包括

    千次阅读 2021-07-18 03:26:35
    --------------------------------------------------正文内容开始----... 一完整的计算机系统包括A). 主机、键盘和显示器 B). 计算机与外部设置C). 硬件系统和软件系统 D). 系统软件与应用软件2). 在多媒体计算...
  • 计算机网络(第7版) - 第章 网络层 - 习题

    万次阅读 多人点赞 2019-07-28 12:54:35
    网际协议TCP、IP是TCP/IP体系中两个最主要的协议之一,与IP协议配套使用的还有四个协议。  ARP协议:是解决同一个局域网上的主机或路由器的IP地址和硬件地址的映射问题。 RARP协议:是解决同一个局域网上的主机或...
  • 深度学习之数学基础(数值计算)

    万次阅读 2017-10-26 23:27:53
    信息论是应用数学的一分支,主要研究的是对一信号能够提供信息的多少进行量化。如果说概率使我们能够做出不确定性的陈述以及在不确定性存在的情况下进行推理,那信息论就是使我们能够量化概率分布中不确定性的...
  • 例题:请补充main函数,该函数的功能是:计算四个学生各科的平均分。 例如,当score[N][M]={{83,65,63},{89,93,95},{90,63,80},{56,75,77}}时,则平均分为:79.5,74.0,78.8 注意:仅在横线上填写所需的若干表达式...
  • 万字文肝Java基础知识(一)

    千次阅读 多人点赞 2021-07-05 18:59:27
    目录一、 前言、入门程序、常量1.1 Java ...软件的安装和配置三、注释和关键字3.1 注释3.2 关键字3.3 标识符、常量与变量4.1 常量的概念和分类4.2 打印不同类型的常量4.3 变量和数据类型【重要】4.4 数据类型转换4.

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 34,200
精华内容 13,680
关键字:

下列四个数值最小