利用传统的数学建模和卡尔曼滤波器对大型土木结构抗震性能的高频数据在线去噪（一）

       正如前面总结的那样，我们先从实际的应用出发，来探索传统的算法，并对传统算法进行改进，我们首先选取几个完全可以利用传统算法解决的例子，到最后我在慢慢的讲到利用传统算法很难去完成的例子，其实通过下面的几章的介绍大家可以发现，传统的算法，通过数学建模可以去解决很多复杂的问题，并不一定要利用人工智能，我的观点还是那样，如果可以利用传统算法去解决的，我们尽量通过建模去完成，因为“白箱子”总会比“黑箱子”更加直观（还有一些原因在上一章介绍过来，这里就不在提了），如果确实无法利用传统算法了，我会在考虑引入深度学习的知识，去解决，还希望有大牛提出更加深刻的意见。
       去年有幸接触到这样的一个项目，就是大型土木结构的抗震性能分析，具体的就是一个振动台，实物图不太好发~见谅哈，我们会利用传感器测试到试件的摆动位移，但是，这个数据还有巨大的噪声，而且由于地震过后争分夺秒，我们必须在线的对其数据进行去噪，然后拿到这些比较可靠的信息去做分析，对其算法的速度和精度都有着很高的要求（必定关乎人命呀）。 


       接下来先看看测量数据（实测数据）和参考数据吧~如下图，可以看得出来，信噪比不忍直视啊，而且参考数据都已经是高频数据了，噪声还是有色噪声，就算没有在线处理的要求，也很难用频域的方法解决（我也参阅了大量的文献（上百篇），最近两年的SCI检索文献也复现过很多（10篇左右），频域的方法很容易丢波峰，而且很难做到在线处理）。


       再看看最后在线处理的结果，可以看得出来，在完全保留了波峰的情况下，对有色噪声的处理也非常的到位，而且可以做到完全的在线，几乎没有任何延迟。


       这是如何实现的呢，在这里我介绍一下我的思路。卡尔曼滤波器最核心的就是通过过程模型得到预测值以及通过预测值和测量模型得到最后的估计值。首先，对于最基本的卡尔曼滤波器来说，是完全无法解决上述问题的，因为经典的卡尔曼滤波器只能对线性的含有白噪声的信号进行滤波，而要想从含有大量有色噪声的信号中提取出有用的信息，我们需要至少考虑两点：第一，高频的参考数据的统计特性；第二，有色噪声的统计特性。但是，其实这两点对于我们来说都是未知的，我们并不知道信号的统计特性，更不知道有色的统计特性，因此，我们需要对其进行建模。（在这里稍微提一下，卡尔曼滤波算法是全局最优的，也就是说在模型足够准确的情况下，卡尔曼滤波可以说是最准确的估计算法，因此建模才是最关键的）。        在这一章我会给出最简单的基于自适应过程模型的算法流程（后面还会有很多的改进，这里并没有涉及到对测量噪声的建模，只涉及到了对过程模型中状态转移矩阵和过程噪声方差矩阵进行建模），这里我们基于“数据中含有模型而模型融于数据”的基本思想，一边利用自适应模型得到更准确的估计值，一边又利用更准确的估计值得出更准确的模型，使其忽略初值的影响，慢慢收敛。其算法主要由两个部分组成，状态估计和参数更新，流程图如下：

后面给出了详细的公式（上图的结果运用的算法还有比较多的改进，这里第一章还是介绍简单一点的，详细的介绍请参考‘Yi S L, Jin X B, Su T L, et al. Online Denoising Based on the Second-Order Adaptive Statistics Model[J]. Sensors, 2017, 17(7).‘）

1.1状态估计

第一步，根据维纳-辛欣定理对有色噪声做白化处理，得到：    

                                (1.1-1)
将上式写为连续时间系统状态方程为： 

     (1.1-2)
即状态转移矩阵和噪声转移矩阵 分别为：   

               (1.1-3)
     第二步，以周期 k采样，离散化后系统目标运动满足以下方程：    

       (1.1-4)                           
随着系统自适应参数的变化而变化；状态转移矩阵的表达式如下：   

             (1.1-5)
过程噪声的方差矩阵的表达式如下：    

    (1.1-6)
第三步，根据建立的具有系统自适应参数的运动模型和初始值完成目标状态的一步预测，预测方程式如下：    

                   (1.1-7)         
接着，按照下式完成目标状态协方差的一步预测：

          (1.1-8)
    第四步，根据目标的状态预测值、观测数据值和目标的状态协方差预测值对目标状态进行更新。
首先，根据目标状态协方差预测值、测量矩阵及测量噪声方差按照下式计算滤波器增益：  

   (1.1-9)
然后，利用目标状态预测值和观测数据值计算目标当前状态估计值：  

     (1.1-10)
    最后，按照下式计算目标状态协方差的估计值：   

         (1.1-11)

1.2参数更新

    第一步，根据目标状态估计值计算目标加速度估计值
    首先，利用下式获取系统k时刻和k+1时刻加速度估计值     

                     (1.2-1)     

            (1.2-2)
    第二步，根据加速度估计值对系统自适应参数进行修正
根据采样时刻k值的大小，选择修正系统自适应参数方法，若k小于等于4进入第三步，若k大于4进入第四步，
第三步，当采样时刻k小于等于4时，因为采样数据较少,按下式计算系统自适应参数，



第四步，当采样时刻k大于4时，利用Yule-Walker方法，按下式计算系统自适应参数，  

            (1.2-3)
根据系统方程得到目标运动加速度，方程满足如下一阶马尔科夫随机序列：       

                 (1.2-4)
一阶马尔科夫时间加速度序列满足以下参数关系：      

              (1.2-5)
其中， 与 分别为加速度的机动频率及其离散化后加速度序列的机动频率，自适应参数 和 可按照下式计算得到：                                   

                   (1.2-6)    

              (1.2-7)

具体还想了解具体算法的小伙伴可以去查阅我的论文（程序比较复杂，如果需要相关程序的也可以找我）。