最近邻滤波快速方法

㈠ 怎么通过系统的频率响应判断系统是什么类型的滤波器

可以用以下方法
根据系统函数快速判断滤波器类型 (1）死办法，用傅里叶变换求出H(f)，在画出幅频特性曲线，看高频部分是不是“通”
（2）用拉氏变换求出H(s)，然后记住一句话：分子上有什么就通什么！
举个例子：
H(s)=as/(bs+c)
分子上有“举锋高次”，所以是高通。
这里的“高次”是这个意思：
分母上有s的0次和1次，分子是s的1次，所以是较高的那个，简数答姿称“高次”。
H(s)=a/(bs+c)
分子上有“低次”，所以是低通。
H(s)=as^2/(bs^2+cs+d)
分子上有“高次”，所以是高通。
H(s)=a/(bs^2+cs+d)
分子上有“低次”，所以是低通。
H(s)=as/(bs^2+cs+d)
分子上有“中间次”，所以是带通。
第（2）种方法还没找到理论根据，如果将分子分母都除以“高次”，在判断频率从薯绝小变化到无穷的情况能理解
如果只有一个零极点，可以根据复平面上零极点位置来判断。

㈡ *维纳滤波与匹配滤波

我们知道，浅部地质体所产生的重磁异常比深部地质体产生的重磁异常要尖锐得多。一个尖锐的异常其幅值从异常中心向外快速下降，以具有很大的高频成分为特征。另一方面，宽缓的异常从中心向外是缓慢的衰减，具有集中于低频端的谱。异常频谱特征的这种差异，提供了分离浅部场和深部场的可能性。1966年，Bhattacharyya详细研究了矩形棱柱体总磁异常的连续谱，1970年Spector与Grant运用统计结构的基本假设，引入“总体平均”的概念，推导分析了航磁图的能谱公式，把关于矩形棱雹此柱体的谱的某些性质推广到块状体，讨论了块状体的水平尺寸、深度和厚度对谱的影响，提出了用能谱分析来粗略估计块状体的埋深、延深的方法。1975年，Spector运用上述方法，提出了“匹配滤波”方法，并用此方法处理科迪雷拉山区的航磁图，消除了火山岩覆盖的干扰，从而得到与成矿有关的火成岩引起的异常图。

图3-7-6 湖北铁山、鄂城岩体的ΔT异常

（单位nT）

图3-7-7 湖北铁山、鄂城岩体的ΔT异常化向地磁极

（单位nT）

（一）最小均方差滤波器与维纳滤波器

最小源缓迅均方差滤波又称为维纳滤波，其原理是设计一个滤波器，使滤波后的输出与期望输出之间的均方差为最小。设g（x）为f（x）滤波后的输出函数，s（x）为期望输出，h（x）为脉冲响应，如使总均方差：

地球物理勘探概论

则相应的滤波器称为最小均方差滤波器。根据最小均方差滤波的理论，由式（3-7-21）我们可以得到维纳滤波器的频率响应H（ω）：

地球物理勘探概论

其中：H（ω）是维纳滤波器的频率响应；P_fs（ω）为输入函数f（x）与期望输出信号s（x）的互功率谱。若假定信号s（x）与干扰n（x）彼此不相关，即Rns（x）＝0，Pns（ω）＝0，则

P_fs（ω）＝P_s（ω）＋P_ns（ω）＝P_s（ω）

P_f（ω）＝P_s（ω）＋P_n（ω）

则式（3-7-22）变为

地球物理勘探概论

因为P_s（ω）＝S^*（ω）·S（ω），P_n（ω）＝N^*（ω）·N（ω），其中*表示共轭。得：

地球物理勘探概论

这样就从一般形式的维纳滤波器得到特殊形式的维纳滤波器。

（二）分离浅源场和深源场

1.维纳滤波器

对于式（3-7-23），为了求出

和

，我们不妨假设：

地球物理勘探概论

地球物理勘探概论

即有用信号及干扰（或称区域场及局部场）分别由埋深h₁和h₂（h₁＞h₂）的地质体所引起，当地质体形态相近时：

F₁（ω）＝F₂（ω）

由此可得分离深源场（区域场）的频率响应：

地球物理勘探概论

分离浅源场（局部场）的频率响应为

地球物理勘探概论

式中的A、B、h₁、h₂四个值由实测数据的对数功率谱曲线上求得。

根据实测数据的对数功率谱曲线lnE（ω），取低频段斜率绝对值较大的直线段作为深源场的反映，并且这段直线的纵轴截距为A²，斜率一半的负数为h₁；中高频段斜率较小的直线段为浅源场的反映，并用其截距求出B²，用斜率之半的负数求h₂，如图3-7-8所示。

图3-7-8 对数功率谱曲线

2.匹配滤波器

如果把功率谱写成如下形式：

P_fs（ω）＝F^*（ω）S（ω），

P_f（ω）＝F^*（ω）F（ω），

代入（3-7-20）式便有：

地球物理勘探概论

如果令S（ω）与N（ω）相同相位（水平位置重合，深度不同的物体相位可以相同），则式（3-7-23）可以得到另哪拆一种特殊形式的滤波器：

地球物理勘探概论

即有分离深源场（区域场）的频率响应：

地球物理勘探概论

以及分离浅源场（局部场）的频率响应为

地球物理勘探概论

（三）实现方法

（1）利用傅里叶变换，由实测异常求频谱：

地球物理勘探概论

（2）由傅里叶变换的实部与虚部求对数功率谱lnE（ω）。

E（ω）＝Re²（ω）＋Im²（ω）

（3）根据对数功率谱曲线lnE（ω）-ω求h₁，h₂，B/A等参数，构制匹配滤波因子。

（4）把实测异常频谱乘以相应滤波因子，得到浅源场（或深源场）的频谱。

（5）反傅里叶变换得到分离的浅源场与深源场。

上述过程可以一次在计算机上实现：对于算出的功率谱，利用可视化技术显示对数功率谱曲线并用鼠标画出深源场与浅源场的回归直线，自动计算其斜率和纵轴截距，即可构制出匹配滤波因子，进行滤波分离不同深度的场源。

（四）应用效果分析

为了检验匹配滤波方法的有效性，设计了水平圆柱体理论模型，剖面长256个测点，点距10m，浅部场源的水平圆柱体中心埋深100m，截面有效磁矩500×10^－3A·m²，深部场源的水平圆柱体中心埋深300m，截面有效磁矩10000×10^－3A·m²，分别正演计算后再相加作为检验该方法的观测值。计算出对数功率谱后，人工用鼠标在屏幕上选择深源场（即低频段）和浅源场（即中高频段）的斜率和截距之比，得h＝75.85m，H＝147.48m，B/b＝10.0。由这些参数构制的匹配滤波器分离出的浅源场和深源场如图3-7-9所示。不难看出，匹配滤波法在一定条件下，能较好地分离出浅深不同地质体产生的场。

图3-7-9 匹配滤波法分离水平圆柱体理论模型的场

1—深、浅两个水平圆柱体的场；2—浅部水平圆柱体的场；3—深部水平圆柱体的场；4—分离后浅源场；5—分离后的深源场

㈢ 给定系统函数怎么在MATLAB中在fitter design中设计滤波器

基于MATLAB信号处理工具箱的数字滤波器设计与仿真
摘要：传统的数字滤波器的设计过程复杂，计算工作量大，滤波特性调整困难，影响了它的应用。本文介绍了一种利用MATLAB信号处理工具箱（Signal Processing Toolbox）快速有效的设计由软件组成的常规数字滤波器的设计方法。给出了使用MATLAB语言进行程序设计和利用信号处理工具箱的FDATool工具进行界面设计的详细步骤。利用MATLAB设计滤波器，可以随时对比设计要求和滤波器特性调整参数，直观简便，极大的减轻了工作量，有利于滤波器设计的最优化。本文还介绍了如何利用MATLAB环境下的仿真软件Simulink对所设计的滤波器进行模拟仿真。
关键词：数字滤波器 MATLAB FIR IIR
引言：
在电力系统微机保护和二次控制中，很多信号的处理与分析都是基于对正弦基波和某些整次谐波的分析，而系统电压电流信号（尤其是故障瞬变过程）中混有各种复杂成分，所以滤波器一直是电力系统二次装置的关键部件【1】。目前微机保护和二桥悉滑次信号处理软件主要采用数字滤波器。传统的数字滤波器设计使用繁琐的公式计算，改变参数后需要重新计算，在设计滤波器尤其是高阶滤波器时工作量很大。利用MATLAB信号处理工具箱（Signal Processing Toolbox）可以快速有效的实现数字滤波器的设计与仿真。
1 数字滤波器及传统设计方法
数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表陆迹输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。
IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配。所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。FIR数字滤波器的单位脉冲响应是有限长序列。它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。
在对滤波器实际设计时，整个过程的运算量是很大的。例如利用窗函数法【2】设计M阶FIR低通滤波器时，首先要根据（1）式计算出理想低通滤波器的单位冲激响应序列，然后根据（2）式计算出M个滤波器系数。当滤波器阶数比较高时，计算量比较大，设计过程中改变参数或滤波器类型时都要重新计算。
（1）
（2）
设计完成后对已设计的滤波器的频率响应要进行校核，要得到幅频相频响应特性，运算量也是很大的。我们平时所要设计的数字滤波器，阶数和类型并不一定是完全给定的，很多时候都是要根据设计要求和滤波效果不断的调整，以达到设计的最优化。在这种情况下，滤波器的设计就要进行大量复杂的运算，单纯的靠公式计算和编制简单的程序很难在短时间内完成设计。利用MATLAB强大的计算功能进行计算机辅助设计，可以快速有效的设计数字滤波器，大大的简化了计算量，直观简便。
2数字滤波器的MATLAB设计
2.1 FDATool界面设计
2.1.1 FDATool的介绍
FDATool（Filter Design Analysis Tool）是MATLAB信号处理工具箱里专用的滤波器设计分析工具，MATLAB6.0以上的版本还专门增加了滤波器设计工具箱（Filter Design Toolbox）。FDATool可以设计几乎所有的基本的常规滤波器，包括FIR和IIR的各种敏腊设计方法。它操作简单，方便灵活。
FDATool界面总共分两大部分，一部分是Design Filter，在界面的下半部，用来设置滤波器的设计参数，另一部分则是特性区，在界面的上半部分，用来显示滤波器的各种特性。Design Filter部分主要分为：
Filter Type（滤波器类型）选项，包括Lowpass（低通）、Highpass（高通）、Bandpass（带通）、Bandstop（带阻）和特殊的FIR滤波器。
Design Method（设计方法）选项，包括IIR滤波器的Butterworth（巴特沃思）法、Chebyshev Type I（切比雪夫I型）法、 Chebyshev Type II（切比雪夫II型） 法、Elliptic（椭圆滤波器）法和FIR滤波器的Equiripple法、Least-Squares（最小乘方）法、Window（窗函数）法。
Filter Order（滤波器阶数）选项，定义滤波器的阶数，包括Specify Order（指定阶数）和Minimum Order（最小阶数）。在Specify Order中填入所要设计的滤波器的阶数（N阶滤波器，Specify Order＝N-1），如果选择Minimum Order则MATLAB根据所选择的滤波器类型自动使用最小阶数。
Frenquency Specifications选项，可以详细定义频带的各参数，包括采样频率Fs和频带的截止频率。它的具体选项由Filter Type选项和Design Method选项决定，例如Bandpass（带通）滤波器需要定义Fstop1（下阻带截止频率）、Fpass1（通带下限截止频率）、Fpass2（通带上限截止频率）、Fstop2（上阻带截止频率），而Lowpass（低通）滤波器只需要定义Fstop1、Fpass1。采用窗函数设计滤波器时，由于过渡带是由窗函数的类型和阶数所决定的，所以只需要定义通带截止频率，而不必定义阻带参数。
Magnitude Specifications选项，可以定义幅值衰减的情况。例如设计带通滤波器时，可以定义Wstop1（频率Fstop1处的幅值衰减）、Wpass（通带范围内的幅值衰减）、Wstop2（频率Fstop2处的幅值衰减）。当采用窗函数设计时，通带截止频率处的幅值衰减固定为6db，所以不必定义。
Window Specifications选项，当选取采用窗函数设计时，该选项可定义，它包含了各种窗函数。
2.1.2 带通滤波器设计实例
本文将以一个FIR 滤波器的设计为例来说明如何使用MATLAB设计数字滤波器：在小电流接地系统中注入83.3Hz的正弦信号，对其进行跟踪分析，要求设计一带通数字滤波器，滤除工频及整次谐波，以便在非常复杂的信号中分离出该注入信号。参数要求：96阶FIR数字滤波器，采样频率1000Hz，采用Hamming窗函数设计。
本例中，首先在Filter Type中选择Bandpass（带通滤波器）；在Design Method选项中选择FIR Window（FIR滤波器窗函数法），接着在Window Specifications选项中选取Hamming；指定Filter Order项中的Specify Order＝95；由于采用窗函数法设计，只要给出通带下限截止频率Fc1和通带上限截止频率Fc2，选取Fc1＝70Hz，Fc2＝84Hz。设置完以后点击Design Filter即可得到所设计的FIR滤波器。通过菜单选项Analysis可以在特性区看到所设计滤波器的幅频响应、相频响应、零极点配置和滤波器系数等各种特性。设计完成后将结果保存为1.fda文件。
在设计过程中，可以对比滤波器幅频相频特性和设计要求，随时调整参数和滤波器类型，
以便得到最佳效果。其它类型的FIR滤波器和IIR滤波器也都可以使用FDATool来设计。
图1 滤波器幅频和相频响应（特性区）
Fig.1 Magnitude Response and Phase Response of the filter
2.2 程序设计法
在MATLAB中，对各种滤波器的设计都有相应的计算振幅响应的函数【3】，可以用来做滤波器的程序设计。
上例的带通滤波器可以用程序设计：
c=95; ％定义滤波器阶数96阶
w1=2*pi*fc1/fs;
w2=2*pi*fc2/fs; %参数转换，将模拟滤波器的技术指标转换为数字滤波器的技术指标
window=hamming(c+1); %使用hamming窗函数
h=fir1(c,[w1/pi w2/pi],window); ％使用标准响应的加窗设计函数fir1
freqz(h,1,512); ％数字滤波器频率响应
在MATLAB环境下运行该程序即可得到滤波器幅频相频响应曲线和滤波器系数h。篇幅所限，这里不再将源程序详细列出。
3 Simulink仿真
本文通过调用Simulink中的功能模块构成数字滤波器的仿真框图，在仿真过程中，可以双击各功能模块，随时改变参数，获得不同状态下的仿真结果。例如构造以基波为主的原始信号，，通过Simulink环境下的Digital Filter Design（数字滤波器设计）模块导入2.1.2中FDATool所设计的滤波器文件1.fda。仿真图和滤波效果图如图2所示。
图2 Simulink仿真图及滤波效果图
Fig.2 Simulated connections and waveform
可以看到经过离散采样、数字滤波后分离出了83.3Hz的频率分量（scope1）。之所以选取上面的叠加信号作为原始信号，是由于在实际工作中是要对已经经过差分滤波的信号进一步做带通滤波，信号的各分量基本同一致，可以反映实际的情况。本例设计的滤波器已在实际工作中应用，取得了不错的效果。
4 结论
利用MATLAB的强大运算功能，基于MATLAB信号处理工具箱（Signal Processing Toolbox）的数字滤波器设计法可以快速有效的设计由软件组成的常规数字滤波器，设计方便、快捷，极大的减轻了工作量。在设计过程中可以对比滤波器特性，随时更改参数，以达到滤波器设计的最优化。利用MATLAB设计数字滤波器在电力系统二次信号处理软件和微机保护中，有着广泛的应用前景。
参考文献
1． 陈德树. 计算机继电保护原理与技术【M】北京：水利电力出版社，1992.
2． 蒋志凯. 数字滤波与卡尔曼滤波【M】北京：中国科学技术出版社，1993
3． 楼顺天、李博菡. 基于MATLAB的系统分析与设计－信号处理【M】西安：西安电子科技大学出版社，1998.
4． 胡广书. 数字信号处理：理论、算法与实现【M】.北京：清华大学出版社，1997.
5． 蒙以正. MATLAB5.X应用与技巧【M】北京：科学出版社，1999.

㈣ 在数字滤波器中,对于变化较慢的参数如温度,应采用什么滤波方法

应采用中值滤波，中值滤波可以对某一被测参数连续采样n次（一般n应为奇数），然后将这些采样值进行排序，选取中间值为本次采样值。

中值滤波对于去掉偶然因素引起的波动或者采样器不稳定而造成的误差所引起的脉冲性干扰比较有效，如电网的波动、变送器的临时故障等。对温度、液位等缓慢变化的被测参数，采用中值滤波法一般能收到良好的滤波效果。但对流量、速度等快速变化的被测参数，一般不宜采用。

(4)最近邻滤波快速方法扩展阅读

中值滤波可以把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替，让周围的像素值接近的真实值，从而消除孤立的噪声点。

在中值滤波窗口内各点有相同的输出作用，若强调中间点或离该点较近的作用点，可改变窗口中变量个数，使多个变量值等于一个点值，再对扩展后的灰度值数字序列求中值。缺点是对边缘像素进行窗口扩展后，将超出图像边界，引起边界效应

㈤ 基于频域的滤波方法有哪些

频域滤波一般都将信号变换到频域，再同所设计的窗函数项乘，然后反变换到时域。窗函数是根据所需滤除的频率分量所决定的。数字域上常用的滤波是FIR，IRR。还有一种同态滤波是分离穗姿两樱坦个卷积分量的。
这些都是经典的频域滤波方法，还有一些现代的滤波方式，一般都是针对具体研究方向才适用，你可以参考一下这里：
http://ke..com/view/162707.htm
频域滤波的意义
1.在时域进行滤波，需要作卷积运算，而猜颂绝转化到时域做的事乘法运算，计算量小，而且有FFT和IFFT的快速算法，所以工程实现上频域滤波容易一些。
2.在频域设计滤波器，滤波器参数的物理意义很明确，分析起来很直观。

㈥ 振幅处理及提高信噪比、分辨率的处理方法

在地震资料处理中，高度保持地震波的真振幅特征，尽量提高地震记录的信噪比和分辨率，即称为“三高”处理，这一直是地震资料处理人员追求的目标。因为“三高”处理的质量直接影响到岩性参数提取以及地震勘探的精度和效果。

10.3.1 真振幅恢复

保持地震波的真振幅特征（简称保幅处理），从广义讲应包含两大方面内容：即真振幅恢复（或称振幅补偿）和其他各项处理中的振幅保持问题。本节主要讨论真振幅恢复的方法，而对其他各项处理中凡要影响到振幅特征的处理方法，则要采取相应的措施，尽可能的使振幅的相对关系保持不变。

地震记录经增益恢复处理后，其振幅特征已与地表检波器所接收到的地震波振幅特征一致。这种振幅仍不称为真振幅，我们所谓的真振幅是指由地层波阻差而产生的反射波振幅，即能反映地层岩性变化的振幅。在地表所接收到的振幅除有地层波阻抗的变化因素外，还有球面扩散因素以及非弹性衰减的因素，因此需要消除球面扩散和非弹性衰减的影响，恢复地震波的真振幅特征。

球面扩散是当波离开震源传播时由于波前扩展造成的振幅衰减。这样的振幅衰减（A）与传播距离r成反比

勘查技术工程学

其中v是界面上覆介质的平均速度；t是反射的记录时间。对球面扩散作校正需要用时变函数vt乘以数据。

非弹性衰减是弹性波能量在岩石中传播时，由于内摩擦而耗散为热被地层吸收的结果。原理部分已说明这种衰减是频率和传播距离的指数形式的函数

勘查技术工程学

其中α为非弹性衰减系数（吸收系数）

勘查技术工程学

所以，用函数e^αvt乘以数据就可校正非弹性衰减。至此，真振幅恢复处理完成。

系数α可从增益恢复及球面扩散校正后的振幅-时间函数来测定。为了得到α的较好统计估计，要用一组地震道测定能量来求得衰减曲线。

还有另一种真振幅恢复的方法，这时不需要速度信息。在增益恢复之后，假设振幅衰减是指数函数。因此，按照最小平方法，用指数函数拟合增益校正后的记录，就得到真振幅校正函数（即包括球面扩散和非弹性衰减校正两者）。

前述已知，波前发散因子K为

勘查技术工程学

式中r和t分别为波的传播距离和传播时间，C和a为与地层速度有关的常数。

吸收衰减因子是

勘查技术工程学

式中α是吸收系数；b是待定的常数。波前发散和吸收衰减总的影响是

勘查技术工程学

求取a和b的方法如下。

从地震记录上读取反射波的振幅极值（波峰或波谷），以（10.3-5）为回归方程，得

勘查技术工程学

式中：u_t=lnA_i-lnt_i；A_i，t_i为振幅极值及其对应的时间；N为振幅极值的点数。校正函数是a^-1te^bt。

为了获得有代表性的真振幅恢复参数，所选的地震道应是没有多次波及有较高的信噪比。对地质条件稳定地区，一组参数就可代表全区。在工区内地质条件有较大变化时，这些参数要重新计算。

10.3.2 提高信噪比的数字滤波处理

在地震勘探中，用于解决地质任务的地震波称为有效波，而其他波统称为干扰波。压制干扰，提高信噪比是一项贯穿地震勘探全过程的任务。除在野外数据采集中采取相应措施压制干扰外，在地震资料数字处理中数字滤波也是一项非常重要的提高信噪比的措施。

数字滤波方法是利用有效波和干扰之间频率和视速度方面的差异来压制干扰的，分别称为频率滤波和视速度滤波。又因频率滤波只需对单道数据进行运算，故称为一维频率滤波。实现视速度滤波需同时处理多道数据，故称为二维视速度滤波。本节主要介绍这两种滤波方法。

10.3.2.1 一维频率滤波

所谓一维数字滤波是指用计算机实现对单变量信号的滤波，该单变量可以是时间或频率，也可以是空间或波数。以时间或频率为例讨论一维数字滤波，其他原理相同。

（1）一维数字滤波原理

设地震记录x（t）是由有效波s（t）和干扰波n（t）组成，即

勘查技术工程学

其频谱为

勘查技术工程学

式中：X（f）为 x（t）的频谱；S（f）、N（f）分别为 s（t）、n（t）的频谱。如果 X（f）的振幅谱|X（f）|可用图10-6表示。说明有效波的振幅谱|S（f）|处在低频段，而干扰波的振幅谱处于高频段。

图10-6 有效波和干扰波频谱分布示意图

若设计一频率域函数 H（f）的振幅谱为|H（f）|，

勘查技术工程学

其图形为图10-7（a）所示。

令

勘查技术工程学

及

勘查技术工程学

在时间域有（利用傅里叶变换的褶积定理）

勘查技术工程学

称 H（f）为一维滤波器频率响应，（10.3-9）式为频率域滤波方程，h（t）为 H（f）的时间域函数，称为一维滤波器滤波因子（图10-7（b））。（10.3-11）为时间域滤波方程，y（t）和 Y（f）分别为滤波后仅存在有效波的地震记录及频谱，φ_x（f）、φ_y（f）、φ_h（f）分别为滤波前、滤波后地震记录及滤波器的相位谱，以上滤波主要是利用了有效波和干扰波的频率差异消除干扰波，故也称为频率滤波。

图10-7 滤波频率响应及滤波因子

以上所述的滤波器称为理想低通滤波，根据有效波和干扰波的频段分布不同，还可将滤波器分为理想带通滤波器、理想高通滤波器等。所谓理想是指滤波器的频率响应是一个矩形门，门内的有效波无畸变地通过，称为通频带，而门外的干扰波全部消除。在数字滤波中这一点实际是做不到的。因为数字滤波时所能处理的滤波因子只能是有限长，而由间断函数组成的理想滤波器的滤波因子是无限长的。实际应用中只能截断为有限长，截断后就会出现截断效应，即截断后的滤波因子所对应的频率响应不再是一个理想的矩形门，而是一条接近矩形门，但有振幅波动的曲线，这种现象称为吉普斯现象。

由于频率响应曲线在通频带内是波动的曲线，滤波后有效波必定会发生畸变。另外，在通频带外也是波动的曲线，必定不能有效地压制干扰。为了避免吉普斯现象，可采用若干方法，其中之一是镶边法。

10.3.2.2 二维视速度滤波

（1）二维视速度滤波的提出

在地震勘探中，有时有效波和干扰波的频谱成分十分接近甚至重合，这时无法利用频率滤波压制干扰，需要利用有效波和干扰波在其他方面的差异来进行滤波。如果有效波和干扰波在视速度分布方面有差异，则可进行视速度滤波。这种滤波要同时对若干道进行计算才能得到输出，因此是一种二维滤波。

地表接收的地震波动实际上是时间和空间的二维函数g（t，x），即是振动图和波剖面的组合，二者之间通过

勘查技术工程学

发生内在联系。式中k为空间波数，表示单位长度上波长的个数，f为频率，描述单位时间内振动次数，v为波速。

实际地震勘探总是沿地面测线进行观测，上述波数和速度应以波数分量k_x和视速度v^*代入。则有

勘查技术工程学

既然地震波动是空间变量x和时间变量t的二维函数，且空间和时间存在着密切关系，无论单独进行哪一维滤波都会引起另一维特性的变化（例如单独进行频率滤波会改变波剖面形状，单独进行波数滤波会影响振动图形，产生频率畸变），产生不良效果。那么只有根据二者的内在联系组成时间空间域（或频率波数域）滤波，才能达到压制干扰，突出有效波的目的。因此，应该进行二维滤波。

（2）二维视速度滤波的原理

二维滤波原理是建立在二维傅里叶变换基础上的。沿地面直测线观测到的地震波动g（t，x）是一个随时间和空间变化的波，通过二维正、反傅里叶变换得到其频率波数谱G（ω，k_x）和时空函数。

勘查技术工程学

上式说明，g（t，x）是由无数个角频率为ω=2πf、波数为k_x的平面简谐波所组成，它们沿测线以视速度v^*传播。

如果有效波和干扰波的平面简谐波成分有差异，有效波的平面谐波成分以与干扰波的平面谐波成分不同的视速度传播如图10-8，则可用二维视速度滤波将它们分开，达到压制干扰，提高信噪比的目的。

（3）二维滤波的计算

图10-8 有效波和干扰波以不同成分平面简谐波的传播

二维线性滤波器的性质由其空间-时间特性h（t，x）或频率-波数特性H（ω，k_x）所确定。同一维滤波一样，在时-空域中，二维滤波由输入信号g（t，x）与滤波

算子h（t，x）的二维褶积运算实现，在频率-波数域中，由输入信号的谱G（ω，k_x）与滤波器的频率波数特性H（ω，k_x）相乘来完成。

勘查技术工程学

由于地震观测的离散性和排列长度的有限性，必须用有限个（N个）记录道的求和来代替对空间坐标的积分。

勘查技术工程学

式中，n为原始道号；m为结果道号。

由（10.3-15）式可见，二维褶积可归结为对一维褶积的结果再求和。故测线上任一点处二维滤波的结果可由N个地震道的一维滤波结果相加得到。这时每一道用各自的滤波器处理，其时间特性h_m-n（t）取决于该道与输出道之间的距离。沿测线依次计算，可以得到全测线上的二维滤波结果（图10-9）。

与理想一维滤波一样，理想二维滤波也要求在通放带内频率-波数响应的振幅谱为1，在通放带外为0，相位谱亦为0，即零相位滤波。因此，二维理想滤波器的频率-波数响应是正实对称函数（二维对称，即对两个参量均对称），空间时间因子必为实对称函数。二维滤波同样存在伪门现象和吉普斯现象，也可采用镶边法和乘因子法解决。因是二维函数，情况复杂得多，通常只采用减小采样间隔（包括时间采样间隔Δt 和频率采样间隔Δf）和增大计算点数（包括时、空二方向上的点数 M 和N）的方法。

图10-9 二维滤波计算示意图（N=5）

（4）扇形滤波

最常用的二维滤波是扇形滤波。它能滤去低视速度和高频的干扰。其频率波数响应为

勘查技术工程学

图10-10 扇形滤波器的频率波数响应

通放带在f-k_x平面上构成由坐标原点出发，以f轴和k_x轴为对称的扇形区域（图10-10）。因此这种滤波器称为扇形滤波器。

利用傅里叶反变换可求出其因子为

勘查技术工程学

当在计算机上实现运算时，需要离散化。对时间采样：t=nΔ，n=0，±1，±2，……，Δ为时间采样间隔，Δ=1／2f_c。空间采样间隔即输入道的道间距Δx。

由标准扇形滤波器可以组构出既压制高视速度干扰，又压制低视速度干扰的切饼式滤波器，进而还可组构出同时压制高、低频干扰的带通扇形滤波器和带通切饼式滤波器。

在叠加前应用扇形滤波，压制的目标可以是面波、散射波、折射波或电缆振动产生的波。至于在叠加后的应用，则可压制从倾斜界面上产生的多次反射或侧面波。

10.3.3 提高纵向分辨率的反滤波处理

由地震波的传播理论可知，在介质中地震波是以地震子波的形式在地下传播。地面接收到的反射波地震记录是地层反射系数与地震子波的褶积。因此，地层相当一个滤波器，使反射系数序列变成了由子波组成的地震记录，降低了地震勘探的纵向分辨率。反滤波的目的就是要设计一个反滤波器，再对地震记录滤波，消除地层滤波的作用，提高地震记录的纵向分辨率。

由前所述，地震记录是地层反射系数序列r（t）与地震子波b（t）的褶积，即

勘查技术工程学

由于子波的问题，使高分辨率的反射系数脉冲序列变成了低分辨的地震记录，b（t）就相当地层滤波因子。为提高分辨率，可设计一个反滤波器，设反滤波因子为a（t），并要求a（t）与b（t）满足以下关系

勘查技术工程学

用a（t）对地震记录x（t）反滤波

勘查技术工程学

其结果为反射系数序列。以上即为反滤波的基本原理。

反滤波在具体实现时，核心是确定反滤波因子a（t）。由于地震子波的不确定性以及地震记录中噪音干扰的存在，实际中要确定精确的a（t）是非常困难的，甚至是不可能的。为此在不同的近似假设条件下，相继研究了很多种确定反滤波因子a（t）的方法，这些方法基本可以分为两大类：一类是先求取地震子波b（t），再根据b（t）求a（t）；另一大类是直接从地震记录中求a（t）。每一类中又有很多不同的方法（就仅反滤波方法之多，说明了反滤波处理的难度）。下面就反滤波方法中具有代表性的几种反滤波进行讨论。

10.3.3.1 地层反滤波

地层反滤波属于先求子波b（t），再求a（t）的方法。该方法要求有测井资料以及较好的井旁地震记录道。首先由声波测井资料转换与井旁地震记录道x（t）相匹配的地层反射系数序列r（t），对r（t）及x（t）求其频谱可得频率域方程为

勘查技术工程学

即有

勘查技术工程学

式中B（ω）为子波b（t）的谱，再由子波与反滤波因子的关系有

勘查技术工程学

经反傅里叶变换得 a（t）。式中 A（ω）为反滤波因子的频谱。写成 z 变换，为 A（z）=，可见A（z）是一个有理分式，要使A（z）具有稳定性，分母多项式B（z）的根必须在单位圆外，即要求子波b（t）为最小相位。

利用测井和井旁地震道求取子波及反滤波因子，即可用该反滤波因子对测线的其他道进行反滤波。

10.3.3.2 最小平方反滤波

最小平方反滤波是最小平方滤波（或称维纳滤波、最佳滤波）在反滤波领域中的应用。

最小平方反滤波的基本思想在于设计一个滤波算子，用它把已知的输入信号转换为与给定的期望输出信号在最小平方误差的意义下是最佳接近的输出。

设输入信号为x（t），它与待求的滤波因子h（t）相褶积得到实际输出y（t），即y（t）=x（t）*h（t）。由于种种原因，实际输出y（t）不可能与预先给定的期望出（t）完全一样，只能要求二者最佳地接近。判断是否最佳接近的标准很多，最小平方误差准则是其中之一，即当二者的误差平方和为最小时，则意味着二者为最佳地接近。在这个意义下求出滤波因子h（t）所进行的滤波即为最小平方滤波。

若待求的滤波因子是反滤波因子a（t），对输入子波b（t）反滤波后的期望输出为d（t），实际输出为y（t），按最小平方原理，使二者的误差平方和为最小时求得的反滤波因子称为最小平方反滤波因子。用它对地震记录x（t）进行的反滤波为最小平方反滤波。

设输入离散信号为地震子波b（n）=｛b（0），b（1），…，b（m）｝，待求的反滤波因子a（n）=｛a（m₀），a（m₀+1），a（m₀+2），……，a（m₀+m）｝，m₀为a（t）的起始时间，（m+1）为a（t）的延续长度，b（n）与a（n）的褶积为实际输出y（n），即

勘查技术工程学

为地震子波与期望输出的互相关函数。

根据最小平方原理，经推导即可得到最小平方反滤波的基本方程：

勘查技术工程学

式中，

勘查技术工程学

为地震子波的自相关函数，

勘查技术工程学

为地震子波与期望输出的互相关函数。

（10.3-24）式是一个线性方程组，写成矩阵形式为

勘查技术工程学

式中利用了自相关函数的对称性。该方程中，系数矩阵为一种特殊的正定矩阵，称为一般的托布里兹矩阵，该矩阵方程可用莱文森递推算法快速求解。

式（10.3-27）适应子波b（n）为最小相位、最大相位和混合相位。式中反滤波因子a（n）的起始时间m₀与子波的相位有关，其取值规则由子波及反滤波因子的z变换确定。

10.3.3.3 预测反滤波

预测问题是对某一物理量的未来值进行估计，利用已知的该物理量的过去值和现在值得到它在未来某一时刻的估计值（预测值）的问题。它是科学技术中十分重要的问题。天气预报、地震预报、反导弹的自动跟踪等都属于这类问题。预测实质上也是一种滤波，称为预测滤波。

（1）预测反滤波原理

根据预测理论，若将地震记录x（t）看成一个平稳的时间序列，地震子波b（t）为物理可实现的最小相位信号，反射系数r（t）为互不相关的白噪声，由地震记录的褶积模型，在（t+α）时的地震记录x（t+α）为

勘查技术工程学

分析（10.3-28）式的第一项

勘查技术工程学

可见这一项是由反射系数r（t）的将来值决定的。若令第二项为

勘查技术工程学

^x（t+α）是 t 和t 以前时刻的r（t）值决定的，也就是说（t+α）可由现在和过去的资料预测，称（t+α）为预测值。求 x（t+α）与（t+α）的差值为

勘查技术工程学

ε（t+α）称为预测误差，或称为新记录。比较（10.3-28）及（10.3-29）两式，当预测值已知时，从原记录x（t+α）中减去预测值（t+α）后形成的新记录ε（t+α）中比原记录中涉及的反射系数少，与子波褶积后波形的干涉程度轻，波形易分辨，即分辨率提高了。

在上式中α称为预测距或预测步长。当α=1时，

勘查技术工程学

即有

勘查技术工程学

这时（t+1）时刻的预测误差与反射系数之间仅差一个常数b（0）。

因此，选预测距α=1，预测误差为反射系数，达到了反滤波的目的，此时称为预测反滤波。

当α>1时，预测误差为预测滤波结果，预测滤波主要用于消除多次波，尤其是消除海上鸣震。

（2）计算预测值（t+α）的方法

在预测滤波及预测反滤波中，关键是计算预测值（t+α），其方法如下。

由反滤波方程

勘查技术工程学

代入预测值（t+α）的表达式

勘查技术工程学

式中令τ=s-j，c（s）=b（j+α）a（s-j）称为预测因子。a（t）为反滤波因子。预测值（t+α）为预测因子 c（s）与地震记录的褶积。

现在需设计一个最佳预测因子c（s），使求取的预测值（t+α）与x（t+α）最接近，即使预测误差的平方和（误差能量）

勘查技术工程学

为最小。根据最小平方原理，可得线性方程组

勘查技术工程学

式中R_xx（τ）为地震记录的自相关函数

勘查技术工程学

T为相关时窗长度，m+1是预测因子长度。将（10.3-34）写成矩阵形式为

勘查技术工程学

解此方程组即可求得预测滤波因子c（t），用它对地震记录x（t）褶积可以求出未来时刻（t+α）时的最佳预测值（t+α）。

㈦ 如何进行软件滤波

1、限幅滤波法（又称程序判断滤波法）
A、方法：
根据经验判断，确定两次采样允许的最大偏差值（设为A），每次检测到新值时判断：
如果本次值与上次值之差<=A,则本次值有效
如果本次值与上次值之差>A,则本次值无效,放弃本次值,用上次值代替本次值
B、优点：
能有效克服因偶然因素引起的脉冲干扰
C、缺点
无法抑制那种周期性的干扰
平滑度差

2、中位值滤波法
A、方法：
连续采样N次（N取奇数）
把N次采样值按大小排列
取中间值为本次有效值
B、优点：
能有效克服因偶然因素引起的波动干扰
对温度、液位变化缓慢的被测参数有良好的滤波效果
C、缺点：
对流量、速度等快速变化的参数不宜

3、算术平均滤波法
A、方法：
连续取N个采样值进行算术平均运算
N值较大时：信号平滑度较高，但灵敏度较低
N值较小时：信号平滑度较低，但灵敏度较高
N值的选取：一般流量，N=12；压力：N=4
B、优点：
适用于对一般具有随机干扰的信号进行滤波
这样信号的特点是有一个平均值，信号在某一数值范围附近上下波动
C、缺点：
对于测量速度较慢或要求数据计算速度较快的首坦谨实时控制不适用
比较浪费RAM

4、递推平均滤波法（又称滑动平均滤波法）
A、方者基法：
把连续取得的N个采样值看成一个队列
队列的长度固定为N
每次采样到一个新数据放入队尾,并扔掉原来队首的一次数据.(先进信陪先出原则)
把队列中的N个数据进行算术平均运算,就可获得新的滤波结果
N值的选取：流量，N=12；压力：N=4；液面，N=4~12；温度，N=1~4
B、优点：
对周期性干扰有良好的抑制作用，平滑度高
适用于高频振荡的系统
C、缺点：
灵敏度低
对偶然出现的脉冲性干扰的抑制作用较差
不易消除由于脉冲干扰所引起的采样值偏差
不适用于脉冲干扰比较严重的场合
比较浪费RAM

5、中位值平均滤波法（又称防脉冲干扰平均滤波法）
A、方法：
相当于“中位值滤波法”+“算术平均滤波法”
连续采样N个数据，去掉一个最大值和一个最小值
然后计算N-2个数据的算术平均值
N值的选取：3~14
B、优点：
融合了两种滤波法的优点
对于偶然出现的脉冲性干扰，可消除由于脉冲干扰所引起的采样值偏差
C、缺点：
测量速度较慢，和算术平均滤波法一样
比较浪费RAM

㈧ 多维数字信号处理滤波方法有哪些

数字信号处理是把信号用数字或符号表示成序列，通过计算机或通用（专用）信号处理设备，用数值计算方法进行各种处理，达到提取有用信息便于应用的目的。例如：滤波、检测、变换、增强、估计、识别、参数提取、频谱分析等。
一般地讲，数字信号处理涉及三个步骤：⑴模数转换(A/D转换)：把模拟信号变成数字信号，是一个对自变量和幅值同时进行离散化的过程，基本的理论保证是采样定理。⑵数字信号处理(DSP)：包括变换域分析(如频域变换)、数字滤波、识别、合成等。⑶数模转换(D/A转换)：把经过处理的数字信号还原为模拟信号。通常，这一步并不是必须的。 作为DSP的成功例子有很多，如医用CT断层成像扫描仪的发明。它是利用生物体的各个部位对X射线吸收率不同的现象，并利用各个方向扫描的投影数据再构造出检测体剖面图的仪器。这种仪器中fft(快速傅里叶变换)起到了快速计算的作用。以后相继研制出的还有：采用正电子的CT机和基于核磁共振的CT机等仪器，它们为医学领域作出了很大的贡献。
信号处理的目的是：削弱信号中的多余内容；滤出混杂的噪声和干扰；或者将信号变换成容易处理、传输、分析与识别的形式，以便后续的其它处理。 下面的示意图说明了信号处理的概念。

㈨ 频率域重磁异常滤波

在频率域中浅部地质体引起的异常较深部地质体引起的异常要尖锐得多。一个尖锐的异常其幅值从异常中心向外快速下降，以具有较强的高频成分为特征；而宽缓异常从中心向外是缓慢的衰减，具有集中于低频端的谱。异常频谱特征的这种差异，提供了分离浅部场和深部场的可能性。频率域重磁异常滤波的作用在于利用重磁异常的频谱特征，区分区域异常与局部异常，即分离叠加异常。目前有许多种滤波方法，这里介绍几种常用的滤波方法。

（一）维纳滤波与匹配滤波

维纳滤波与匹配滤波方法是根据不同埋深的场源在对数功率谱上的不同特征来构制滤波器，通过滤波达到分离区域场与局部场的目的。

1.维纳滤波器

最小均方差滤波也称为维纳滤波，其基本思路是：设计一个滤波器，使其输出与希望输出之间的均方差为最小。

设I（x）为输入信号，F（x）为I（x）滤波后的输出函数，E（x）为期望输出，ε为总均方误差，则有

勘探重力学与地磁学

相应的滤波器为最小均方差滤波器。

具体应用时常做一些假定，如设有用信号是由深部场源引起，干扰信号是由浅部场源引起。并设S_I（ω）为I（x）的频谱，S_E（ω）为E（x）的频谱，在频率域中根据巴什瓦等式可将上式写成

勘探重力学与地磁学

设深部场源、浅部场源的埋深分别为h₁和h₂（h₁＞h₂)，ω为角频率，将整个场向下延拓h₂，使浅源成为干扰。异常其振幅谱为B，令深源异常的频谱服从下列关系：

勘探重力学与地磁学

则求ε极小公式可改写为

勘探重力学与地磁学

为求频率响应ϕ（ω），用［ϕ（ω）＋λG（ω）］代替ϕ（链衫ω），代入上式。这里λ为任意数，G（ω）为ω的函数，棚孙腔其性质同ϕ（ω），求使

时ϕ（ω）值，即

勘探重力学与地磁学

由于G（ω）的任意性，最后得到

勘探重力学与地磁学

故有

勘探重力学与地磁学

式中：

。

对（10-184）式积分，选择合适的脉冲响应使ε最小（即求积分式右边的泛函的极值），经过变换得到Wiener-Hopf积分方程。对其作傅里叶变换，并假定信号s（x）与干扰n（x）彼此不相干。经推导可以得到特殊形式的维纳滤波器：

勘探重力学与地磁学

为了求出|S（ω）和|N（ω）|，不妨假设

勘探重力学与地磁学

勘探重力学与地磁学

即有用信号及干扰信号（或称区域场及局部场）分别由埋深为h₁和h₂（h₁＞h₂)的地质体所引起，l为浅部场源下延深度。上式中A，B为物性参数；H（ω）是维纳滤波器的频率响应；S（ω）为有用信号s（x)（即区域场）的频谱；N（ω）为干扰信号n（x)（即局部场）的频谱。

在高频端（10-187）式可近似写为

勘探重力学与地磁学

当地质体形态相近时，有

F₁（ω）＝F₂（ω）

由此可得分离区域场的频率响应

勘探重力学与地磁学

式中的A，B，h₁，h₂四个值由实测数据的对数功率谱曲线上求得。

根据实测数据的对数功率谱曲线lnE（ω)，取低频段斜率绝对值较大的直线段作为深部场源的反映，并切这段直线的纵轴截距为A²，斜率一半的负数为h₁，中高频段斜凯销率较小的直线段为浅部场源的反映，并用其截距求出B²，用斜率一半的负数求h₂。如图10-20所示。

2.匹配滤波器

如果令S（ω）与N（ω）相同相位（深度不同但水平位置重合的地质体其相位可以相同)，则由（10-185）可以得到另一种形式的特殊滤波器：

图10-20 对数功率谱曲线

勘探重力学与地磁学

即有分离场的频率响应

勘探重力学与地磁学

从（10-186）式可以看出，为实施匹配滤波必须先求得h₁、h₂和

等值。这些值可以从径向平均对数功率谱求得。从（10-186），（10-187）式出发，则应有

勘探重力学与地磁学

在实际工作中，可以用ω＝0时的E（0）值A²对E（ω）作规格化处理。若记作E_n（ω），则

勘探重力学与地磁学

当ω很大时，

，则

勘探重力学与地磁学

这表明在lnE_n（ω）～ω曲线的高频段的拟合直线斜率为-2h₂，而直线的截距为

。

而当ω很小时，由（10-191）式可得

勘探重力学与地磁学

lnE_n（ω)＝-2h₁ω （10-193）

因此在lnE_n（ω）～ω曲线的低频段的拟合直线斜率为-2h₁。同时由于ω＝0时，lnE_n（ω）＝0，即此拟合曲线与纵轴的交点相当于坐标原点。这就表明两条拟合直线与纵轴交点之间的距离就是

。由此可见，从径向平均对数功率谱曲线上可以取得为实施匹配滤波所需要的参数。应该指出，对重力异常，为确定上述参数，首先需要进行一次垂向一阶导数的运算，然后再作出径向平均对数功率谱曲线，确定上述参数。匹配滤波应用于磁异常的解释时，其应用条件为：①深浅源场既要同相位，又要不相干；②场源应近似于球体与下延很大棒状体叠加的情况（管志宁、安玉林，1984)。

3.实现步骤

（1）利用傅里叶变换，由实测异常求频谱。

（2）由傅里叶变换的实部和虚部求对数功率谱lnE（ω）。

E（ω）＝Re²（ω）＋Im²（ω）

（3）根据对数功率谱曲线lnE（ω）～ω，求h₁，h₂，B/A等参数，构制匹配滤波因子。

（4）将实测异常频谱乘以相应滤波因子，得到浅源场（或深源场）的频谱。

（5）反傅里叶变换得到分离的浅源场与深源场。

（二）消除高频干扰的正则化方法

在A.H.吉洪诺夫的《不适定问题的解法》一书，涉及了研究解决迅速形成“具有理想低通滤波特性和较强适应能力”的滤波因子问题。并将所构成的滤波因子称为正则化稳定因子。基于该书的理论，结合磁异常的特点安玉林等提出了正则化稳定因子（安玉林、管志宁，1985)。现给出一种正则化稳定因子：

勘探重力学与地磁学

式中：β≥2；f₀为要消除的高频干扰信号的最小频率（波数)，等于其最大水平尺度的倒数

；

；λ_x，λ_y为基波波长；频率（波数）

，

，

。

上列正则化稳定因子经理论模型检验和实际资料处理，效果较好。经实践证明，正则化参数α可直接取2≤α≤3，如取

勘探重力学与地磁学

即可直接滤除高频干扰。

正则化稳定因子中，参数f₀与λ₀具有重要意义，它们表明要消除的局部磁异常的尺度。这两个参数可以直接从原始磁异常剖面图或平面等值线图上量取，这是该方法易于应用的主要原因。图10-21是原始观测场ΔZ的剖面曲线图，ΔZ曲线即包含有区域场③，又包含大小不同的局部异常①，②。如果要消除所有局部异常，则选取

。如果仅消除局部异常②，则选

，被保留的部分将作为区域磁异常。

当欲被分离的场即不是高频，也不是低频，而是中频时，则可采用带通滤波器。设计带通滤波器的方法之一，是将低通滤波器与高通滤波器串联，其频率响应等于低通与高通滤波器频率响应的乘积。

图10-21 ΔZ剖面曲线图

已知低通滤波器的频率响应为H₁（f)，则高通滤波器的频率响应H₂（f）为：H₂（f）＝1-H₁（f）。这里给出正则化带通滤波因子为

其中

勘探重力学与地磁学

式中：f₁是欲分离出的场的最小波数，近似等于最大水平尺寸λ₁的倒数1/λ₁；f₂是欲滤除的干扰场的最小波数，近似等于最大水平尺寸λ₂的倒数1/λ₂。λ₁，λ₂可以从叠加场平面等值线图上量出，一般可取α₁＝α₂＝2.5。

若取β₁＝β₂≥2，则频率特征曲线有理想带通滤波器频率特征（图10-22)。若取β₁＝β₂＜2，则f₁，f₂处频率特征曲线变缓。若β₁≠β₂，则曲线不对称。

图10-22 带通正则化滤波因子曲线图