递推最小二乘法

2024-06-04

递推最小二乘法（共8篇）

递推最小二乘法篇1

带观测噪声的自回归 (AR) 模型广泛应用于通讯、信号处理、语音增强、目标跟踪等领域。周知, 当自回归 (AR) 模型的输入噪声是白噪声时, 用普通递推最小二乘法可得到AR参数的一致估计;但当输入噪声为有色噪声时, 用普通递推最小二乘法将引出有偏估计[1]。为此, 在假设观测白噪声方差已知的情况下, 文献[2,3]提出了带白色观测噪声的AR信号参数估计的偏差补偿递推最小二乘法。此后文献[4,5]又在假设有色观测噪声相关函数已知的基础上, 提出了带有色观测噪声的AR信号参数估计偏差补偿递推最小二乘法, 并用随机压缩映像原理严格证明了算法的收敛性, 推广了文献[2,3]的结果。后来文献[6,7,8]又在观测噪声方差未知的条件下给出了偏差补偿最小二乘算法, 但其中AR参数估计与观测噪声方差估计是互耦的, 这给收敛性证明带来了很大困难, 因而均没有给出严格的参数估计收敛性证明。以上文献都是针对单变量系统辨识的, 对于多变量系统, 在新近的文献[9]中也只是研究了多输入单输出 (MISO) 系统, 且没有给出参数收敛性证明。

所以本文针对多变量AR信号, 在观测白噪声方差阵假设已知的条件下, 提出了AR信号参数估计的多变量偏差补偿递推最小二乘法, 推广了单变量AR信号参数估计的偏差补偿递推最小二乘法[2,3]。并用动态误差系统分析的方法[10]严格证明了该算法的收敛性, 即偏差补偿递推最小二乘参数估计以概率1收敛于相应真实值。

1 多变量偏差补偿递推最小二乘法

考虑带白色观测噪声的多变量平稳AR (n) 信号

其中y (t) ∈Rm为真实AR信号, z (t) ∈Rm是对y (t) 的观测信号, 且ε (t) ∈Rm和v (t) ∈Rm是零均值、方差阵各为Qε和R的相互独立白噪声。已知R, 但Ai和Qε未知, 问题是基于观测 (z (1) , …, z (t) ) 用偏差补偿原理求它们的最小二乘估值。

将式 (2) 代入式 (1) 引出

式 (3) 有LS结构

于是由多维RLS算法[11]可得有偏的Θ (t) 的估值 $\hat{Θ}_{b} (t)$ 为

并有非递推LS估值的正规方程为

将式 (4) 代入式 (10) 并两边除以t整理得

由y (t) 的平稳性假设及遍历性[12], 当t→∞时以概率1 (w.p.1) 有

由式 (2) , 式 (5) , 式 (7) 有

由平稳过程的遍历性[11], 当t→∞时以概率1有

于是当t→∞时有

定义

所以

由式 (12) , 式 (13) , 式 (18) 得出渐近偏差为

由P (t) 的定义式 (11) , 当t→∞时以概率1有

当t充分大时, 可近似用tP (t) 代替M-1, 式 (19) 右边的Θ可用估值 $\hat{Θ} (t - 1)$ 近似代替, 于是由式 (19) 可得多变量偏差补偿递推最小二乘估值为

上述推导可概括为如下定理。

定理1 带白色观测噪声的平稳AR (n) 模型式 (1) 、式 (2) 的多变量偏差补偿递推最小二乘法AR参数估计为

其中取初值 $\hat{Θ} (0) = 0$ , $\hat{Θ}_{b} (0) = 0$ , P (0) =αImn, α>0。

下面估计噪声方差阵Qε。计算式 (7) 两边随机过程的方差得

方差阵Qe的估值 $\hat{Q}_{e}$ 可由递推公式计算得到

将 $\hat{Q}_{e} (t)$ 和 $\hat{A}_{i} (t)$ 代入式 (25) 即可得到Qε在t时刻的估值

2 收敛性

用动态误差系统分析[10]的方法严格证明偏差补偿递推最小二乘法的收敛性, 首先给出以下几个引理。

引理1[13] 设A, B均为n阶对称矩阵, 且A≥0, B>0, 则B>A的充分必要条件是谱半径ρ (AB-1) <1。

引理2[10] 考虑一个时变动态误差系统

式 (29) 中t≥0, 输出 (动态误差) δ (t) ∈Rp×q, F (t) ∈Rp×p, 输入u (t) ∈Rp×q。假设F (t) →F (t→∞) , F是一个稳定矩阵, 且假设u (t) 是有界的, 则有δ (t) 是有界的。

引理3[10] 考虑一个稳定的动态误差系统

式 (30) 中t≥0, 输出 (动态误差) δ (t) ∈Rp×q, 输入u (t) ∈Rp×q。假设F∈Rp×p是一个稳定矩阵, 且设u (t) →0 (t→∞) , 则有δ (t) →0 (t→∞) 。

引理2和引理3的证明可用完全类似于文献[10]的方法推导, 详细推导从略。

定理2 带白色观测噪声的系统的多变量偏差补偿递推最小二乘法参数估值以概率1收敛于真实值, 即

证明由式 (20) 得

且

由式 (19) , 式 (24) , 式 (32) 得

记

则式 (34) 可表示为

由式 (2) 引出

对式 (39) 取方差运算引出

定义随机过程y (t) 的相关函数为R (k) =E[y (t) y (t-k) Τ], 则有R (-k) =RΤ (k) , 且

因方差阵是对称非负定阵, 故对任意mn×1实向量有

为了证明Q是正定的, 只需证明

的充要条件是x=0。即

充分性是显然的。现证明必要性, 首先构造函数

且由y (t) 的平稳性知y (t) 可表示为均方收敛的级数[11]

式 (46) 中规定G0=Im, 系数阵Gi可递推计算为Gi=-A1Gi-1-…-AnGi-n, i=1, 2, …, Gi=0 (i<0) 。于是有

式 (47) 中*代表相同于左边中标号内的项, 常数Δ≥0, 假定Qε>0, 由式 (47) 左边为零可推出x1=0。重复上述步骤可推出x2=…=xn=0, 即证明了必要性。故有

所以由式 (40) 有

因为矩阵Rd和M均为mn×mn阶对称阵, 且有Rd≥0, M>0, 由引理1引出ρ (RdM-1) <1, 因为矩阵RdM-1和M-1Rd具有相同的特征值, 于是有ρ (M-1Rd) <1, 即F的谱半径ρ小于1, 即F为稳定矩阵。

由式 (22) 得

注意到

由式 (19) 知 $\lim_{t \to \infty} \hat{Θ}_{b} (t)$ 极限存在, 所以 $\hat{Θ}_{b}^{Τ} (t)$ 有界。由引理2引出 $\hat{Θ}^{Τ} (t)$ 有界。又因为式 (33) 和式 (36) 及 $\hat{Θ}_{b} (t) \to \lim_{t \to \infty} \hat{Θ}_{b} (t) (t \to \infty)$ , 所以

故应用引理3到式 (38) 有

根据式 (35) 即有

假设白噪声ε (t) 和v (t) 以概率1有界, 则由y (t) 的平稳性假设引出y (t) 和z (t) 以概率1有界。从而由式 (4) 和式 (26) 引出

注意式 (27) 的非递推表达式为

而由遍历性[12]有

故也有[10]

由式 (58) 及式 (25) , 式 (28) 可得

证毕。

3 仿真例子

考虑带白色观测噪声的AR信号

其中y (t) 为真实AR信号, z (t) 是对y (t) 的观测信号, 且ε (t) 和v (t) 是零均值、方差阵各为Qε和R的相互独立白噪声。已知R, 但A1, A2和Qε未知。仿真中取模型参数为和, 噪声方差阵为和, 问题是基于观测 (z (1) , …, z (t) ) 如何辨识A (q-1) 和Qε。

仿真结果如图1至图3所示, 其中直线表示真实值, 实曲线表示估值。由图1看到普通递推最小二乘法参数估计是有偏的, 由图2看到偏差补偿递推最小二乘法参数估计是一致的, 且由图3看到噪声方差估计也是一致的。在t=8 000处普通递推最小二乘算法模型参数估值为

偏差补偿递推最小二乘算法模型参数估值为

噪声方差阵估值为

4 结论

本文推广单变量偏差补偿递推最小二乘算法到多变量情形, 提出了带白色观测噪声的多变量AR信号参数估计的偏差补偿递推最小二乘算法。利用动态误差系统分析方法严格证明了用该方法得到的参数估计和噪声方差估值的强一致性。

摘要：对于带白色观测噪声的多变量自回归 (AR) 信号, 提出了未知模型参数和噪声方差估计的偏差补偿递推最小二乘算法, 用动态误差系统分析方法严格证明了所得到的模型参数和噪声方差估值是强一致的, 即它们以概率1收敛于相应真实值。一个仿真例子说明了其有效性。

关键词：多变量 (AR) 信号,参数估计,偏差补偿递推最小二乘法,收敛性,强一致性,动态误差系统分析方法

递推最小二乘法篇2

【关键词】BP神经网络最小二乘法盾构故障

一、故障预测常用方法种类及特点

常用的故障预测种类有以下几种：

（一）曲线拟合法：本方法简单实用，好理解，较多的应用在监控机械装置的系统，缺点是预测误差结果比较大。

（二）卡尔曼滤波法：计算量小、预测精度高，主要用于线性系统，在非线性系统中需要进行扩展，但其模型的不确定性很差。

（三）灰色预测：该方法需将两个部分进行置换，即“随机过程”转换为“灰色过程”，把“随机量”处理成“灰色量”，然后在理论系统的模型GM（1，1）中进行处理。

（四）神经网络：具有学习记忆功能，能很好解决非线性系统问题，不对预测模型做限制，同时还能把以前的历史数据完整映射到未来的数据库中，可广泛用在故障检测中。实际使用广泛的神经网络预测有：BP神经网络和自组织特征映射网络。

二、最小二乘法与BP神经网络相结合的预测方法

基于本文采用最小二乘法与BP神经网络相结合的方法来预测故障，因此需要对所选取的数据样本集进行直线最小二乘法的拟合，然后用BP神经网络算法进行计算，最后将这些值作为样本集，输入到BP神经网络模型来进行训练预测。

（一）BP神经网络模型的构建

在构建神经网络之前，需将数据进行归一化处理，将归一化的数据样本构建成BP神经网络模型。在本论文的神经网络模型中，训练函数，学习函数，传读函数，性能函数，训练次数，训练误差，学习率分别为：traingdx，learngdm，logsig，mse，1000次，0.0001，0.08.输入层，输出层及隐层分别为：6，1，5。

（二）最小二乘法的直线拟合过程

根据最小二乘法的理论[4]及BP神经网络的模型可知，要预测故障，既是对出现故障的时间序列的预测，而预测时间序列X={Xt︱t=0，1，2…m}，实际上就是利用已知的{xk，xk-1…xk-n+1}来求xk+1。

将得到的数据样本再通过BP神经网络的算法对样本计算，得到其对应的各个参变量值，即给出Xk，得到期望值N（Xk），其中k=n-1，n…m，然后把这些值作为样本集，用在BP神经网络中来训练整个网络。

三、盾构机PU电流序列预测应用

（一）单一的BP神经网络的预测结果

单一的BP神经网络采用以1小时为间隔所提取得 100个PU电流数据作为样本集对进行训练，采用输入连续的6个PU电流值，紧接着输出相连续的一个PU电流值，即输入与输出层节点数比例为6：1，隐层节点数为5，将网络进行收敛，当达到预定精度以后，预测其结果如图1中曲线。

（二）最小二乘法及BP神经网络相结合的预测结果

用上面的样本集数据，对最小二乘法与BP神经网络的合成的样本进行训练，并进行网络收敛，当达到预定精度，预测其结果如图2中曲线所示。

（三）两者结果比较

图1和图2中的两条预测线分别代表PU实际电流变化曲线（实线）以及用神经网络预测的电流变化曲线（虚线）。可以看到，在图2中，这两种曲线的变化趋势比较接近，仅有极少的时间序列点有较大偏差，而图1中的两条曲线则在多个时间序列段出现较大的偏差，两条曲线走势明显偏离，所以，可以认为，采用合成的BP神经网络算法预测到的PU电流时间序列精度较高，比采用单一的BP神经网络更能客观反映盾构机的实际电流变化情况。虽然有个边点出现误差，也是未考虑机组负荷的影响。因此，整体上说合成的BP神经网络算法对盾构机的PU电流还是具有很高的预测能力。

图1 图2

四、总结

上述内容主要采用最小二乘法与BP神经网络结合的方法对盾构机的故障进行预测，将二者结合起来，可以充分利用二者的优点，优势互补。笔者将其应用到了实际的工作中，起到了良好的预测效果。

参考文献：

[1]蒋瑜，杨雪，阮启明.机械设备故障规律及运行趋势预测方法综述[J].机电一体化，2001（3）：14-17.

[2]陈敏泽，周东华.动态系统的故障预报技术[J].控制理论与应用，2003，20（6）：820-822.

[3]苏春华，罗雷，海军，梅检民，肖云魁.基于BP神经网络的汽车发动机寿命预测[J].军事交通学院学报，2009，11（4）：50-51.

递推最小二乘法篇3

桥梁线形预测与调整是桥梁施工监控计算和分析的主要内容之一。在桥梁施工监控中, 对于设计参数误差的调整就是通过量测施工过程中实际结构的行为, 分析结构的实际状态和理想状态的偏差, 用误差分析理论来确定或识别引起这种偏差的主要设计参数误差, 来达到控制桥梁结构的实际状态与理想状态的偏差, 使结构的成桥状态与设计相一致。因此, 根据当前施工阶段结构的实际状态进行正装计算至成桥状态, 预告今后施工可能出现的应力和变形状态, 形成施工监控的2大任务:即结构的前期预报和后期调整。目前, 桥梁线形预测与调整主要是立模标高的预测与调整, 确定立模标高的理论和方法主要有Kalman滤波法、灰色理论法和最小二乘法等[1]。

1795年, K.F.Gauss提出了最小二乘法并应用于天文计算的实践中, 指出对于未知的但要求估计的参数的最适宜的值是最可能的值, 定义“未知量的最可能值是这样的一个值, 它使得实践值与计算值的差的平方乘以测量精度后所求得的和最小”。后来, 在控制系统的参数估计领域内也发现和采用了这种方法。以前, 在稳态系统数学模型的回归分析方面用得比较成熟和广泛。在20世纪60年代, 瑞典学者K.J.Astrom把这个方法用于动态系统的辨识中, 并取得了许多成果。最小二乘法在我国桥梁施工监控中的应用始于20世纪80年代后期, 并取得了较好的成果[2], 计算公式和方法详见文献[2]。

2 工程概况

某桥位于江苏省宿迁市, 为跨京杭运河特大桥, 主桥总长300 m, 为80 m+140 m+80 m三跨预应力混凝土变截面连续箱梁桥。桥梁设计荷载为:汽车-超20级, 挂车-120, 设计地震裂度为8级[1]。

主桥箱梁采用直腹板形式的双箱单室结构, 两箱之间由60 cm桥面板湿接缝连接。单箱顶板宽11.5 m, 底板宽度为6 m。边、中跨合龙段和边跨现浇段箱梁中心高度为3 m, 主墩墩顶箱梁中心高度为8 m, 相应位置处的底板厚度分别为32 cm和80 cm。单“T”箱梁梁高和底板厚度均按1.5次抛物线变化。箱梁采用50#混凝土。箱梁在0#梁段和中跨跨中分别设置了横隔板, 以增强桥梁的横向整体性。箱梁悬臂浇注节段划分长度共分为3 m和4 m 2种, 共19个悬臂节段, 83个施工节段。纵向预应力分为顶板束、腹板下弯束、边跨底板束、中跨底板束、顶板合龙束, 在空间上即有平弯又有竖弯。锚具采用OW15-12型锚具, 预应力孔道采用金属波纹管;横向预应力钢束主要布置于0#节段;竖向预应力布置间距为0.75~1.00 m, 采用冷拉Ⅳ级钢筋, 锚具采用轧丝锚。

在节段浇注循环过程中, 挂篮前移、混凝土浇注、预应力筋张拉等工况, 悬浇各节段均产生竖向线位移与截面内转角位移, 悬臂端竖向与转角位移与悬臂各个节段的材料性能与应力水平有关;由于混凝土材料的收缩与徐变, 随着时间的变化悬臂端部竖向与转角位移均会产生变化;梁体由于温度差, 热胀冷缩变形不均, 易导致梁端较大的竖向、转角及横向位移, 如环境温度变化、太阳辐射, 梁体内外、顶底板、上下游均会产生温度差, 从而产生位移。悬臂浇注施工挂篮就位之前, 需要分析与判断已浇梁体节段所处位置 (高程) 与理想状态的偏离量, 然后确定下一段的立模标高, 以期成桥后, 控制点高程与设计高程偏差在允许范围之内。

3 节段立模标高预测

某特大桥主桥上部箱梁共19个悬臂节段, 根据施工工序、力学性状和监控的难度, 一般分3个阶段进行有侧重的监控:1#~7#节段施工为施工监控前期;8#~16#节段施工为施工监控中期;17#~19#节段施工为施工监控后期。在施工监控前期, 由于主梁悬臂矩还比较短, 其刚度相对较大, 节段自重、预应力张拉和其它施工荷载对主梁挠度、应力影响较小;而在施工后期, 随着悬臂变长, 主梁的相对刚度逐渐变小, 各种施工荷载的力学效应也就趋于明显, 尤其是对立模标高的调整比重加大, 所以在施工监控中期最适合对箱梁的力学特性进行反演评估, 并进一步预测成桥线形、调整后期施工立模标高。

以施工监控中期的10#节段、16#节段施工阶段做为算例来说明递推最小二乘法在连续梁桥挂篮施工监控中的应用, 并通过施工监控后期19#节段施工实际成桥线形的误差分析来评价递推最小二乘法预测、调整标高的效果。

首先进行桥梁施工力学的前进分析, 评估箱梁力学参数对监控目标的敏感性, 其次根据各参数的敏感性分为主要设计参数和次要设计参数, 并结合实际桥梁施工过程中各参数发生偏差大小的可能性, 以选取合适最小二乘法辨识的力学参数[3]。参考桥梁监控相关文献和本桥施工的实际情况, 选取混凝土弹性模量、混凝土徐变系数、有效预应力系数作为重点辨识参数[2,4], 由此建立的参数调整向量为

式中:θ (1) 为混凝土弹性模量的调整系数;θ (2) 为混凝土徐变的调整系数;θ (3) 为有效预应力的调整系数。

计算三者偏差±10%对预拱度的影响, 发现正负偏差对预拱度的影响是对称的, 单位偏差 (1%) 对监控目标 (预拱度) 的敏感程度不足3 mm, 因而一般假定辨识参数偏差10%范围内, 参数误差向量测向量误差的线性变换矩阵是线性的。

在10#节段混凝土浇注后, 以4#节段、6#节段、8#节段、9#节段、10#节段监测点信息为量测向量, 见式 (2) 和表1:

则实测误差向量见式 (3) :

在10#节段预应力张拉阶段后, 以4#节段、6#节段、8#节段、9#节段、10#节段监测点信息为量测向量, 见表2和式 (4) 。

则实测误差向量见式 (4)

在16#节段混凝土浇注后, 以8#节段、10#节段、14#节段、15#节段、16#节段监测点信息为量测向量, 见表3和式 (5) 。

则实测误差向量见式 (6)

在16#节段预应力张拉阶段后, 以8#节段、10#节段、14#节段、15#节段、16#节段监测点信息为量测向量见表4。

则实测误差向量

初步分析4个实测误差向量, 在混凝土浇注后, 节段自重的增大导致梁体发生的变形为主要变形, Y1、Y3的实测值比计算值偏大, 可能是实际混凝土综合弹性模量较计算设计参数值大造成的。在预应力张拉后, 预应力荷载的施加导致梁体发生的变形为主要变形, 并伴随着徐变的发生, 由于Y2、Y4的实测值比计算值偏小, 可能是箱梁混凝土的实际有效预应力较计算设计参数值小以及徐变的不定性造成的。

在实际监控中, 由于新的监测信息随着工程进展而不断反馈, 所以一般的监控系统都采用最小二乘法的递推算法, 不必重新对大量旧数据算一遍, 而是设法在原先预测的基础上加入新信息以修正, 从而得到新的参数预计值。

Y1对应的线性变换矩阵φ1

Y2对应的线性变换矩阵φ2

Y3对应的线性变换矩阵φ3

Y4对应的线性变换矩阵φ4

对于加权系数矩阵, 其对角矩阵值的确定牵涉很多因素, 其中包括结构本身的各类客观因素、设计参数的可能变动幅度, 各监测点、监测工况信息可靠度的评估, 施工工况后线形的光顺程度以及工程技术人员的工程经验等主观因素。因此难以准确推导出加权系数矩阵。

考虑到节段推进, 箱梁相对刚度降低, 施工的挠度效应较大, 更易于被监测设备所观测, 即观测误差占总的监测量的比重下降, 因而在同等条件下, 应加重后期监测信息的权重。根据结构模型计算可知, 实际发生的挠度线很可能是二阶连续可导的光滑平顺曲线, 一般可加重悬臂距长的监测量的权重, 计算结果见表5。

在预测终期 (17#节段、18#节段、19#节段) 施工时3个参数的调整值和取值见表6。

根据预测的参数结果, 进一步指导调整19#节段施工立模标高, 采用新参数值预测挠度值和实测数据比较见表7和表8。从表可见, 经过模型参数的调整之后, 19#节段施工的阶段误差不足5 mm, 所占比例多数低于10%, 满足了悬臂施工之后体系转换的合龙精度要求。

4 线形控制精度

采用递推最小二乘法对某桥施工监控工作始于2004-11, 至2005-05下旬两边跨成功合龙, 2005-06上旬中跨顺利实现合龙。大桥边、中跨合龙前后线形良好, 悬臂端相对高差及高程偏位见表9, 均满足规范要求。

5 结论和建议

(1) 桥梁线形预测与调整是桥梁施工监控的主要工作, 递推最小二乘法数学基础坚实, 物理意义明确, 结合某桥施工监控中的应用表明, 该法能够满足精度要求, 可用于施工监控中的线形调整;

(2) 影响桥梁线形的因素较多, 本文采用的递推最小二乘法以混凝土弹性模量、混凝土徐变系数、有效预应力系数作为重点辨识参数, 工程应用表明, 参数选择合理有效, 但对于不同桥型应根据实际情况确定辨别参数;

(3) 加权系数矩阵的选取是递推最小二乘法中的关键因素, 应综合考虑桥梁结构、设计参数, 各监测点、监测工况信息可靠度的评估, 施工工况后线形的光顺程度等, 合理调整该矩阵以提高预测精度。

参考文献

[1]唐云清.大跨度连续箱梁桥施工监控技术研究[D].南京:河海大学, 2007, 4.

[2]向中富.桥梁施工监控技术[M].北京:人民交通出版社, 2001.

[3]李国平, 刘键.大跨连续梁桥线形最优施工监控的理论方法[J].华东公路, 1992, 15 (2) :66-69.

递推最小二乘法篇4

最小二乘算法是系统辨识中用得最广泛的估计方法之一。标准的递推最小二乘算法是通过极小化关于输入输出数据的一个二次准则函数,即极小化估计残差平方和而得到的算法[1]。本文基于卡尔曼滤波原理[2],从一个新的角度来推导递推最小二乘算法:定义一个参数估汁误差协方差矩阵,通过极小化该协方差矩阵而推导出递推最小二乘算法。

1最小二乘辨识算法

考虑下列单输入单输出(SISO)系统的参数估计问题,

A(z)y(t)=B(z)u(t)+v(t)(1)

(1)式中{u(t)}和{y(t)}分别是模型的输入和输出序列,{v(t)}为零均值、方差为σ2的不相关随机白噪声序列A(z)和B(z)均为单位后移算z-1的多项式[z-1y(t)=y(t-1)],且

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

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

定义参数向量θ和信息向量φ(t)分别为

θ:=[a1,a2,…,ana,b1,b2,…,bnb]T∈Rn,

φ(t):=[-y(t-1),-y(t-2),…,-y(t-na),u(t-1),u(t-2),…,u(t-nb)]T∈Rn。

那么式(1)写可成下列最小二乘辨识模型,

y(t)=φT(t)θ+v(t)(2)

本文的目标是,利用系统的输入输出数据{u(t),y(t)}或{y(t),φ(t)},基于卡尔曼滤波原理,极小化参数估计误差协方差阵,而推导辨识参数向量θ的递推最小二乘算法。

很多文献给出了计算θ估计白 $\hat{θ}$ (t)的递推最小二乘算法[1,3]:

$\hat{θ} (t) = \hat{θ} (t - 1) + L (t) [y (t) - φ^{Τ} (t) \hat{θ} (t - 1)]$ (3)

L(t)=P(t-1)φ(t)[1+φT(t)P(t-1)φ(t)]-1(4)

P(t)=[I-L(t)φT(t)]P(t-1)(5)

其中 $\hat{θ}$ (t)为t时刻的估计,L(t)称为增益向量,P(t)称为协方差矩阵。

2算法推导

2.1单输入单输出系统的递推最小二乘算法

仍考虑上述单输入单输出系统式(1)或式(2)。参照式(3),设t寸刻的参数估计 $\hat{θ}$ (t)足由前一时刻参数估计 $\hat{θ} (t - 1)$ 加上一个补偿项得到,且具有下列形,

$\hat{θ} (t) = \hat{θ} (t - 1) + L (t) [y (t) - φ^{Τ} (t) \hat{θ} (t - 1)]$ (6)

那么问题就是如何找到增益向量L(t)∈Rn。方法如下,定义参数估计误差为

$\hat{θ}$ (t): $= \hat{θ} (t) - θ$ (7)

将式(6)代入式(7),并利用式(2)展开可得

$\hat{θ} (t) = [Ι - L (t) φ^{Τ} (t)] \hat{θ} (t - 1) + L (t) v (t)$ (8)

(8)式中I是一个适当维数的单位阵.目标是确定一个最优增益向量L(t)使参数估计误差 $\hat{θ}$ (t)最小。由于v(t)是白噪声,假设L(t)和φ(t)与v(t)独立,对式(8)两边取期望得 $E [\hat{θ} (t)] = 0 (E$ 为期望算子)。定义参数估计误差协方差阵.

P(t): $= E [\hat{θ} (t) \hat{θ}^{Τ} (t)]$ (9)

将式(8)代入式(9)可得

P(t)=[I-L(t)φT(t)]P(t-1)×

[I-L(t)φT(t)]T+L(t)σ2LT(t)(10)

因为P(t)为非负定矩阵,将P(t)配成下列形式:

把P(t)看作待定增益向量L(t)的函数。通过极小化估计误差协方差矩P(t),可求得最优增益向量L(t).上式中矩阵P(t)包含3项,第一、二项与L(t)无关,第二项中,因为P(t-1)是非负定矩阵,所以σ2+φT(t)P(t-1)φ(t)≥0。如果选择增益L(t)使得第三项为零,即取

L(t)=P(t-1)φ(t)[σ2+φT(t)P(t-1)φ(t)]-1(12)

那么协方差阵P(t)最小,也就是参数估计误差最小。

$E [∥ \hat{θ} (t) ∥^{2}] = t r [Ρ (t)] = \min$ 。

将式(12)代入式(11)得到

P(t)=P(t-1)-P(t-1)φ(t)×

[σ2+φT(t)P(t-1)φ(t)]-1φT(t)P(t-1)=

[I-L(t)φT(t)]P(t-1)(13)

式(6),式(12)和式(13)构成了基于卡尔曼原理的单输入单输出系统的最小方差递推最小二乘算法它与标准递推最小二乘算法式(3)-式(5)的差别在于式(12)-式(13)中取σ2=1.

下面讨论一种特殊情况。对于确定性系统,噪声方差σ2=0,则算法式(6),式(12)-式(13)退化为

$\hat{θ} (t) = \hat{θ} (t - 1) + L (t) [y (t) - φ^{Τ} (t) \hat{θ} (t - 1)]$ ,(14)

L(t)=P(t-1)φ(t)[φT(t)P(t-1)φ(t)]-1(15)

P(t)=[I-L(t)φT(t)]P(t-1)(16)

为了防止式(15)中分母为零,可加上一个小常数ε(如取ε=10-6),则算法式(14)-式(16)可修改为

$\hat{θ} (t) = \hat{θ} (t - 1) + L (t) [y (t) - φ^{Τ} (t) \hat{θ} (t - 1)]$ (17)

$L (t) = \frac{Ρ (t - 1) φ (t)}{ε + φ^{Τ} (t) Ρ (t - 1) φ (t)}, 0 ＜ ε \leq 1$ (18)

P(t)=[I-L(t)φT(t)]P(t-1)(19)

在后面的例子中,将研究ε取不同值时,参数估计精度的变化。

这是文献[2]状态估计算法推导过程,推导递推最小二乘参数估计的一种新的简便方法。从上述推导看:这种参数状态估计误差具有最小方差性质。也可理论上证明这个递推最小二来算法的参数估计误差收敛于零。

2.2多输入多输出系统的递推最小二乘算法

上述方法也叮推广到多输入多输出(MIMO)系统.考虑下列多输入多输出系统,

A(z)y(t)=B(z)u(t)+v(t)(20)

A(z):=I+A1z-1+A2z-2+…+Anaz-na,

B(z):=B1z-1+B2z-2+…+Bnbz-nb,

u(t)∈Rr为系统输入向量,y(t)∈Rm为系统输出向量,u(t)∈Rm为零均值、协方差阵为Rv=E[v(t)vT(t)]∈Rm×m的随机噪声向量。

定义参数矩阵ϑ和信号向量φ(t)如下

ϑT:=[A1,A2,…,Ana,B1,…,Bnb]∈Rm×m,φ(t):=[-yT(t-1),…,-yT(t-ma),tT(t-1),…,uT(t-nb)]T∈Rn。

则式(20)可等价写为

y(t)=ϑTϕ(t)+v(t)(21)

参照式(6),由上式可知ϑ的递推估计应取为下列形式:

$\hat{ϑ} (t) = \hat{ϑ} (t - 1) + L (t) [y^{Τ} (t) - ϕ^{Τ} (t) \hat{ϑ} (t - 1)]$ (22)

则参数估计误差可写为

$\begin{array}{l} \hat{ϑ} (t) ∶ = \hat{ϑ} (t) - ϑ = [Ι - L (t) ϕ^{Τ} (t)] \times \\ \hat{ϑ} (t - 1) + L (t) v^{Τ} (t) . \end{array}$

定义参数估计误差协方差短阵P(t): $= E [\hat{ϑ} (t) \times \hat{ϑ}^{Τ} (t)]$ 。令ϵ:=tr[Rv]=E[‖v(t)‖2].采用与上节相同的方法,可以得到多输入多输出系统的最小方差递推最小二乘算法:

$\hat{ϑ} (t) = \hat{ϑ} (t - 1) + L (t) [y^{Τ} (t) - ϕ^{Τ} (t) \hat{ϑ} (t - 1)]$ (23)

$L (t) = \frac{Ρ (t - 1) ϕ (t)}{ϵ + ϕ^{Τ} (t) Ρ (t - 1) ϕ (t)}$ (24)

P(t)=[I-L(t)ϕT(t)]P(t-1)](25)

如果式(24)中取ϵ=1,就得到多输入多输出系统的标准递推最小二乘算法.

3仿真例子

例1 考虑下列随机系统,

A(z)y(t)=B(z)u(t)+v(t),

A(z)=1+a1z-1+a2z-2=1-1.60z-1+0.80z-2,

B(z)=b1z-1+b2z-2=0.412z-1+0.309z-2,

θ=[a1,a2,b1,b2]T=[-1.60,0.80,0.412,0.309]T。

仿真时,输入u(t)采用零均值单位方差不相关可测随机序列,v(t)采用零均值方差为σ2的白噪声序列.当噪声方差为σ2=0.102和σ2=1.002时,对应的噪信比分别为δns=14.26%和δns=142.63%.应用提出的算法式(6),式(12)-式(13)估计这个系统参数,不同噪声方差下的参数估计及其误差如表1所示,参数估计误差δ: $= ∥ \hat{θ} (t) - θ ∥ / ∥ θ ∥$ 随t变化曲线如图1所示.

当噪声方差为σ2=0.102,噪信比为δns=14.26%时,利用算法式(17)-式(19)来估计这个系统参数,不同ϵ下的信真结果如表2所示.

从表1和图1可知,在不同噪声水平下,噪信比越小(噪声方差小),参数估计精度高,参数估计误差随数据长度增加总的趋势是减小的。从表2可知,在同一噪声水平下,算法参数估计对ϵ取值大小不敏感。

例2 考虑下列确定性系统,即例v(t)=0情形。输入v(t)采用零均值单位方差不相关随机序列,应用算法式(17)-式(19)估计其参数,不同ϵ值时的仿真结果如表3所示。

4结论

本文基于卡尔曼滤波原理,通过极小化参数估计误差协方差矩阵的方法,分别推导了单输入单输出系统和多输入多输出系统的递推最小二乘算法,并用仿真例子研究了参数估计的性质。

摘要：基于卡尔曼滤波原理,通过极小化参数估计误差协方差矩阵,导出了递推最小二乘辨识算法。仿真例子说明了算法的有效性。

关键词：最小二乘,参数估计,参数估计误差,协方差

参考文献

[1]谢新民,丁锋.自适应控制系统.北京:清华大学出版社,2002

[2]str om KJ.Introduction to stochastic control theory.NewYork:Ac-ademic Press,1970

基于最小二乘法离场航迹构造方法篇5

飞机离场航迹是飞机起飞过程的形象化体现。目前,飞行航迹的构造方法通常采用两类方式,第一类是在飞行计划确定及气象资料完整的情况下,结合飞行动力学和运动学模型正向推导的方法。飞机起飞航迹计算方法研究提出了对不同机型统一的离场航迹构造方法,该方法主要针对离场航迹剖面进行了构造,缺少对水平投影面构造的方法;基于ANP数据库的飞机起飞仿真研究是基于详细的飞行资料和性能参数的前提下,提出了飞机离场剖面航迹构造方法;离场航迹降噪优化设计的多目标智能方法是一种利用航段飞行特征约束求解离场航迹的方法。第二类是在拥有较为准确的雷达位置信息点的情况下对雷达数据去噪,拟合出最佳函数匹配从而得到平滑的航迹。经纬仪目标交汇测量及航迹曲线拟合文中提出根据不同时刻的坐标,用最小二乘法对目标航迹进行拟合,从而推测下一时刻的位置速度及加速度;三维航迹的B样条曲线拟合算法利用B样条曲线的几何性质,解决了飞行器三维航迹拟合中的边界条件等约束问题。第二类多用于飞机离场结束后航迹的拟合。上述提出的两类方法用于离场航迹的构造存在以下三种问题:一是由于飞机离场属于低空飞行,雷达捕捉飞行器在低空飞行的位置信息不准确,飞机离场的雷达点相比于真实点误差较大,且飞机离场的方式不同,导致无法单一的利用拟合离场雷达数据的方法确定离场航迹。二是现有方法多为对离场剖面航迹进行构造,忽略由于离场方式的不同导致水平面航迹存在误差。三是在离场数据、飞行计划及气象资料缺失的情况下,无法对飞机离场航迹进行构造的问题。由此本文为了解决上述问题,采用最小二乘法结合两种离场方式特征,提出离场航迹逆向构造方法。

递推最小二乘法篇6

目前,我国10~35 k V系统以中性点不接地或经消弧线圈接地(称为小接地电流系统)的方式为主。中性点不接地方式的显著优点是,系统在发生单相接地故障(约占60%以上)时,可以短时继续运行。但是,此时由于电容电流的存在,就有可能引发间歇性弧光接地、铁磁谐振、中性点位移过电压等,甚至发展为事故。随着城市建设及城市电网改造,城市变电站大量采用了电力电缆送电线路,而在相对落后的农村地区,虽然农网改造已经大大地改善了当地的供电条件,但仍存在配电系统供电半径过大等问题。这些情况都会造成变电站不接地系统电容电流过大,弧光过电压时有出现,严重影响了系统和人身安全以及对用户的可靠供电。

电网电容电流的大小是决定是否装设消弧线圈以及消弧线圈调谐的依据,为了掌握电容电流增长变化情况以及正确设定消弧线圈补偿的档位,准确测量系统电容电流是非常必要的。

电容电流的测量可采用直接法和间接法。直接法由于对试验人员和配网系统均存在很多不安全的因素,现在几乎不再采用了。目前广泛采用的是分相外加电容法和变压器中性点外加电容法等间接法。但是这些测量方法都要接触到高压一次侧,且存在操作繁琐、准备工作耗时长、测量工作效率低等缺点[1]。

基于上述原因,本文介绍一种采用参数辨识原理测量电容电流的新方法。新方法采用注入信号法与递推式最小二乘法辨识相结合的方案,从PT二次侧进行采样,通过输入和输出间的数学模型,辨识计算出电网电容电流。

1 辨识算法的测量原理

1.1 注入信号法

首先,需要明确电容电流辨识的关键——信号注入法与最小二乘法辨识之间的关系。最小二乘法是一种系统辨识方法,用在这里是为了辨识电力系统的接地电容。然而,只有在一个系统模型中包括进去电力系统的接地电容,才能运用最小二乘法进行包括接地电容在内的系统参数的辨识。这里面的电力系统模型就是在信号注入法的实现方案中得到体现的。信号注入法运用到电压互感器PT,就包括了对电压互感器PT回路的电路等效与建模。运用这个电路模型作为目标系统模型,用注入信号角频率ω作为系统输入,用测量到的回路阻抗Z和相角θ作为系统输出,那么最小二乘法就能够对此系统发挥作用,有效地辨识到电压互感器PT回路的电路参数,得到接地电容,进而求得电容电流[2]。具体方案设计如图1所示。

图中WA、WB、WC分别为电压互感器(PT)三相的高压绕组,二次绕组Wa、Wb、Wc组成开口三角形;CA、CB、CC为三相导线对地电容。若在PT开口三角端注入一恒定电流i0,就会有3个大小相等、相位相同的电流i1、i2、i3从PT的高压侧流出,这3个电流将分别在PT三相的绕组电阻R、漏抗XL和导线对地电容上产生压降。根据图2所示的PT等值电路,可得:

一般地,三相PT的参数(绕组电阻R和漏抗XL)是对称的,而且三相导线对地电容CA、CB、CC也是基本相等的,因此三相电流i1、i2、i3分别在三相PT与导线对地电容中产生的压降是基本相同的,即uA=uB=uC,这时在PT开口三角端就可以测到一零序电压u0。

将式(2)代入式(1),得:

式(3)中:i0、u0为输入和输出变量;R、L和C为待辨识参数。这样,问题就转化为参数辨识问题[3]。

1.2 递推式最小二乘法

递推式最小二乘法的计算方法是这样的。一般最小二乘法计算系统参数的公式为

那么,在进行第N+1次数据观测后,有

经过推导可以将递推算法的公式总结为:

其中:Pi=(XiTXi)-1,γi+1=(1+XTi+1Pi Xi+1)-1,KN+1=γN+1PNxN+1为增益矩阵[4]。

式(3)描述了系统输入i0,系统输出u0与系统参数R,L,C的关系。但是其中电气量的微分,积分算子,不易进行数值化和系统实现[5]。所以,这里对公式(3)进行了变形。首先,对式(3)由时域变换到频域,微分算子d/dt变换为jω,积分算子∫dt变换为1/jω,得到:

式中:是PT开口端的电压频域变量;是PT开口端的电流频域变量;ω是注入信号的角频率。接着,对上式移到左边,公式两边取虚部,进一步变形为

其中:为阻抗模,θ为阻抗相角,n1、n2分别是电压互感器PT的高、低压绕组匝数。

这里信号注入法的测量对象是ω,U0和I0。通过测量PT开口端特定频点处的U0、I0就可以求得Z和sinθ。获得PT等效回路的电路模型,测出了ω,Z和sinθ,接下来的工作就是运用最小二乘法对电路模型中的L,C系统参数进行辨识,得到电力系统接地电容。递推式最小二乘法辨识的目标是获得系统参数L,C,记为Φ=(L,C),系统输入为注入信号的角频率ω和其倒数1/ω,记为X=(ω,1/ω),系统输出为等效阻抗的模Z和相角θ,记为Y=Zsinθ。递推式最小二乘法的算法流程如图3所示。

递推过程是这样的:在PT开口三角侧注入角频率ω逐渐变化的电流信号。在10~200 Hz的范围里,经过大量仿真实验和计算比较,最终选择扫频范围10~181 Hz之间的20个频点,各频点间隔9 Hz,作为注入信号的频率。不断地读取系统输入X和相应的输出Y,i时刻读取的一对输入、输出记为(Xi,Yi)。每次读取一对(Xi,Yi),用最小二乘法计算出当前迭代的系统参数iΦ。如此循环,直到系统参数Φn与系统参数Φn+1之差的绝对值小于一个非常小的正数ε,那么迭代结束,取得系统辨识参数Φn=(Ln,Cn)。最后,计算出电网电容电流Ic=3ωCU相。

2 电容电流辨识算法实现

最小二乘法辨识电容电流的实现主要是完成辨识算法程序,并将其加载到DSP芯片中。辨识电容电流程序的目标器件选用TMS320LF2407芯片,采用C语言进行开发,不仅能够方便地实现算法的编程,而且也能够达到较高的辨识速度,不会影响电容电流辨识的实时性。在前面设计的辨识算法的基础上,可以比较容易地制定出辨识电容电流的程序流程,如图4所示。

其中,程序运行的终止条件可以根据系统要求和实际情况而定。如果对辨识电力系统电容电流的反应时间有规定,那么可以设定迭代次数不超过N。经过多次仿真实验和计算,这里N一般可取10,就可以使系统辨识结果具有较高精度。如果对辨识电力系统电容电流的精度有较高要求,那么可以设置辨识系统参数Φn和下一次的参数Φn-1之差的绝对值小于一个非常小的正数ε,那么迭代结束。一般地,ε可以取10-4,所用的辨识迭代次数通常在10次以内。另外,最小二乘法程序中涉及到对矩阵的较多操作,比如辨识算法中PiXi+1是n×n方阵和1×n向量相乘,(I-Ki+1XTi+1)P是两个n×n方阵的乘积,Pi=(XiTXi)-1需要进行矩阵求逆运算。程序中对各种类型的矩阵运算进行了单独编程,并且用函数进行实现,从而能将辨识过程进行明确划分,更好地执行辨识流程。

3 电容电流辨识算法仿真

这里,我们给出CCS2.20中TMS320LF2407芯片完成电容电流辨识仿真的过程和结果。CCS2.20是TI推出的用于开发其DSP芯片的集成开发环境。它采用Windows风格界面,集编辑、编译、链接、软件仿真、硬件调试及实时跟踪等功能于一体,极大地方便了DSP程序的设计与开发。对电压互感器PT回路的建模采用Matlab实现,仿真并获得PT开口端的等效回路的各次电流Ii和电压Ui,从而计算得到各个角频率相应的电路阻抗值Zsinθ。根据各次电流Ii和电压Ui,在CCS2.20中仿真运行递推式最小二乘法算法程序,然后我们就可以在CCS2.20中设置断点,监测到程序辨识到的接地电容,从而验证DSP的电容电流辨识程序的正确性。

在PT开口端输入10~181 Hz,频率间隔为9Hz,幅度为1 A,相位为0°的20个电流信号Ii。在每输入一个电流信号之后,根据PT开口端测量到的电压值,就可以计算等效阻抗值Zsinθ。各个角频率对应的电路阻抗值Zsinθ见表1所示。

然后,把电流注入法检测到的各次电路等效阻抗Z和相角θ作为输入在CCS2.20中运行递推式最小二乘法程序,检测接地电容C。表1中列出了辨识的结果。由表1可见系统对电感L的辨识误差从第三迭代开始都在1%以内,而对电容C的辨识误差则在第一次辨识就达到1%以内。

由此,我们可以看到,这里实现的递推式最小二乘法辨识电容电流具有较高的辨识效率,而且相比一般最小二乘法,可以使系统计算开销降低,辨识速度提高。

4 结论

本文实现的接地电容电流辨识的实现方案具有一定创新性。这里辨识算法首先对电力系统加入了电压互感器PT,并且建立了电力系统的等效模型和系统方程。辨识所需的电力系统状态的获取是通过信号注入法获得的,而对系统参数的辨识中使用的是递推式最小二乘法。本文在接地电容电流辨识中创新性地对PT等效回路的电路方程进行了变换,使得其中的系统输入、输出变量的采集、获取比较容易实现。此外,本文电路系统参数的求解运用了递推式最小二乘法,它的特点是能够实时、在线辨识系统参数,满足消弧线圈实时补偿电力系统电容电流的要求。

参考文献

[1]孙结中,蒋学军.配电系统电容电流测试仪的研制[J].广西电力,2006(2):1-4.SUN Jie-zhong,JIANG Xue-jun.Development of Capacitive Current Tester in Distribution System[J].Guangxi Electric Power,2006(2):1-4.

[2]赵正军,姜新宇.信号注入法在配电网电容电流测量中的研究[J].广东电力,2005,17(6):25-28.ZHAO Zheng-jun,JIANG Xin-yu.Research of Signal Injection Method used in the Measurement of Capacitive Current in Distribution Network[J].Guangdong Electric Power,2005,17(6):25-28.

[3]徐宁寿.系统辨识技术及其应用[M].北京:机械工业出版社,2003.XU Ning-shou.System Identification Technology and Its Application[M].Beijing:China Machine Press,2003.

[4]陈自宽,母国光.相关数据集的最小二乘处理方法[J].数据采集与处理,2007,11(1):66-68.CHEN Zi-kuan,MU Guo-guang.Least Square Method Sets to Deal with the Related Data[J].Data Acquisition and Processing,2007,11(1):66-68.

递推最小二乘法篇7

1 北票市气候背景

北票市属温带半湿润半干旱的大陆性季风气候, 年平均降水量450~500 mm, 无霜期130~150 d。由于受蒙古高原大风的侵袭, 干旱、高温现象十分严重。北票市气象局多年统计资料显示, 历年气温极端最高值出现在2000 年7 月14日, 为42.3 ℃, 气温极端最低值出现在2001 年1 月13 日, 为-28.2 ℃。

2 统计学角度订正方法

2.1 最小二乘法的数据处理

在利用日本传真图预测本地乡镇温度中, 如何进行数据处理, 减小系统性误差, 是实现气象业务精准化的关键所在。最小二乘法通过最小化误差的平方和寻找数据的最佳函数匹配, 是一种数学优化技术, 可以利用计算机编程的形式处理大量数据, 免除了人工计算的繁琐, 结果也更为准确。目前北票市乡镇预报对温度的订正主要采用最小二乘法[1]。

2.2 气温订正方案

本文使用的数据为2013 年5—6 月北票市下府开发区的预测与实况资料, 通过最小二乘法对6 月日最高气温和最低气温进行订正, 进而对5 月日最高气温与最低气温进行效果检验。在利用最小二乘法对数据进行处理中, 把预报值选作x, 实测值作为y, 因此所有的误差只是y的误差[2]。

将数据导入spss分析软件中, 对自变量x与因变量y做线性回归分析并绘制散点图 (图1) , 可知气温预报值与实况值之间存在一定的线性关系, 可以通过最小二乘法对其进行订正。

设一元线性拟合公式:y=a0+a1。现在已知m个试验点xi, yi (i=1, 2, ..., m) , 求2 个未知参数a0, a1。

由最小二乘法原理, 参数a0, a1应使

取得极小值。根据极小值的求法, a0和a1应满足:

这就是含有2个未知数和2个方程的正规方程组。

从中解得a0, a1, 即:

其中

将数据代入公式可得最高温度与最低温度的订正公式:

最高温度:Ymax=0.886Xa+3.981

最低温度:Ymin=0.570Xi+6.640

2.3 乡镇预报订正效果的检验

为了检验线性回归统计方法来制作城镇气温预报结果的实际效果, 将订正后的气温预报结果与预报员预报结果做比较。最高气温订正前后与实况值的效果对比如图2 所示, 可以看出, 订正结果对预报效果的改善是明显的。最低气温订正前后与实况值的效果对比如图3 所示, 可以看出, 订正后的结果对预报预测效果有较好的可信度[3,4]。

3 结论

采用线性回归统计方法订正未来24 h乡镇气温预报方法, 是具体天气预报业务中应用的方法之一。应用结果表明其可用性较强, 但与预报员的水平还存在一定的差距。该文所用资料是夏季气温, 其偶然性较多, 冬季的效果可能更好。

参考文献

[1]余建英, 何旭宏.数据统计分析与SPSS应用[M].北京:人民邮电出版社, 2001.

[2]黄治勇, 张文, 陈璇, 等.湖北省乡镇温度预报方法初探[J].气象, 2011 (12) :1578-1583.

[3]吴丹, 王优, 冯雪君, 等.基于统计分析的龙城地区乡镇温度预报方法探讨[J].现代农业科技, 2015 (4) :341-342.

递推最小二乘法篇8

由于岩土介质都是不均匀、不连续的,在其中传播的波动现象也非常复杂,波的反射与折射、介质体内的内摩擦导致能量不断被吸收而发生能量散耗现象,使得地震波向外传播的过程是一个指数不断衰减的过程。

由于爆源的复杂性,以及传播介质体物理力学性质和地质结构构造的多样性,使得爆破地震波具有随机波非重复性的特性,因此,爆破地震波是一种随机波。爆破地震波不仅在振幅上随时间做复杂的变化,而且波的频率构成和持续时间等也随环境条件、爆心距、爆破规模以及地层地质等的影响表现为极为复杂的现象。

1 回归分析原理

利用《爆破安全规程》预测爆破地震动强度的萨道夫斯基的经验公式,并利用最小二乘法对其进行回归分析,其公式的具体形式是:

对式(1)两边同取自然对数得:

得到:

式(3)中x是可控或能够精确测量的变量,若使x取一组不完全相同的值xi(i=1,2,3,…,n),进行计算分析,就得到与之对应的观测值yi(i=1,2,3,…,n)。(x1,y1),(x2,y2),…,(xn,yn)就是一组样本。显然可以通过这组样本估计y的数字特征。其中数学期望μ=μ(x)是重要的数字特征,函数μ(x)称为y对x的回归。如果y与x的关系是线性的,则称为一元线性回归。一元线性回归的数学模型:

其中,εi(i=1,2,3,…,n)为一组相互独立且服从于同一正态分布N(0,σ2)的随机变量。在上述条件下,y是服从正态分布N(βo+βx,σ2)的随机变量,在一元回归模型下有:

称它为y对x的线性回归方程,其图形称为回归直线。其中βo和βx称为回归方程的回归系数。求βo,βx的估计量bo和bx,也就是要确定一条经验回归直线y°=bo+bx来近似表达y与x的关系以下利用最小二乘法进行分析

对于每一个xi,都可由式y°=bo+bxi确定一个回归值y°i,回归值与观察值之差yi-y°i=yi-bo-bxi刻画了与回归直线的偏离程度。显然,对于所有的xi,若y°i与yi的偏离越小,则就认为回归直线和所有试验点拟合的越好。显然,全部观察值yi与回归值y°i的偏差平方和刻画了全部观察值与回归直线的偏离程度。

最小二乘法是使得偏差平方和Q(bo,b)达到最小的一种确定bo和b的方法。因此求bo和b就成为求Q(bo,b)最小值的问题。由于Q(bo,b)是bo和b的非负二次函数,所以它存在最小值。根据微积分学中的极值原理,要求的估计值bo和b是下列方程组的解:

式(6)方程组称为正规方程组,它可转化为:

解方程组得:

2 监测数据分析

以下利用所测得的数据对公式进行回归分析,建立爆破地震振动速度与爆心距之间的线性关系。选取某地下工程施工时对现场进行的部分监测数据。

根据表1中的数据,利用最小二乘法回归分析的原理进行线性回归,求得:K=186.49,α=1.434 7。故可得出爆破时振动速度沿不同方向传播衰减规律,如下所示:

其中,V为测点处的地表面质点振动速度,cm/s;Q为单段最大装药量为传播距离

在这一地下工程爆破中,利用最小二乘数法能够容易得出振动速度与装药量、传播距离之间的具体关系。并且利用式(9)对以后爆破进行指导都能起到较理想效果。故利用最小二乘数回归分析的爆破经验公式的方法是可行的

3 结语

本文通过对某地下工程监测数据利用最小二乘法对萨道夫斯基的经验公式进行回归分析得出其具体参数值,并且能够在具体工程中得到很好的利用。因此说明该方法能很好地解决爆破震动参数的回归分析进一步为工程提供爆破理论依据

摘要：针对爆破的不确定性与难估量性以及传播介质特性的复杂性,利用最小乘法原理对爆破震动进行回归分析,得出爆破震动速度与装药量、传播距离的关系,实践证明该方法能很好地解决爆破震动参数的回归分析。

关键词：爆破,最小二乘法,回归分析,爆破震动