本文目录一览：
- 1、 随机噪声的滤波
- 2、 实现方法
- 3、AvaReader如何采用滤波,过滤掉噪声
- 4、平滑滤波的滤波方法
- 5、Gabor滤波器
- 6、plc模拟量输入滤波程序和方案？
随机噪声的滤波
在观测的结果中，不可避免地会有随机噪声。产生噪声的原因有三方面，第一是仪器本身读数不稳定引起的仪器噪声；第二是观测者肉眼的分辨力不够所引起的观测误差，这种误差具有随机性质；第三是仪器所处的环境引起的干扰，这里所指的环境包括自然因素和人文因素，干扰则可分为随机干扰和非随机干扰，后者往往具有局部场的特点，而随机干扰则是这里所指的随机噪声。
不失一般性，假定观测结果是单变数x的函数，x可以是位置变量，也可以是时间变量。设观测值用t1（x）表示，区域场用t0（x）表示，噪声用n（x）表示，则有
信息论、系统论与地质找矿工作
假定n（x）具有分段均值为零的性质，即
信息论、系统论与地质找矿工作
式中x0为计算平均值的点的坐标，△x为采样点点距，（2n+1）△x为计算平均值的窗口。
选择削弱噪声滤波方法的原则是：依据区域场的性质，选择方法，使滤波所引起的区域场的畸变尽可能地小，而噪声的削弱程度则尽可能地大。因此，在（7—2）式成立条件下，根据区域场的特点，常用来削弱噪声的滤波方法有以下五种：
第一种方法，滑动窗口平均法。这种方法适用于区域场在滑动窗口内呈线性变化的情况，即
信息论、系统论与地质找矿工作
将坐标原点取在窗口中心，则有xi=i△x及
信息论、系统论与地质找矿工作
对（7—1）式在（-n,n）窗口范围内取平均，并令平均值为A，则有
信息论、系统论与地质找矿工作
信息论、系统论与地质找矿工作
根据（7—2）及（7—4）式，上式化简为：
A=a=t0（0）
即平均值即为在窗口中心的区域异常值。图7—2是一个例子。图中区域场用下式表示：
信息论、系统论与地质找矿工作
式中x=i△x；△x=1;i=1,2,3，…，100。
图7—2 叠加正负等幅噪声的信号
1—信号；2—叠加噪声后的信号；3—噪声
随机噪声为正负等幅，振幅为10，加在8到90号点上，滤波窗口取11△x。由于用3单数点求平均值，故（7—2）式不成立，而有噪声的平均值为A：
信息论、系统论与地质找矿工作
如果对平均的结果再取一次平均，设其均值为A1，则有
信息论、系统论与地质找矿工作
图7—3是对图7—2中的t（x）两次取平均的结果。从图看出，消除随机噪声的效果是很好的。但仔细看图，还可以发现平均的结果在实际值的上下有小的波动。为了图面清晰，滤波结果每隔5个点取1个点在图上表示，而区域场则用曲线表示。
第二种方法，自动趋势面拟合法。这种方法适用区域场可用变数x的低价多项式近似时的情况。例如对（7—5）式，令y=x-50，则有
t0（y）=260000/（10000+y2）=260000（10000+y2）-1
由于-50≤y≤50，故y2/10000的最大值为2500/10000=0.25。因此，在（7—5）式中x的取值范围内，令x1=y2/1000,x1＜1，故可将上式展开为
t0（y）=260000（10000+y2）-1=26（1+y2/10000）-1
信息论、系统论与地质找矿工作
已知
的最大值为0.25×0.25=0.0625，故上述展开式可取到x1的1次项，即y的平方项。
图7—3 图7—2曲线2两次滑动窗口平均的结果（窗口大小为11）
1—信号；2—平均法的结果
图7—4图 7—2曲线2用自动趋势面拟合的结果
1—信号；2—自动趋势面拟合的结果
图7—5 图7—2曲线2用高频截断法滤波结果
1—信号；2—高频截断滤波结果
因此，对图7—2中的t（x），也可用自动趋势面拟合法直接求出区域场。图7—4是用这种方法对t（x）的拟合结果，计算程序自动确定拟合多项式的阶次为2。从图看出，计算的结果是令人满意的。
第三种方法，频率域的高频截断法。由于区域异常以低频成分为主，而随机噪声则以高频为主，故高频截断能很好地消除信号中的噪声。图7—5是用这种方法对图7—2中t（x）的滤波结果。从图可以看出，滤波结果是好的。
高频截断所用的滤波算子或频率域中的响应函数L（ω）为
信息论、系统论与地质找矿工作
式中△ω＝2π/[dx×（2N+1）];dx＝取样点距；N＝采样点数的1/2（采样点数应为偶数）；k为滤波器参数，应根据t（x）的频谱特点选定。此次滤波选k=0.5N，由于采样点数为100，即N=50，故有k=0.5N=25。
第四种方法，斯特拉霍夫滤波器[32]。这种滤波器应用的条件是：噪声n（x）应是各态遍历的、平均值为零的随机函数；
，h为区域场源的埋深，t0（ω）为T0（x）的傅里叶变换。图7—6是用这个方法对图7—2中t1（x）的滤波结果。滤波所用公式是：
信息论、系统论与地质找矿工作
式中C0=0.4452,C±1=0.0944,C±2=0.0932,C±3=0.0911.
上述系数适合于△x/h=0.1的情况。由于此处区域场源的埋深h=100，即△x/h=0.01，故滤波效果不如前述三种方法的滤波效果好。
图7—6 图7—2曲线2用斯特拉霍夫滤波器的滤波结果
1—信号；2—斯特拉霍夫滤波器的滤波结果
第五种方法，向上延拓。从数学上讲，向上延拓高度h相当于对t1（x）的频谱T1（ω）乘一个因子exp（-hw）。因此，以高频成分为主的噪声衰减比以低频成分为主的区域异常衰减快得多。图7—7是用这种方法对图7—2中t1（x）的滤波结果，上延高度h=4△x。从图看出，滤波结果是好的。这是因为区域场源埋深为100△x，故上延4△x高度将使区域极大值下降3.8%，影响不大。对噪声而言，从频谱图7—8可以看出，t1（x）中噪声的谱主要集中于ω＞π/2△x部分，而当ω＝π/2△x时，exp（-4△xπ/2△x）=exp（-2π/△x）＝exp（-2π）＝1.87×10-3。即噪声的衰减将大于99.8%。
图7—7 图7—2曲线2向上延拓的结果（延拓高度为4倍点距）
1—信号；2—向上延拓高度4△x的结果（△x为采样点距）
图7—8 功率谱曲线图
因此，用高频截断法或向上延拓法削弱信号中的噪声时，最好先作频谱分析，计算出t1（x）的频谱图，了解噪声谱分布的主要区间，据此确定（7—6）式中的k值或向上延拓高度。
实现方法
4.2.2.1 M周期的谱分析
1.M周期谱分析步骤
M周期的识别与划分用谱分析来实现，具体步骤包括：
①对测井曲线滤波，目的是了解测井曲线的频谱组成和消除各种干扰。
②消除测井曲线的整体漂移，减小长周期信号的影响。
③求取频谱，可用付氏变换直接求取频谱，也可用Burg最大熵谱法求取频谱。
④求取M周期，划分沉积旋回。根据已知的经验，我们可以确定出塔北各时期的最大沉积速率。已知M周期分别为413ka、125ka、95ka和21ka，这样，就可确定各级M周期所对应的最大地层厚度，从而可确定分析时窗的大小（见图4—6），并可保证选取的时窗大于一个M周期。
图4—6 分析时窗选择图解
具体作法是：在横轴上找到已知的沉积速率，投映到待研究的周期的对应的直线上，作纵轴的垂线，读出其厚度值，即为所需窗长。所选择的窗长不应小于该厚度。假定沉积速率为每千年500mm，若要识别413ka的周期，对应的厚度是206.5m，窗长至少要达到这一厚度，而要识别125ka的周期则需要62.5m的窗长，而滤波分析时滑动窗口的步长的选择则很灵活，根据研究的详细程度选2m、5m、10m、20m等。
如果测井曲线中某个沉积旋回正好全部处于分析时窗内，则功率谱上就会显示几个不同波长的谱峰，虽然由于沉积速率不同，谱峰的波长是会变化的，但谱峰的频率的倒数（即其所代表的厚度）的比值，接近于M周期的比，即410:123:41。如果分析时窗跨越两个沉积旋回，则谱峰的数目与波长就会失去规律性。据此，我们可以识别沉积旋回和划分旋回边界。
测井曲线的谱分析，除了用于确认M周期外，还可用于沉积速率分析，计算沉积速率。
2.M周期的确认
识别M周期可用谱峰比值法或谱峰强度。
①谱峰比值：在图4—7中谐波的波长比在5241.7m深度处，第一、第三、第五个谱的频率比为0.23∶0.68∶1.8，其倒数比（厚度比）接近410∶123∶41，即接近M周期的比值。表明在以此深度点为中心的分析窗口内，包含着一个完整的M周期。
图4—7 Milankovitch周期的谱分析
图4—7中深度5293.7m处的第一、第二、第五个大的谱峰的比为0.1:0.3:0.9，其倒数也接近410∶123∶41。
图4—7中深度5189.7m处的第一、第三、第五个谱大的谱峰的比值为0.19∶0.58:1.79，倒数的比值接近M周期410:123:41。
功率谱上极大谱峰的比接近M周期的比是M周期存在着强有力的证据。如有古地磁资料、同位素测年资料、微体古生物定年资料，亦可作为另外的证据。
②谱峰的强度：图4—7中对应谱峰最强的是410ka的M周期，这与大的周期对沉积的控制作用有关，小周期叠加在大周期的变化作用之上，它们的比例分别对应其周期长度。
3.沉积速率变化的谱嶂显示
沉积速率的渐变：在谱峰上显示为峰值向右（沉积速率变小）或向左（沉积速度变大）发生小的移动。如图4—7中5245.7m至5247.7m的第一个最大峰由0.35移到0.19（周/m），并且下一深度还在向左偏移，说明沉积速度向上减小（向下增大）。据沙11井的综合柱状图知道，该井段为浅灰色粉砂质泥岩夹绿灰色泥岩，在51.2m的时窗内顶部偏泥而底部偏砂（在自然伽马曲线上的显示更明显），反应了向上沉积速度减小。
沉积速率突变：5177.7m至5179.7m上下频谱特征的突变使我们了解到沉积速率的突变。该深度上下的谱特征有明显的差别，必然是沉积速率变化的结果.在岩性柱状图上，顶部为绿灰色粉砂质泥岩与泥质粉砂岩互层，而其下有一层绿灰色长石岩屑细砂岩。反应了向上岩性突变，沉积速率突变。
5221.7m上下的谱峰突变也同样反应了沉积速率突变（见图4—7）。
4.计算沉积速率
虽然在几万年的时间里沉积速率不可能保持不变，对应的垂向周期的波长也将有变化，但是由于M周期的时间比恒定，这些周期在窗口内的比也保持恒定，窗长内存在M周期则说明此段地层是在一个完整周期内沉积下来的，这样就不难求得沉积速率。各代表点的沉积速率标在图4—7中。计算公式如下：
V（H）＝1/（F*T）
式中，V（H）表示H深度的沉积速率、F表示由频谱图上读得的周期为T的谐波的频率（单位：周/m，相当于该周期波波长的倒数）。
以下几例分析的结果说明，利用M周期求取沉积速率是有实际意义的。
在以5221.7m为中心的长度为51.2m时窗内，顶部为浅灰白色细砂岩，中间为大段的浅绿灰色泥质粉砂岩夹一层岩屑细砂岩，说明该段为一较快速的沉积段，计算得沉积速率为48.78mm/ka。
在以5229.7m为中心的长度为51.2m的时窗内，是大段的绿灰色泥岩与浅灰色泥质粉砂岩，泥岩发育很好，表示了一种慢速的沉积过程，计算得沉积速率为7.31mm/ka。从5229.7m到5221.7m岩性变粗，由绿灰色泥岩过渡到粉砂质泥岩夹细砂岩，是向上水深变浅、粒度变粗的过程，沉积速率也由小变大。
在以5241.7m为中心的长度为51.2m时窗内，为绿灰色泥岩与粉砂质泥岩，沉积仍为一慢速过程，计算得沉积速率为10.6mm/ka。
在以5247.7m为中心的长度为51.2m时窗内，仍为绿灰色泥岩与粉砂质泥岩，5250m上下自然伽马对应两个降低段，是粉砂质泥岩含砂较多的表现，在该时窗及上下时窗范围内含粉砂较多，计算得沉积速率为16.3mm/ka，大于上下时窗内的沉积速率。
以5255.7m为中心的长度为51.2m时窗内，为泥岩与粉砂质泥岩，沉积较慢，计算得的沉积速率为13.6mm/ka。
以5267.7m为中心的长度为51.2m时窗内，为大段的泥岩与粉砂质泥岩，计算得沉积速率为16.6mm/ka。
Milankovitch周期是地层中普遍存在的，并且是在测井分辨率意义下可体现的周期信号，是利用测井分析地层周期性、进行层序划分与可容纳空间分析的基础。
4.2.2.2 Fischer图解应用条件修正
Fischer图是根据潮坪浅水环境中形成的碳酸盐岩层序来求取海平面变化的一种图解。其理论依据是：在潮坪浅水环境中，碳酸盐生长速率很高，在任何时候，由相对海平面上升所增加的可容纳空间，随时都能被新产生的碳酸盐物质所充填。因此，在相同周期的海平面变化过程中所形成的碳酸盐岩旋回的厚度，间接地反映了该周期相对海平面升降的幅度。这一基本假设前提在硅质碎屑环境中则不满足，在深水中，往往沉积地层厚度小，浅水环境中沉积地层厚度大，这就需要作水深校正，使该旋回的厚度经校正后能反映可容纳空间的大小。
根据统计分析，自然伽马、自然电位和视电阻率等曲线与岩石平均粒度有较好的对应关系（如图4—8,4—9，），黄智辉等人在不同地区的研究也证实了这种规律，如图4—10所示。对应的公式如下：
ΔGR＝（GR-GRmin）/（GRmax-GRmin）
ΔSP＝（SP-SPmin）/（SPmax-SPmin）
MD=Aexp（-CΔGR）
式中，GR、SP分别为自然伽马值和自然电位值，MD表示平均粒度，A.C为常数。
图4—8 自然伽马（GR）与岩性的对应关系（据沙32井）
在一些特别的情况下可能有不同的例子，但从总体来说，自然伽马电位与平均粒度有较好的对应关系。
在静水条件下，随可容纳空间的增大，水深增加，距离岸线也将更远，沉积物的供应必将减少，在相同的时间段内沉积的地层厚度也将减小，粒度将变细，自然伽马曲线值将增大。可容纳空间越大，水深越深，离岸线越远，沉积物的旋回厚度就越小。反之，旋回厚度就变大，而粒度就将变粗，自然伽马曲线也将减小。
图4—9 泥质含量与地层电阻率交会图
图4—10 自然伽马曲线与平均粒度对应关系（据黄智辉，1986）
这样，当一个旋回的可容纳空间大时，应对应较大的沉积物厚度，而实际情况是旋回厚度较小。而可容纳空间变小时，应对应的沉积厚度较小，而实际情况是旋回厚度大。大的沉积可容纳空间对应细粒的沉积物和较小的旋回厚度，而小的沉积可容纳空间对应粗粒的沉积物和较厚的旋回厚度。我们设计了一个指数模型来逼近这种变化，因为指数模型能展开为多项式形式，是众多的逼近函数的首选函数。校正的公式如下：
H′＝Hexp（a·Δ）
式中H′是校正后的旋回厚度，H是校正前的旋回厚度，a是校正的常数，△可表示为：
Δ＝MIN（ΔGR，ΔSP）
式中MIN是取极小的函数，由于各种因素对测井曲线的影响使得相应的测井曲线值升高，取极小可以减小偶发因素的影响。
参数a的选择要谨慎，要依据井区的环境资料或根据经验，目的是尽可能地补偿水深的影响，而又避免引进不合理的变化因素。
图4—11是校正前后的图解曲线，由于旋回多，仅作出包络线，左边第一栏是自然伽马曲线，第二栏是校正前的可容纳空间变化曲线，第三栏是校正后的可容纳空间变化曲线，校正后对应于较好的泥岩段可容纳空间变化增大更明显，更能反映相对海平面的变化，而未校正的曲线中间部分的下降显然不如校正后中—上部的一致偏高更接近于实际。
图4—11 校正前后的可容纳空间变化曲线
4.2.2.3 区域场与局部场的分离
上述处理得到的曲线含有各种变化成份，而且长波信号的幅度远大于短波信号，变化快的低幅信号显示不清或根本显示不出来。从图4—11中也可以清楚地看到这一点。为了提高分析信号的能力，采用了提取区域场与局部场的方法。
设待处理信号为Y（X），有下列实现模型：
新疆塔里木盆地北部应用层序地层学
式中，Ai是系数，N是多项式的次数，已有的序列作为它的解，用最小二乘法求解超定方程，得到系数，就可用上述公式计算区域场和局部场。
Z（X）＝Y′（X）-Y（X）
式中Z（X）表示X点的局部场值，Y′（X）表示X点的实际信号值，Y（X）是X点区域场值。
图4—12是经过各阶多项式提取后的区域场与局部场，右边第一栏为原始曲线，左边一、二、三、四栏分别是二、三、四、五次多项式代表的区域场与局部场。从图中可以清楚地看到，各短波信号明白无遗。而未处理前的信号由于长波信号太强使短波显示不出来。这种方法实现了在同一个曲线上显示各信号的能力。
图4—12 区域场与局部分离
AvaReader如何采用滤波,过滤掉噪声
1、给定均值滤波窗口长度，对窗口内数据求均值。
2、作为窗口中心点的数据的值，之后窗口向后滑动1，相邻窗口之间有重叠。
3、边界值不做处理，即两端wid_length2长度的数据使用原始数据即可过滤掉噪声。
平滑滤波的滤波方法
图像的噪声滤波器有很多种，常用的有线性滤波器，非线性滤波器。采用线性滤波如邻域平滑滤波，对受到噪声污染而退化的图像复原，在很多情况下是有效的。但大多数线性滤波器具有低通特性，去除噪声的同时也使图像的边缘变模糊了。而另一种非线性滤波器如中值滤波，在一定程度上可以克服线性滤波器所带来的图像模糊问题，在滤除噪声的同时，较好地保留了图像的边缘信息。
邻域平滑滤波原理
邻域平均法[2]是一种利用Box模版对图像进行模版操作（卷积运算）的图像平滑方法，所谓Box模版是指模版中所有系数都取相同值的模版，常用的3×3和5×5模版如下：
邻域平均法的数学含义是：
（式4-1）
式中：x,y=0,1,…,N-1;S是以（x,y）为中心的邻域的集合，M是S内的点数。
邻域平均法的思想是通过一点和邻域内像素点求平均来去除突变的像素点，从而滤掉一定噪声，其优点是算法简单，计算速度快，其代价会造成图像在一定程度上的模糊。
中值滤波原理
中值滤波[2]就是用一个奇数点的移动窗口，将窗口的中心点的值用窗口内的各点中值代替。假设窗口内有五点，其值为80、90、200、110和120，那么此窗口内各点的中值及为110。
设有一个一维序列f1,f2,…,fn,取窗口长度（点数）为m（m为奇数）,对其进行中值滤波，就是从输入序列中相继抽出m个数fi-v,…,fi-1,fi,fi+1,…,fi+v(其中fi为窗口中心值,v=(m-1)/2),再将这m个点按其数值大小顺序排序，取其序号的中心点的那个数作为滤波输出。数学公式表示为：
Yi=Med{fi-v,…,fi-1,fi,fi+1,…,fi+v} i∈N v=(m-1)/2 （式4-2）
Yi称为序列fi-v,…,fi-1,fi,fi+1,…,fi+v的中值
例如，有一序列{0，3，4，0，7}，重新排序后为{0，0，3，4，7}则Med{0，0，3，4，7}=3。此列若用平滑滤波，窗口也取5，那么平滑滤波输出为（0+3+4+0+7）/5=2.8。
把一个点的特定长度或形状的邻域称作窗口。在一维情况下，中值滤波器是一个含有奇数个像素的滑动窗口。中值滤波很容易推广到二维，此时可以利用二维形式的窗口。
对于平面图像采用的二维中值滤波可以由下式表示：
（式4-3）
式中：A为窗口，{fij}为二维数据序列,即数字图像各点的灰度值。
对于本系统，由于采集到的是24位真彩色图像，每个像素点分别有R、G、B三个灰度分量，故要在窗口内分别找到这三个分量的中值，分别用这三个中值去代替窗口中心像素点的R、G、B三个灰度分量的值。
Gabor滤波器
Gabor函数是一个用于边缘提取的线性滤波器，十分适合纹理表达和分离，其频率和方向表达同人类视觉系统类似。另外从生物学实验发现，Gabor滤波器而已很好地近似单细胞的感受野函数。
另外，Gabor函数分为实部和虚部，用实部进行滤波后图像会平滑；虚部滤波后可用来检测边缘（也许这是一个传说），不过的确可以实践一下。
Gabor变换时基于傅里叶变换的局限性提出来的，见下图：
其较为详细的说明可参见博文
为解决傅里叶变换中存在的局部性问题，其变换的基本思想是：把信号划分成许多小的时间间隔，用傅里叶变换分析每一个时间间隔中的频域信息，其处理方法是对傅里叶变换函数设置一个滑动窗口，再对其进行傅里叶变换。其表达如下图所示：
原傅里叶变换：
比较可见，Gabor变换通过增加一个窗函数来实现了局部化的操作，充分体现了时域与频域的局部化信息，同时，能够在整体上提供信号的全部信息。
plc模拟量输入滤波程序和方案？
1，硬件配置滤波， 如果是200PLC打开系统块，再Analog里设定滤波时间和频率 如果是300400PLC打开硬件配置，再相关模块里设定滤波时间和频率，这个一般是过滤高频的杂波 2，然后再程序里，编程实现： 均值滤波：我一般用最后五次采样的平均值，采样时间间隔和几次求平均值可以自己定。 中值滤波：我没用过，可以尝试。 峰值滤波：直接取多次采样的最高或最低值，也是特殊情况有用的。 总结：你首先要观察你的测量量的特性，否则滤波是低效、盲目的。
