最小二乘法曲线拟合

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

最小二乘法曲线拟合 篇1

在科学实验和生产实践中, 我们常常需要从一组测定的数据 (例如N个点 (xi, yi) ) 去求得自变量x和因变量y的一个近似解析表达式y=φ (x) , 最小二乘法是经常使用的方法.最小二乘法适用于做线性拟合, 此时拟合曲线的一般式为y=a+bx, 而在许多实际问题中, 变量之间内在的关系呈非线性关系, 此时我们可以把拟合曲线y=a+bx中的自变量x和因变量y看成是其他变量的函数, 例如原来的函数关系是f (y′) =a+bg (x′) , 令x=g (x′) , y=f (y′) , 则经过变换后即得y=a+bx的形式, 于是许多非线性问题就转换成了线性问题而得到解决.

如何找到更符合实际情况的数据拟合, 一方面要根据专业知识和经验来确定经验曲线的近似公式, 另一方面要根据散点图的分布形状及特点来选择适当的曲线拟合这些数据.

例1在某化学反应里, 测得生成物浓度y (%) 与时间t (min) 的数据见下表, 试用最小二乘法建立x与y之间的经验公式.

解将已知数据点 (ti, yi) (i=1, 2, …, 16) 描绘在坐标纸上 (图略) , 观察可知拟合曲线y=φ (t) 应具有如下特点:

(1) 曲线随着t的增加而上升, 但上升速度由快到慢.

(2) 当t=0时, 反应尚未开始, 即y=0;当t→∞时, y趋于某一常数, 故曲线通过原点且有一水平渐近线.

具有上述特点的曲线很多, 选用不同的数学模型, 可以获得不同的拟合曲线和经验公式.

把这两个经验公式进行比较:

从均方误差和最大偏差两个不同的角度看, 后者优于前者.因此在解决实际问题时, 常常要经过反复分析, 多次选择、计算和比较, 才能获得较好的数学模型.

最小二乘法的应用很广, 下面我们利用最小二乘法来解决一个实际问题.

例2用电压U=10伏的电池给电容器充电, 电容器上t时刻的电压为, 其中U0是电容器的初始电压, τ是充电常数, 试由下面一组数据确定U0和τ.

解电容器上t时刻的电压为:U (t) =U- (U-U0) ·

式变为y=a+bx, 求解得a=1.4912, b=-0.2857, 进而求得V0=5.5577, τ=3.5002.

Matlab程序如下:

最小二乘法曲线拟合 篇2

工程施工中, 我们会经常取得一些相关的数据, 这些数据往往来自与施工密切相关的测量或试验中, 比如石灰剂量的滴定试验中, EDTA用量与灰土剂量存在某种关系;又如预应力千斤顶与油表的配套校验中, 油表读数与千斤顶实际张拉力又有一种关系, 这些原始数据一般是5组以上。

我们可以通过作图或多段插值取得变量之间的联系, 但作图和插值查图往往误差较大, 这时可采用最小二乘法先拟合出一个多项式 (也称经验公式) , 再根据此多项式求解任一自变量所对应的因变量较精确的结果, 据此绘图可得到较精确、较合理的曲线。

2 最小二乘法原理

现实中通过测量或试验取得的各组数据 (xi, yi) 其本身不可避免的地带有测试误差, 如果构造一个较为简单的插值法P (x) 来逼近真实函数f (x) , 当个别数对误差影响较大时就会引起插值函数发生严重波动, 从而影响逼近精度, 因为插值法要求插值函数在节点 (已知数对) 处满足标准条件, 即P (xj) =f (xj) , (j=0, 1, …, n) 。这时候, 为尽可能减小测试误差对逼近精度的影响, 我们可以用另一种方法构造一个经验公式, 使得该公式在每一个节点上所求得的结果与原测试结果的差的平方和最小, 即曲线拟合的误差最小, 精度最高, 这就是最小二乘法原理, 用定义表述为:

设有n对数据 (xi、yi) , (j=0, 1, …, n) , 通过这些数据找一个m次近似多项式P (x) =a0+a1x+…+amxm (m

使得为最小值,

则称P (x) 为最小二乘拟合多项式, 或称x、y之间的经验公式。

3 最小二乘原理应用的具体做法

根据定义, 我们最终目的是根据已知数据组求出合适的系数ai, 使多项式φ (ai) 取得最小值, 运用多元函数求极值的方法可导出如下方程组:

这是以ai为未知数的m+1阶线性方程组, 写成矩阵形式有

因此, 我们只要求出S0, S1, S2, …, Sm, …, S2m以及T0、T1、…、Tm, 建立正规方程组即可求出系数ai, 代入多项式P (x) 即得拟合曲线方程。

4 应用实例

4.1 灰剂量EDTA滴定法标准曲线绘制备5种试样共10个样品的灰土混合料, 其石灰剂量为已知, 在6%~18%范围内按3%剂量递增, 该范围已涵盖2:8灰土到3:7灰土 (V:V) 的石灰剂量范围。分别对样品进行滴定, 测得EDTA标准液耗量 (ml) 如下:

步骤1:作草图, 先将数据点画在坐标纸上, 连线分析其变化的粗糙趋向, 确定多项式的大致形式, 可以看出, 本例为一缓和抛物线。

步骤2:设拟合曲线方程为y=a+bx+cx2, 列表求和计算Sk、Tk (见下表)

步骤3:建立正规方程组

解得c=-0.0013b=0.382a=0.253

因此得方程y=0.253+0.382x-0.0013x2

该经验公式求得的结果与测试数据的最大离差Δ=-0.2, 离差的平方和最小为0.06, 据此可绘制石灰剂量 (%) 与EDTA耗量 (ml) 关系的标准曲线, 应用中当已知某灰土样品液滴定的EDTA耗量x时, 直接代入公式或查标准曲线便可得到较为精确的待测灰土的石灰剂量, 方便、快速而又准确。

4.2 求解张拉千斤顶与油表读数的回归方程预应力千斤顶与油表的配套校验中, 分级张拉数据可达到5~20组, 而张拉力与油表读数实际为线性关系, 一般只需两组数据便可确定其关系式, 但数据越多, 回归方程越真实, 越精确。此时采用最小二乘法可使每一组数据参与回归。

对线性方程y=ax+b

用最小二乘法公式求系数a、b (斜率和截距) , 公式简化为:

某千斤顶校验数据见下表 (n=7)

由于第一数据 (0, 0) 将可能使计算失真, 误差较大, 回归时舍去该组数据, 取n=6组, 计算得:

则有a=0.035 b=0.68,

故得经验公式y=0.035x+0.68

回归相关系数r=0.99993, 精度较高。应用时将所需张拉力的值代入公式便可得到相应张拉力下的油表读数, 以控制张拉力。

5 结语

最小二乘法作为函数逼近的一种重要方法, 在工程技术中的曲线拟合、求取经验公式等方面有着广泛的应用, 但一般来说拟合式不高于三次多项式, 否则在计算上将十分繁琐, 且拟合精度将受到影响甚至出现病态。而实际工程应用中拟合的曲线多为一次线性或二次抛物线型, 因此, 最小二乘法拟合曲线在工程应用中仍占有着重要的位置。

摘要：施工中常遇到一些相关数据的分析处理, 并要求拟合曲线以便反映数组规律和扩大应用范围, 本文介绍了目前应用较广、较为精确的曲线拟合方法, 以及回归方程的建立。

关键词：最小二乘法,拟合曲线,施工,应用

参考文献

最小二乘法曲线拟合 篇3

传统的分段曲线拟合根据主观经验和绘制数据散点图来确定拟合的经验函数和分段点。文献[5] 提出分段区间重合的拟合方法, 由每4个数据点决定一个三次曲线, 但分段区间太密, 不适用于密集的数据拟合。文献[6]提出的多项式基函数的全局连续拟合方法, 只限于2个分段区间。文献[7]提出多分段区间全局连续的曲线拟合方法, 但基函数只限于一次多项式。文献[5—7]提出的方法都是根据主观经验来确定经验函数和分段区间, 然后进行拟合。文献[8]引入拟合误差限度 εmax, 在误差大于该限度的点处分段, 但该限度 εmax不容易界定。

提出一种自动分段多项式曲线拟合方法, 该方法在实际运用中具有如下有点: 1提供几种不同的经验函数, 根据不同经验函数拟合的数据和实测数据的误差的方差, 自动确定较优的经验函数; 2自动确定数据的分段区间, 在各个区间进行最小二乘拟合。通过数值实验, 证明该方法有较高的拟合精度。

1经验函数的确定及拟合步骤

1. 1经验函数的确定

数据拟合方法有很多, 例如对数曲线拟合, 反函数曲线拟合, 二次曲线拟合, 三次曲线拟合, 幂函数曲线拟合, 指数曲线拟合等。一般先观察散点图来确定曲线的类型, 不过散点图都是相关关系的粗略表示, 有时候散点图可能与几种曲线都很接近, 这时建立相应的经验函数可能都是合理的, 但由于选择不同的曲线, 得到同一个问题的多个不同经验函数, 怎样从这些经验函数中选择最优的一个。文献[1] 提出的, 用几种函数进行拟合, 计算历史数据点实测值和拟合值的误差平方和最小的为最优经验函数这一方法, 可能存在这样的问题: 误差平方和最小, 但误差波动较大, 即一些点误差很小, 一些点误差相对较大。针对这种情况, 本文提出一种新的确定经验函数的方法, 用几种不同的函数进行拟合, 从中选取最优的经验函数, 最优经验函数确定的条件如下:

1) 历史数据点误差为正和误差为负的个数之差小于适应性参数, 在本文试验中选用的适应性参数为3。

2) 计算误差的方差, 方差最小的为最优的拟合函数。

方差是各个数据与平均数之差的平方和的平均数。通俗点讲, 就是各点与平均数偏离的程度, 来衡量一批数据相较于平均数的波动的大小。方差的数学描述为

式 ( 1) 中把x－视为第一个历史点的误差, xi为后面各历史点的误差, 该公式即可理解为各点误差相较于第一个点的误差的波动大小, 波动越小, 说明各个点的误差相差越小, 分布较均匀, 那么该函数即为最优经验函数。经验函数确定的流程图如图1所示。

1. 2自动分段拟合的步骤

本文选取三种较为典型的回归方程类型

分段拟合步骤如下, 分段拟合流程图如图2所示。

1) 根据1. 1节提出的方法拟合所有的历史数据, 在上述的三种拟合函数中选择最优的拟合函数。

2) 由第1) 步选取的经验函数, 计算历史数据点拟合值和实测值的误差, 计算历史数据点误差绝对值的均值S。

3) 比较各历史点误差与误差均值S, 若连续3个点的误差绝对值大于均值S, 则从第一个误差大于均值的点处分段, 执行第4) 步, 若不存在拟合值和实测值误差大于均值的点, 或者这种点只有一个或两个, 则不进行分段, 执行第5) 步。

4) 从分段点处到历史数据最后一个点重新拟合, 选择最优的函数, 重复以上步骤。

5) 根据以上步骤得到的分段, 以及各段选出的最优经验函数对历史数据进行拟合。

2数值实验

2. 1自动分段拟合实验

为检验所论述的自动分段多项式曲线拟合方法, 采用两组数据, 厦门高崎机场1990 ~ 2012年的起降架次的数据［9］和1990 ~ 2004年中国人均GDP的数据来进行拟合。

提出的自动分段多项式曲线拟合方法对厦门高崎机场1990 ~ 2012年起降架次统计数据和1990 ~ 2004年中国人均GDP数据进行拟合, 用JAVA编程实现, JFree Chart绘制曲线。表1和图3是厦门高崎机场起降架次的数据与经过分段曲线拟合后的数据的比较, 表2和图4是中国人均GDP数据与拟合后的数据的比较。

对于厦门高崎机场起降架次的数据, 提出的该方法自动将数据分为3段, 即数据1990 ~ 1994, 1994 ~ 2003, 2003 ~ 2012, 这三段数据选取的最优经验函数都是指数函数形式。对于中国人均GDP的数据, 程序自动将数据分为三段, 即数据段1990 ~ 1995, 1995 ~ 1999, 1999 ~ 2004, 这三段数据选取的最优经验函数分别是多项式函数形式, 直线形式和多项式函数形式。由表1、表2和图3、图4可以看出, 提出的自动分段拟合方法拟合精度较高, 能够很好的拟合历史数据。

2. 2对比实验

本文选取1990 ~ 2004年中国人均GDP数据, 分别用本文提出的自动分段拟合方法和传统直线拟合、指数拟合、多项式拟合方法对该组数据进行实验, 对比实验各种方法的相对误差如表3所示。

由表3可以看出, 自动分段拟合方法相对误差比传统拟合方法相对误差小, 且相对误差波动比传统方法波动小。自动分段拟合方法不仅在拟合精度上比传统方法高, 而且分段方便, 不需用户根据主观经验和绘制散点图来分段, 方便用户使用。

3结论

提出的自动分段曲线拟合方法采取对历史数据进行拟合, 自动选取较优的经验函数, 自动进行分段, 使得在进行数据拟合时, 不需要人为的绘制出历史数据散点图来选取经验函数和分段。通过编程实现和数值实验验证, 该方法不仅便于在计算机上编程实现, 而且拟合效果较好。

摘要：针对传统的分段曲线拟合方法在选择拟合函数和确定分段区间时经验成分较多的不足, 提出一种自动分段多项式曲线拟合方法, 根据误差方差和误差均值, 自动确定经验函数和分段区间。通过实际数据的检验, 验证了该方法的拟合效果。

关键词：数据拟合,分段拟合,多项式曲线,最小二乘法

参考文献

[1] 吴宗敏.散乱数据拟合的模型、方法和理论.北京:科学出版社, 2007

[2] 杨维, 张晓明.利用最小二乘法进行回归分析及经验公式的确定.沈阳工业大学学报, 1991;13 (2) :1—6

[3] 侯超钧, 曾艳姗, 吴东庆, 等.全局连续的分段最小二乘曲线拟合方法.重庆师范大学学报 (自然科学版) , 2011;28 (6) :44—48

[4] Grama L, Rusu C.Phase approximation by divide-and-conquer piecewise linear fitting of gain.International Symposium on Signals, Circuits and Systems, 2009:1—4

[5] 蔡山, 张浩, 陈洪辉, 等.基于最小二乘法的分段三次曲线拟合方法研究.科学技术与工程, 2007;7 (3) :352—355

[6] 刘晓莉, 陈春梅.基于最小二乘原理的分段曲线拟合方法.伊犁教育学院学报, 2004;17 (3) :132—134

[7] Gluss B.An alternative method for continuous line segment curve-fitting.Information and Control, 1964;7 (2) :200—206

[8] 王成山, 黄碧斌, 李鹏, 等.燃料电池复杂非线性静态特性模型简化方法.电力系统自动化, 2011;35 (7) :64—69

最小二乘法曲线拟合 篇4

MLS法是对传统加权最小二乘法的进一步推广, 克服了最小二乘法的不足, 突出了加权函数的紧支性, 据此建立的拟合函数, 具有局部近似的特点。MLS法的基本理论本文不再赘述详见参考文献[4]。本文针对具体二维和三维算例分析了MLS的局部近似性, 为其他研究领域中引入该方法解决复杂问题提供了一种新的思路。

2 算法流程

结合移动最小二乘法的特点, 给出了基于MATLAB软件编写的MLS法的曲线曲面拟合程序, 程序设计流程如图1所示。

3 权函数的影响分析

由图2所示:权函数为恒定常数, 样本点的权函数值均为1, 在求解过程中系数a (x) 不再是自变量x的函数, 而是在整个区域内a (x) 为常数, 可以看出系数a (x) 与权函数的选取有关, 因此, 曲线拟合就退化成最小二乘法的曲线拟合。这样权函数形成的近似为线性拟合。

由图3所示:权函数为紧支常数, 只有在中心点x (28) 5.0权函数存在非零区域, 非零区域尺寸即为节点权函数支持域尺寸, 这样在任意一点x的支持域内均有两个节点对其有贡献 (n (28) m (28) 2) , 权函数为1, 其他节点的权函数为零, 权函数构造的近似就变成了分段线性插值, 类似有限元法的函数逼近。由图4所示:权函数为紧支且光滑的函数 (三次样条函数) , 计算点x的支持域内包含的节点数大于基函数的个数 (n (29) m) , 为使局部近似误差最小, 加权平方和J (x) 对a (x) 求偏导, 最终得到MLS的近似位移, 这样的求解过程用于EFG方法中, 即使是线性基函数, 近似值也是连续光滑的, 因为MLS近似继承了权函数的连续光滑性。因此, 本部分的分析, 选取三次样条函数。

4 曲线拟合

采用常数基函数]1[, 拟合结果如图5所示。

采用一次基函数[1x], 拟合结果如图6所示。

采用二次基函数[1xx2]拟合结果如图7所示。

由图7可以看出, MLS拟合函数及其一阶导数拟合的效果比较好, 二阶导数拟合呈现稍微的震荡。为了更加精确的分析拟合结构, 引入相对误差范数公式:

精确解公。式中:N为样本点个数, yN和yE分别为对应的数值解和

在移动最小二乘法中基函数的选取也有一定的规律, 从以上分析, 结合表1可以看出:同一个函数采用不同基函数, 对计算结果有直接的影响, 线性基的基函数的次数越高, 逼近效果越好且数值计算的相对误差也越小, 二次基函数能够得到良好的拟合效果。

5 曲面拟合

采用MLS法拟合以下函数:

设定区域在区域W上均匀布置8 (8个节点, 图8显示了精确值的曲面分布图, 利用MATLAB中Griddata函数插值得到的。为进行曲面拟合, 采用具有C0 (W) 连续的高斯型权函数 (b=30.) 和线性基函数[1xy], 在区域上均匀划分25 (25个规则的节点网格, 由MLS法公式, 求出各节点的近似函数值, 绘制网格曲面, 如图10所示, 从图中可以看出, 用移动最小二乘法拟合的曲面比较光滑, 并与Griddata插值形成的曲面相比, 如图9所示, 曲面光滑准确, 从而验证该方法的有效性。

为了衡量拟合的质量, 定义拟合的相对误差:

其中:zEXACT (x, y) 是精确解, zMLS (x, y) 是使用移动最小二乘法计算值, 本文得到的拟合误差为:L=2.1526, 这说明MLS拟合具有比较高的精度, 虽然MLS拟合的计算量稍微大一些, 和其精确性相比, 是可以接受的。

6 结语

移动最小二乘法虽然不具备传统有限元或边界元形函数所具有的插值特征, 即不满足Kronecker函数的属性, 给边界条件处理带来了困难, 但是它采用了较低的基函数通过选取适当的权函数, 获得具有较高连续性和相容性的形函数, 数值精度较高, 数据简单等特点, 是其他数值近似方法无法比拟的, 对以后的数值近似的研究具有重要的参考价值。

参考文献

[1]王能超.数值分析简明教程[M].北京:高等教育出版社, 2003.

[2]常青, 余湘娟.软土次固结变形特性的试验研究[J].南水北调与水利科技, 2008 (1) :148-150.

[3]张荣.利用VB实现水泵性能曲线的自动绘制[J].水利科技与经济, 2006 (1) :62-64.

最小二乘法曲线拟合 篇5

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

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

1 最小二乘曲线拟合

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

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

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

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

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

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

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

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

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

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

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

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

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

在MATLAB中输入以下程序:

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

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

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

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

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

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

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

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

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

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

例2假设已知

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

在MATLAB中输入以下命令:

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

运行结果为:

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

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

在MATLAB中输入以下命令:

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

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

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

3 结论

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

参考文献

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

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

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

最小二乘法曲线拟合 篇6

航空运输行业的快速发展, 为人们的出行提供了便利。但是随着空中交通流量的持续增长, 航空延误、空域拥堵等问题时有发生, 使得航空器冲突探测与解脱、进离场航班排序、基于轨迹运行等空管自动化与智能化方法成为空中交通管理领域的研究重点, 而如何快速准确的进行航空器飞行实时轨迹预测是实现上述方法的基础与保障。

目前飞行轨迹预测的方法主要有两种。第一种是基于混合估计理论实现的航迹预测, 利用交互式多模型 (IMM) 算法通过对横向、纵向以及垂直方向三种状态估计加权求和实现随机线性混杂系统的状态估计, 进而预测出航迹。第二种是基于航空器动力学及运动学模型, 利用各类机型的性能参数实现航迹的实时预测。如文献提出根据飞行阶段特点, 用基本飞行模型构建水平航迹、高度剖面和速度剖面, 根据航迹特征点的飞行状态信息拟合生成完整的4D航迹。

上述提出的两类方法用于飞行轨迹的预测存在以下两种问题:一是航空器的飞行模拟状态需要随实际情况的变化而变化, 而IMM算法中假设的模拟转换矩阵为固定值。二是在不考虑气象信息或不能准确判断航空器飞行意图的情况下, 很难保证航迹预测的准确度。

为了解决上述问题, 本文中根据已知的飞机航行轨迹, 判断当前所处阶段为直线航行或曲线航行, 估测后续飞行轨迹点, 进而拟合出行进轨迹的方法。最后, 再对因轨迹预测可能产生的偏差分情况分别进行纠正, 形成一条符合实时环境状况、飞行器飞行特征和飞行意图的平滑的行进轨迹。

方法原理

拟合原理

最小二乘法最早被称为“回归分析法”, 由著名统计学家道尔顿于1806年所创, 当时作为一种数学优化技术, 解决了日常生活中的计算问题。它通过最小化误差的平方和寻找数据的最佳函数匹配。在现在, 通过这种方法可以对未知数据进行预测, 并且极大化地减小误差, 求得与真实数据误差最小的预测数据。常用来解决曲线拟合、最佳逼近、物体运动轨迹等问题。其他一些优化问题也可通过最小二乘法来解决。

当存在大量实验数据时, 不能绝对要求拟合函数ϕ (x) 在数据点 (xi, yi) 处的误差, 即δi=ϕ (xi) -yi (i=, 12....) 为0, 但为了使曲线尽量逼近真实数据点的变化趋势, 将误差的平方和最小, 这就是最小二乘法的原理。

航迹的分段

就典型的航空器飞行航迹来说, 航空器一般都是以直线的飞行轨迹从一个航迹点飞行至另一个航迹点, 或者在空中进行了转弯等飞行动作, 以曲线的飞行轨迹从一个航迹点飞行至另一个航迹点, 这些航路点信息可以从领航计划报 (filed flight plan message, FPL) 中得知。即航迹可由一系列的直线飞行航迹段和曲线飞行航迹段组成。

直线飞行航迹段是飞行过程中的最简单的部分, 如图1所示, P1 (x1, y1, z1) 点和P2 (x2, y2, z2) 点分别为一条直线飞行航迹段的起点和终点, 在空间中两点之间的直线距离为

曲线飞行航迹段是由调整航路或航班避让以及进场盘旋等原因形成。通常, 转弯可分为三种模型:内切转弯, 末端转弯和约束转弯。

航路上最常使用的为内切转弯, 如图2所示。其中点p 1 (x 1, y 1) 为转弯初始点与p 2 (x 2, y 2) 转弯结束点。

故在实际航迹的预测中, 需判断得出一段航迹点中的直线段与曲线段, 将相应段的航迹点进行划分, 从而得到相应的直线段和曲线段拟合方程。具体航迹分段方法如下:

在典型的水平飞行航迹中, 若有A, B, C三个相邻的航迹点, 其中A (Xa, Ya, Za) , B (Xb, Yb, Zb) , C (Xc, Yc, Zc) , 则线段AB的直线方程为, 斜率线段BC的直线方程为, 斜率。故两直线间的夹角θ有, 则θ=arctan P。若此时θ<10°, 则A, B, C三点间用直线连接, 如图3所示。若此时θ>10°, 则A, B, C三点间用曲线连接, 如图4所示。

实际航迹预测过程中的分段过程如图5所示。若此时有a, b, c, d, e, f, g, h, 8个航迹点。第一步首先判断a, b, c三点, 经过计算可得到直线ab与直线bc间的夹角大于170°, 故判断a到c点之间以直线连接;第二步, 计算直线ac与直线cd之间的夹角, 经计算得出, 夹角小于170°, 故此时的处理过程为, 返回c点之前的点b, 此时bd段可判断为曲线;第三步, 计算直线cd与直线de之间的夹角, 经计算得出, 夹角小于170°, 此时ce段可判断为曲线, 又因为bd属于曲线段, 故此时be段属于曲线部分;第四步, 计算直线de与直线ef之间的夹角, 经计算得出, 夹角大于170°, 此时e点为曲线和直线的连接点, 故ef段属于直线, 采用相同的方法计算直线ef和直线fg之间的夹角, 直线fg, 和直线gh间的夹角, 得出eh段属于直线部分。此时如图5所示, 将整条航迹分为了三个部分, 即直线部分ab段, 曲线部分be段和直线部分eh段。

实时航迹矢量化方法

在非实时航迹的矢量化中, 可以通过对已知所有航迹点进行拟合, 从而得到整条非实时的航迹;但是在实时航迹的矢量化中, 仅通过预缓存的若干点进行最初阶段航迹拟合, 然后对下一时间点的实时航迹位置进行预测, 在得到真实的该时间点航迹位置后, 又采用相应偏差纠正方法对航迹进行平滑纠正。航迹点的预测主要分为直线航迹点的预测以及曲线航迹点的预测 (此处我们用椭圆方程进行曲线拟合) 两类。

直线航迹点的预测

方法1:由于所预测的航迹点都是以4s为时间间隔, 故在短时间内的运动皆近似为匀速运动处理。如图6所示。

其中A、B、C、D分别为时间间隔4秒的真实航迹点, E为预测的下一个航迹点。具体预测方法为:采用最小二乘法拟合得到直线方程, 其中以航迹点的纬度和经度分别定义坐标轴的横纵坐标;因该段过程我们采用匀速运动处理, 故 (其中Xi表示i点的横坐标, i为A、B、C、D、E) , 将Xe的值代入拟合得到的直线方程中, 得到E点的纵坐标Ye, 综上我们得到预测点E的坐标 (Xe, Ye) 。

曲线航迹点的预测

此处, 我们采用椭圆方程对曲线进行拟合。

(1) 焦点在坐标轴上

方法2:由于所预测的航迹点都是以4s为时间间隔, 故在短时间内的运动皆近似为匀速运动处理。如图7所示。

其中, A、B、C、D分别为时间间隔4s的真实航迹点, E为预测的下一个航迹点。又因为该椭圆的焦点分别位于X轴Y轴上, 故该椭圆的方程可采用参数方程:

其中a, b分别为椭圆的长短半轴, θi为i点的参数 (i为A、B、C、D、E点) 。原点O的坐标 (0, 0) , 设A点坐标 (Xa, Ya) , D点坐标 (Xd, Yd) , E点坐标 (Xe, Ye) 。若向量, 间的夹角为ϕ, 则, θe=θd+ϕ/4, 将参数θe代入参数方程, 即可得到E点的坐标。

(2) 椭圆焦点不在坐标轴上

方法3:由于所预测的航迹点都是以4s为时间间隔, 故在短时间内的运动皆近似为匀速运动处理。如图8所示。

其中, A、B、C、D分别为时间间隔4s的真实航迹点, E为预测的下一个航迹点, O’点坐标为O’ (m, n) 。其中, 椭圆的焦点不在坐标系XOY上。虚线所表示的为原始坐标系XOY, 实线所表示的是经旋转和平移后的坐标系X’O’Y’, 两坐标系横轴的夹角 (即旋转角) 为ν。故, 两坐标系之间的转化关系为

又因为该椭圆的焦点分别位于旋转平移后的X’轴Y’轴上, 故该椭圆的方程可采用参数方程

实时航迹纠正偏差方法

通过实验发现当飞机行驶方向发生偏转时, 由于预测具有延迟性, 预测点与真实点存在偏差, 如果立即按照新预测的航迹线进行飞行, 就会造成整体的航迹线在转弯处转弯角过大。为了解决上述问题, 从直线段到直线段的过渡、曲线段到曲线段的过渡、直线段到曲线段的过渡及直线段到曲线段的过渡以上四种情况对航迹进行纠偏, 使得经过纠偏处理后的数据拟合出的航迹更加平滑。

(1) 实时航迹直线段到直线段过渡偏差修正

如图9所示, 按照拟合完成的直线A’B’方程, 飞机由A’点飞至B’点, 将A’B’段的航迹预测看做理想状态, 故此处的实际航迹AB段与预测航迹A’B’重合。在B’点之后, 按照本文第三部分中所介绍的实时航迹点预测方法1, 预测得到4s后的航迹点C’, 但飞机经过4s后由B’点飞行至C’点, 此时雷达传来真实的航迹点C点。为了将模拟航迹平滑过渡至真实航迹, 采用以下方法:用C点以及该时刻向后倒推4s, 8s, 12s的航迹点, 采用最小二乘法拟合出直线方程“f1=k1x+b1”;采用相同的方法, 拟合出关于C’的直线方程“f2=k2x+b2”, 当︱k1-k2︱<0.1时, 令k=k1, ;当︱k1-k2︱>0.1时, 令。所预测的下一时刻实时航迹点D’坐标, 由此时所得到的直线D’B’的方程“f=kx+b”按照实时航迹点预测方法1得出。不断地迭代此过程, 直至预测航迹与实际航迹接近乃至重合。故整个预测模拟航迹为点A’B’C’D’E’F’的连线。

(2) 曲线段到曲线段的过渡偏差修正

(3) 曲线段到直线段的过渡偏差修正

如图11所示, 按照拟合完成的曲线AD方程, 飞机由A点飞至D点, 将AD段的航迹预测看做理想状态, 故此处的实际航迹AD段与预测航迹A’D’重合。在D’点之后, 按照本文第三部分中所介绍的实时航迹点预测方法2或3, 预测得到4s后的航迹点E’, 但飞机经过4s后由D’点飞行至E’点, 此时雷达传来真实的航迹点E点。为了将模拟航迹平滑过渡至真实航迹, 采用以下方法:用E点以及该时刻向后倒推4s, 8s, 12s的航迹点, 采用最小二乘法拟合出直线方程“f1=k1x+b1”;采用相同的方法, 拟合出关于E’的直线方程“f2=k2x+b2”, 当︱k1-k2︱<0.1时, 令k=k1, ;当︱k1-k2︱>0.1时, 令。所预测的下一时刻实时航迹点F’坐标, 由此时所得到的直线C’F’的方程f=kx+b按照实时航迹点预测方法1得出。不断地迭代此过程, 直至预测航迹与实际航迹接近乃至重合。故整个预测模拟航迹为点A’B’C’D’E’F’G’的连线。

(4) 直线段到曲线段的过渡偏差修正

如图12所示按照拟合完成的直线AC方程, 飞机由A点飞至C点, 将A’C’段的航迹预测看做理想状态, 故此处的实际航迹AC段与预测航迹A’C’重合。在C’点之后, 按照本文第三部分中所介绍的实时航迹点预测方法1, 预测得到4s后的航迹点D’, 但飞机经过4s后由C’点飞行至D’点, 此时雷达传来真实的航迹点D点。为了将模拟航迹平滑过渡至真实航迹, 采用以下方法:用E点以及该时刻向后倒推4s, 8s, 12s的航迹点, 采用最小二乘法拟合出曲线段C’E’的椭圆方程, 采用相同的方法拟合出曲线段CE的椭圆方程。当时, 令;当时, 令。预测的下一时刻实时航迹水平面坐标点E’由曲线段CE’ (两条曲线中所夹的曲线) 椭圆的方程得出。不断地迭代此过程, 直至预测航迹与实际航迹接近乃至重合。故整个预测模拟航迹为点A’B’C’D’E’F’的连线。

实验验证

本实验以首都机场为例, 采用首都机场提供的2013年6、7月的飞机飞行数据及机场地理位置信息, 随机抽取1000条离场航迹进行实验。首先, 依据实时航迹预测方法对数据进行处理预测出实时航迹, 在平面坐标系内进行效果对比。利用预测出的航迹与实际航迹之间数据的剩余方差, 判断预测航迹与实际航迹之间的误差, 从而验证方法的有效性。

通过实验前后对比图可看出, 该实时航迹预测方法预测效果较好, 可满足实时航迹预测的要求。

航迹的预测精度即航迹的预测位置与实际位置之间的差值。通过计算该预测航迹与实际航迹的稳定程度来近似估算预测方法的精度具体来说主要步骤如下:

1) 对预测航迹或实际航迹S, 假设按时间顺序该航迹由点P1, ..., Pm组成, 对这些点的纬度x, 经度y值分别做差分得到两组序列X={X1, ..., Xm-1}, Y={Y1, ..., Ym-1}, 其中X1=Pi+1 (x) -Pi (x) , Yi=Pi+1 (y) -Pi (y) , 通过序列X, Y把目标在x, y轴上的运动自动分解成一系列的匀速运动和匀加速运动区间U1, ..., Uk;

2) 对每个区间Ui根据目标在该区间x, y轴上的运动模型用对应的线性或二次函数对该区间的运动轨迹进行回归分析, 得到目标在该区间的运行轨迹函数x (t) 和y (t) ;

4) 对航迹S上各区间的剩余方差取均值得到各航迹的剩余方差S2 (x) 和S2 (y) ;

5) 对计算得到的近1000条真实航迹的剩余方差取平均值。

通过实验结果得到, 四种情况下实际航迹与预测航迹的最大剩余方差平均值与最小剩余方差平均值的差值为0.03211, 该数量级的误差, 在航迹矢量化处理后显示在雷达, 将会是2至6个像素的误差, 故在允许的误差范围之内。故说明该实时航迹矢量化方法是一种正确可靠, 有效的实时航迹矢量化方法。

结语

最小二乘法曲线拟合 篇7

强震后余震是大地震发生后的常见灾害之一,震后余震对救灾活动、难民安置以及灾后重建等工作有着严重危害。所以,对震后余震的分析及预测就显得尤为重要。

文献[1]在研究了我国20世纪30年代到70年代接近50年内发生的28次地震的强余震序列后,发现了“等待时间”与余震“发生时刻”的对数线性关系,并对其作了初步理论解释。文献[2]分析了广州邻近区域的几次典型的余震序列,发现余震“发生时刻”的对数与“等待时间”对数具有较好线性关系,进一步印证了文献[1]的理论。文献[3]系统地分析了1966~1996年间70余次Ms5.5以上的地震的余震序列,发现强余震“等待时间”模型对6级以上主余震型序列预测成功率较高。文献[4]对多次大震和中强震进行比较研究,发现“等待时间法”对大震的强余震预测效果较好。文献[5]利用混沌理论对汶川地震余震“等待时间”序列进行研究,发现“等待时间”序列存在混沌动力学特征。文献[6]以“等待时间”作为权值弥补双对数变化对余震预测误差的影响,发现经过加权改正后,模型的预测精度有所提升。然而,现有的利用“等待时间法”进行余震预测的研究都是基于文献[1]的理论,仅考虑了“等待时间”的随机误差,运用最小二乘法(LS)或加权最小二乘法对文献[1]的模型进行求解,未能顾及“等待时间”模型中由“发震时刻”构成的系数矩阵的随机误差。

本文根据“等待时间”与余震“发生时刻”的对数线性关系,利用总体最小二乘方法(TLS),同时考虑“等待时间”观测向量的随机误差和以“发震时刻”构成的系数矩阵的随机误差来进行拟合;并将该方法应用在尼泊尔地震的强余震序列预测和拟合中,得到了有益的结果。

1 利用总体最小二乘法拟合的“等待时间”强余震预测模型

1.1“等待时间”强余震预测模型

连续发生的两个强余震之间的时间间隔为后一个强余震的“等待时间”。即:

文献[1]发现“等待时间”与“发震时刻”存在以下对数关系:

在式[1]及式(2)中i为第i次强余震的序号,ti为第i次强余震“发生时刻”(距主震时间),Δti则表示第i次强余震的“等待时间”,以上时间单位均为小时。β0,β1为拟合参数。式(2)中的β0,β1值可通过总体最小二乘方法拟合已知的强余震序列求出。根据得到的β0,β1值可以得到拟合直线,再通过最近一次余震的“发震时刻”对数即可预测出下次余震的时间对数,进而得到下次余震的“发震时刻”。

顾及观测向量和系数矩阵同时存在的随机误差的情况,建立EIV模型如下:

式中y为观测向量值,ey为观测值误差:

A为系数矩阵,EA为系数矩阵误差:

ξ为未知参数:

总体最小二乘的目标函数:

式(8)中vec为矩阵拉直符号。

利用总体最小二乘方法求解的迭代算法可表示如下[7]:

[1]以最小二乘方法得到的为起始点ξ^[1]:

(2)计算出残差:

(3)将式(9)、(10)得到的结果代入式(11)求出:

(4)重复步骤(2)、(3),当:

时终止迭代,本文中取ε=10-10。

精度评定:

1.2 作穷举线求交点

根据总体最小二乘方法拟合得出的β0,β1作拟合直线,并将求出的直线同穷举线作交点,穷举线绘制是依照以下公式进行的:

穷举线与拟合线的交点即为预测的第N+1个强余震,交点的横坐标即为预测下次强余震“发生时刻”的对数值,经过简单变换即可以得到预测时间[6]。

2 以尼泊尔地震为例的强余震序列拟合和预测

2015年4月25日在尼泊尔博拉克(北纬28.2°,东经84.7°)发生的Mw7.8强震给当地人民造成了巨大的生命和财产损失,接踵而来的余震也对救援和安置工作造成了很大的困难。本文分别应用最小二乘与总体最小二乘方法对尼泊尔地震强余震序列进行拟合及预测,并对结果进行讨论分析。

2.1 数据选择

对于余震目录的筛选,在利用“等待时间”序列进行余震预测时,一般采用主震震级减2~3级为余震目录的震级下限,并把时间间隔小于前一强余震“等待时间”的1/20的余震与前一次余震当作同一次能量释放,可以把它从余震目录中剔除[8]。

本例采用的余震目录是由美国地质勘探局(USGS)在其官网提供的尼泊尔Mw7.8级地震5级以上余震序列,并剔除了时间间隔小于前一强余震“等待时间”的1/20余震序列,时间截至2015年4月27日12:35:53(UTC时)。余震目录见表1。

2.2 震例分析

为了验证总体最小二乘方法的拟合效果,本文首先应用最小二乘法和总体最小二乘法对表1中的前17个样本进行拟合,利用总体最小二乘法得到β0=-0.638,β1=1.087,σ02=0.054。利用最小二乘法得到β'0=-0.612,β'1=0.980,σ'02=0.132。以两种方法求得的单位权中误差σ02和σ'02作为拟合线残差。通过求得的参数即可绘出拟合直线,见图1。

根据式(14)作穷举线,穷举线与拟合直线的交点即为预测的第N+1个强余震。预测残差是下一次强余震的“发生时刻”对数值的预测值和实际值之差。分别用eLS和eTLS来表示最小二乘与总体最小二乘的预测残差。结果见图2。

为了体现样本个数对拟合和预测结果的影响,本文进行了如下实验:采用前n个样本预测第n+1次余震,分别利用最小二乘和总体最小二乘方法进行直线拟合并预测下一次余震,样本数从11个逐次增加到19个,进行9次实验,结果见表2。

通过表2及图3可以看出,在该震例中,总体最小二乘对于余震“等待时间”和“发震时刻”对数模型的拟合结果要优于采用最小二乘拟合方法,预测精度在总体上也比直接采用最小二乘有较好的结果。预测残差的绝对平均值由0.082降为0.071,降幅约为13%。拟合残差的绝对平均值由0.071降为0.059,降幅约为17%。图3中随着样本数的增加,两种方法的预测误差都向0值靠拢,而总体最小二乘的结果在整体上更加接近0值。

3 其它地震实例

为了进一步验证总体最小二乘方法应用于余震“等待时间”和发震时刻对数模型的拟合和预测效果,本文另选了几组国内外典型强震的强余震序列样本进行最小二乘方法与总体最小二乘方法对比实验,结果见表3及图4。

由表3和图4可以看到,利用总体最小二乘方法得到的结果比直接运用最小二乘方法得到的结果在拟合精度上有明显的优势,在预测精度方面也有着不同程度的提高。总体上,拟合残差绝对均值从0.133降低为0.053,降低幅度约为60%,预测残差绝对均值从0.170降低为0.148,降低幅度约为13%。其中日本311地震精度提高幅度最大,拟合残差绝对均值由0.301降低为0.078,预测残差绝对均值由0.101降低为0.066,降低幅度分别为74%和35%。

4 总结

本文提出了同时顾及“等待时间”观测向量的随机误差和以“发震时刻”构成的系数矩阵的随机误差,采用总体最小二乘方法来进行强余震“等待时间”序列的拟合,并将其应用到2015年尼泊尔地震的强余震预测中。再经过对10个典型震例的实验分析后,结果表明,总体最小二乘方法可以用于“等待时间”余震模型的拟合以及余震预测;在拟合方面,因为总体最小二乘方法顾及了系数矩阵的误差,所以较最小二乘而言有较好的提升;在余震预测精度方面也有一定的提高。但是由于本文所涉及的余震震例较少,可能未能全面反映总体最小二乘在该方面的优劣性,利用该方式进行实际的余震预测有待进一步研究。

摘要：本文针对“等待时间法”余震预测模型在拟合过程中未考虑系数矩阵的随机误差的问题,利用总体最小二乘方法进行该模型的拟合,并将该方法应用于尼泊尔地震及其它震例的强余震预测。结果表明,运用总体最小二乘方法进行的“等待时间”余震序列拟合和预测所得结果精度比传统的方法在总体上有一定的提升。