滤波算法影响

1. 什么是滤波算法

卡尔曼滤波器（Kalman Filter）是一个最优化自回归数据处理算法（optimal recursive data processing algorithm）。对于解决很大部分的问题，他是最优，效率最高甚至是最有用的。他的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。

最佳线性滤波理论起源于40年代美国科学家Wiener和前苏联科学家Kолмогоров等人的研究工作，后人统称为维纳滤波理论。从理论上说，维纳滤波的最大缺点是必须用到无限过去的数据，不适用于实时处理。为了克服这一缺点，60年代Kalman把状态空间模型引入滤波理论，并导出了一套递推估计算法，后人称之为卡尔曼滤波理论。卡尔曼滤波是以最小均方误差为估计的最佳准则，来寻求一套递推估计的算法，其基本思想是：采用信号与噪声的状态空间模型，利用前一时刻地估计值和现时刻的观测值来更新对状态变量的估计，求出现时刻的估计值。它适合于实时处理和计算机运算。

现设线性时变系统的离散状态防城和观测方程为：

X(k) = F(k,k-1)·X(k-1)+T(k,k-1)·U(k-1)

Y(k) = H(k)·X(k)+N(k)

其中

X(k)和Y(k)分别是k时刻的状态矢量和观测矢量

F(k,k-1)为状态转移矩阵

U(k)为k时刻动态噪声

T(k,k-1)为系统控制矩阵

H(k)为k时刻观测矩阵

N(k)为k时刻观测噪声

则卡尔曼滤波的算法流程为：

预估计X(k)^= F(k,k-1)·X(k-1)

计算预估计协方差矩阵
C(k)^=F(k,k-1)×C(k)×F(k,k-1)'+T(k,k-1)×Q(k)×T(k,k-1)'
Q(k) = U(k)×U(k)'

计算卡尔曼增益矩阵
K(k) = C(k)^×H(k)'×[H(k)×C(k)^×H(k)'+R(k)]^(-1)
R(k) = N(k)×N(k)'

更新估计
X(k)~=X(k)^+K(k)×[Y(k)-H(k)×X(k)^]

计算更新后估计协防差矩阵
C(k)~ = [I-K(k)×H(k)]×C(k)^×[I-K(k)×H(k)]'+K(k)×R(k)×K(k)'

X(k+1) = X(k)~
C(k+1) = C(k)~

2. 扩展卡尔曼滤波（EKF）算法详细推导及仿真（Matlab）

姓名：王柯祎

学号：20021110373T

转自 ：https://blog.csdn.net/gangdanerya/article/details/105105611

【嵌牛导读】介绍扩展卡尔曼滤波（EKF）算法的详细推导，局限性和MATLAB仿真。

【嵌牛鼻子】扩展卡尔曼滤波（EKF）

【嵌牛正文】

扩展卡尔曼滤波算法 是解决非线性状态估计问题最为直接的一种处理方法，尽管EKF不是最精确的”最优“滤波器，但在过去的几十年成功地应用到许多非线性系统中。所以在学习非线性滤波问题时应该先从EKF开始。

EKF算法是将非线性函数进行泰勒展开，然后省略高阶项，保留展开项的一阶项，以此来实现非线性函数线性化，最后通过卡尔曼滤波算法近似计算系统的状态估计值和方差估计值。

一、EKF算法详细推导

【注】EKF推导参考的是黄蔚的博士论文“CKF及鲁棒滤波在飞行器姿态估计中的应用研究”，论文中EKF，UKF和CKF等算法讲解的都很详细，值得一看。

我们把KF与EKF算法拿出来对比可以发现：

二、EKF算法局限性：

该算法线性化会引入阶段误差从而导致滤波精度下降，同时当初始状态误差较大或系统模型非线性程度较高时，滤波精度会受到严重影响甚至发散。

需要计算雅克比矩阵，复杂，计算量大，影响系统的实时性，还会导致EKF算法的数值稳定性差。

当系统存在模型失配，量测干扰，量测丢失，量测延迟或状态突变等复杂情况时，EKF算法鲁棒性差。

三、Matlab仿真：

clear all;clc;   close all;

tf = 50; 

Q = 10;w=sqrt(Q)*randn(1,tf); 

R = 1;v=sqrt(R)*randn(1,tf);

P =eye(1);

x=zeros(1,tf);

Xnew=zeros(1,tf);

x(1,1)=0.1; 

Xnew(1,1)=x(1,1);

z=zeros(1,tf);

z(1)=x(1,1)^2/20+v(1);

zjian=zeros(1,tf);

zjian(1,1)=z(1);

for k = 2 : tf

%%%%%%%%%%%%%%%模拟系统%%%%%%%%%%%%%%%

    x(:,k) = 0.5 * x(:,k-1) + (2.5 * x(:,k-1) / (1 + x(:,k-1).^2)) + 8 * cos(1.2*(k-1)) + w(k-1); 

    z(k) = x(:,k).^2 / 20 + v(k);

%%%%%%%%%%%%%%%EKF开始%%%%%%%%%%%%%%%

    Xpre = 0.5*Xnew(:,k-1)+ 2.5*Xnew(:,k-1)/(1+Xnew(:,k-1).^2) + 8 * cos(1.2*(k-1));  

    zjian =Xpre.^2/20;

    F = 0.5 + 2.5 * (1-Xnew.^2)/((1+Xnew.^2).^2);

    H = Xpre/10;    

    PP=F*P*F'+Q; 

    Kk=PP*H'*inv(H*PP*H'+R);

    Xnew(k)=Xpre+Kk*(z(k)-zjian);

    P=PP-Kk*H*PP;

end

  t = 2 : tf;  

 figure;   plot(t,x(1,t),'b',t,Xnew(1,t),'r*');  legend('真实值','EKF估计值');

仿真结果：

3. 采样频率对四点加权平均滤波有什么影响

四点加权平均滤波算法是对各次采样输入值取不同的比例后再相加.一般,次数愈靠 后,控制系数(比例)取愈大,这样,最近一次采样输入值影响愈大.

4. 数学滤波算法可以处理三个坐标点吗

滤波算法可以处理三个坐标点的。滤波在三坐标中的应用：
1、粗糙度对测量的影响：测量点也在图中被放大获取到大量的点，表面粗度被认为是，引起“噪点”的原因。
2、探针的机械滤波：
选择探针直径-使用探针测量工件会由于工件表面结构的影响产生机械滤波。
由于探针直径过大精细的工件表面的形状无法捕捉，因此可看作是机械低通滤波。
3、三坐标的滤波：
用同样参数进行低通滤波的扫描线。
如下图所示，描绘出的图形差异并不明显。

4、2 RC滤波：不再使用圆度测量最初的标准化滤波器，但是已被现代滤波计算所取代。
5、高斯滤波：坐标测量技术中标准滤波算法。此滤波方法为标准算法被广泛使用。他使用高斯曲线加权计算测量点得到新的轮廓。
6、样条滤波：基于滤波方程的增强滤波方法（多项式计算），样条滤波更合乎标准，也更优于高斯滤波但并不是标准滤波方法。
(4)滤波算法影响扩展阅读：
图像滤波是一种非常重要的图像处理技术，现在大火的卷积神经网络其实也是滤波的一种，都是用卷积核去提取图像的特征模式。不过，传统的滤波，使用的卷积核是固定的参数，是由经验非常丰富的人去手动设计的，也称为手工特征。而卷积神经网络的卷积核参数初始时未知的，根据不同的任务由数据和神经网络反向传播算法去学习得到的参数，更能适应于不同的任务。
自适应中值滤波
中值滤波器是一种常用的非线性滤波器，其基本原理是：选择待处理像素的一个邻域中各像素值的中值来代替待处理的像素。主要功能使某像素的灰度值与周围领域内的像素比较接近，从而消除一些孤立的噪声点，所以中值滤波器能够很好的消除椒盐噪声。不仅如此，中值滤波器在消除噪声的同时，还能有效的保护图像的边界信息，不会对图像造成很大的模糊（相比于均值滤波器）。
中值滤波器的效果受滤波窗口尺寸的影响较大，在消除噪声和保护图像的细节存在着矛盾：滤波窗口较小，则能很好的保护图像中的某些细节，但对噪声的过滤效果就不是很好，因为实际中的噪声不可能只占一个像素位置；反之，窗口尺寸较大有较好的噪声过滤效果，但是会对图像造成一定的模糊。另外，根据中值滤波器原理，如果在滤波窗口内的噪声点的个数大于整个窗口内非噪声像素的个数，则中值滤波就不能很好的过滤掉噪声。

5. 无味卡尔曼滤波与扩展卡尔曼滤波的具体区别，以及算法

EKF是对非线性系统模型（方程）进行的线性化近似，以利用KF算法进行滤波估计。而UKF是对状态的概率统计近似，即设计少量的σ点，由σ点经由非线性函数的传播，计算出随机向量一、二阶统计特性的传播，对于高斯噪声的假设，UKF能够达到三阶估计精度，而EKF只能达到二阶精度，但其算法仍然是利用KF的算法。
现在国内外的文献大都是对UKF算法的改进和应用进行论述，但对算法的稳定性等没有系统的论述。我了解得沈阳自动化所做的这方面的工作很多。

6. 滤波算法总结

1、限幅滤波法（又称程序判断滤波法）

2、中位值滤波法

3、算术平均滤波法

4、递推平均滤波法（又称滑动平均滤波法）

5、中位值平均滤波法（又称防脉冲干扰平均滤波法）

6、限幅平均滤波法

7、一阶滞后滤波法

8、加权递推平均滤波法

9、消抖滤波法

10、限幅消抖滤波法

程序默认对int类型数据进行滤波，如需要对其他类型进行滤波，只需要把程序中所有int替换成long、float或者double即可。

原文

7. 卡尔曼滤波算法的功能是什么

卡尔曼滤波是用来进行数据滤波用的，就是把含噪声的数据进行处理之后得出相对真值。卡尔曼滤波也可进行系统辨识。卡尔曼滤波一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。由于观测数据中包括系统中的噪声和干扰的影响，所以最优估计也可看作是滤波过程。

8. 突变值对卡尔曼滤波算法的影响

突变值对卡尔曼滤波算法的影响会逐步变小。卡尔曼滤波经过多次的迭代后会逐渐缩小原始值误差，输出结果也更接近真实的值。

9. 减少噪声的匹配滤波算法

（1）传统匹配滤波算法

Rickett et al.（2001）给出了匹配滤波简要的公式及算子长度设计标准，本节给出了更为详细的匹配 滤波公式，并给出推导公式基本条件和结果。

设同一地区不同时期Y1，Y2得到的地震数据分别为G^Y1（t），G^Y2（t），取Y1年份的地震记录为参

考地震道，使Y2年份相应的地震记录与之匹配。选取归一化算子p使得目标泛函：

海上时移地震油藏监测技术

极小。最终得到关于求解匹配滤波器｛P（m），m＝1，2，…，L｝的L个方程的方程组：

海上时移地震油藏监测技术

为意义更明确，对上面的公式进一步简化，令

海上时移地震油藏监测技术

上两式中：R_Y2Y2（m－n）为时间延迟为m－n的时期Y2地震记录在设计窗口中的自相关；R_Y1Y2（n）为时间延迟为n的时期Y1与时期Y2地震记录在设计窗口中的互相关，于是方程（4.8）可以进一步写成：

海上时移地震油藏监测技术

求解方程组（4.11）得到匹配滤波器算子｛P（m），m＝1，2，…，L｝，用

海上时移地震油藏监测技术

校正相应的地震剖面。通过实际数据处理结果验证了上述推导的正确性和方法的有效性。

方程（4.11）写成矩阵形式：

海上时移地震油藏监测技术

式中：M为时期Y2地震记录在设计窗口中的自相关序列组成的Toeplitz矩阵，R为时期Y1与时期Y2地 震记录在设计窗口中的互相关序列向量。求解方程（4.13）可采用Levinson递推算法，计算效率高。

为了减少噪音的影响，通常引入阻尼项，方程（4.13）变为

海上时移地震油藏监测技术

式中：μ为很小的数，通常为可设为0.01或0.001。

实际应用中，可以发现式（4.13）受噪声的影响很大，不稳定。虽然加入阻尼项后结果有所改善，但 如何选取合适阻尼因子又是一个难题。为此推导新的匹配滤波表达形式，寻求更稳健的求解方法。

（2）新匹配滤波公式

同样设同一地区不同时期Y1，Y2得到的地震数据分别为G^Y1（t），G^Y2（t），取Y1年份的地震记录 为参考地震道，使Y2年份相应的地震记录与之匹配。则匹配过程可描述为

海上时移地震油藏监测技术

其中M为G^Y2组成的褶积矩阵。如果设地震道的采样点数为n，设计滤波器f长度为m，M则为（2×n－1）×m矩阵，为保持矩阵维数相同，一种方法是将G^Y1后面补零为（2×n－1）×1向量，另一种方法是取 矩阵M的前n×m项。如果采用第一种方法，可以验证得到的公式与（4.13）式相同。在此采用后一种方 法，得到新的匹配滤波方程。只要设计滤波器f足够长，总能满足能量差e（f）最小，根据范数定义：

海上时移地震油藏监测技术

求解能量差e（f）最小问题可转化为

海上时移地震油藏监测技术

即对滤波因子向量求导，最终可归结为求解线性方程：

海上时移地震油藏监测技术

如果记A＝M^TM，b＝M^TG^Y1，方程（4.18）转化为

海上时移地震油藏监测技术

（4.19）式形式上与（4.13）式类似，内容不同，不再是Toeplitz矩阵，因此不能应用Levinson递推算法求解。因此，引入奇异值分解方法求解方程（4.19）。

（3）基于奇异值分解的匹配滤波算法

矩阵的奇异值分解，是矩阵计算中一套很有用的技术。它可以有效地处理系数矩阵是奇异的或者接 近奇异的方程组。对于矩阵A，如果A∈R^m×n，并且A的秩为r，总有

海上时移地震油藏监测技术

其中， V为正交阵。 ，并且 为A 的奇异值。

公式（4.20）即为矩阵A的奇异值分解，根据正交矩阵的性质：

海上时移地震油藏监测技术

很容易表示出矩阵A的逆矩阵

海上时移地震油藏监测技术

将式（4.22）带入式（4.19）中，得到滤波因子的表达式为

海上时移地震油藏监测技术

实际计算中，当A是奇异阵出现奇异值，或A接近奇异或病态矩阵时，（4.23）式的计算过程就无法进行。这时可将出现的奇异项 （σ_k是零，或者数值很小）简单地替换成零或很小的常数，通过这种方法能得 到方程稳定的解。

对于实际含有噪声的信号，信号能量主要分布在奇异值大的分量上，因此去除小奇异值同时能消除 噪声影响。通常可选取某一能量百分比的奇异值作为去除的阈值，以这种方式既能克服A接近奇异或病 态矩阵的影响，又能减小噪声的影响，使滤波因子稳健。

（4）模拟数据验证

模拟得到一组存在时间、振幅、频率、相位差异的信号，作为基测线与监测测线地震道，对监测测 线地震道加入不同比例的随机噪声，组成验正算法有效性的数据体，如图4.10所示。分别用传统的匹配 滤波方法和重新推导的基于奇异值分解的匹配滤波方法进行匹配处理，比较匹配后基测线与监测测线振 幅差异，结果见图4.11和图4.12。可以看出，传统匹配滤波公式的计算结果受噪声的影响很大，而基于 奇异值分解的匹配滤波方法具有很好的抗噪声能力。

图4.10 模拟地震记录（从上至下依次为加入0％，10％，20％，30％噪声的信号）

图4.11 传统方法匹配结果

图4.12 基于奇异值分解方法匹配结果

（5）实际数据验证

选择一块同一地区两次不同时间测得的两条二维测线；选取油藏上方时间长度为300ms的窗口作为 滤波因子设计窗口，并以抽取其中139道构成验证互均衡算法的数据体（图4.13，图4.14）。分别采用 传统匹配滤波公式与基于奇异值分解的匹配滤波两种方法进行校正。比较差异剖面的平均能量，结果见 图4.15。从图中可知基于奇异值分解的匹配滤波方法具有更好的抗噪声能力，匹配误差远小于传统匹配 滤波。

图4.13 某地区时间1地震记录

图4.14 某地区时间2地震记录

图4.15 两种匹配方法结果误差能量对比图

本节推导了新的匹配滤波方程，提出基于奇异值分解的匹配滤波算法，理论和实际数据都验证了该 方法有效性。这里从计算精度上比较两种匹配滤波算法，实际处理时移地震数据时还要考虑计算时间，此时寻求快速的奇异值分解算法是一种提高处理效率的方式，另外针对不同信噪比，将传统匹配滤波算 法与基于奇异值分解的匹配滤波算法结合应用同样是一种很好的方式。总之，基于奇异值分解的匹配滤 波提高了匹配精度，有利于为时移地震解释提供一致性更好的地震资料。