提高信噪比的数字滤波技术及其VC++实现.docx
《提高信噪比的数字滤波技术及其VC++实现.docx》由会员分享，可在线阅读，更多相关《提高信噪比的数字滤波技术及其VC++实现.docx（14页珍藏版）》请在冰豆网上搜索。
提高信噪比的数字滤波技术及其VC++实现
提高信噪比的数字滤波技术及其VC++实现
论文提要
滤波是地震资料数据处理中最常用的一种方法。
本文首先介绍了地震数据处理中的数字滤波技术的基本原理，然后利用VC++程序设计工具进行滤波程序的编译和调试，在此基础上，选择不同的滤波参数，对滤波结果进行了对比和分析。
在软件设计的整个过程中，融入了软件工程的基本思想。
所编译的程序能够在一定程度上对地震数据进行简单快速的滤波处理。
并具有一定的滤波效果，但由于一维滤波本身的局限性，滤波虽然有一定效果，但存在不足，需要发展更加优越的滤波技术，如二维视速度滤波、非线性滤波、最小平方滤波等。
正文
一、前言
（一）选题依据
随着能源储量的逐渐减少，人类对能源需求的不断增加，能源问题日益成为了每个国家的首要问题。
石油天然气等能源的勘探开发成为了各个国家政府工作中的重中之重。
人工地震是一种重要的勘探手段，但所接受到的地震波并非完全都是有用的信号。
在记录来自地下地层的有效波的同时也记录了来自地上和地下的各种各样的与地下地层无关的干扰波。
干扰波和有效波混在一起，不仅增加了识别有效波的困难，而且也很难提取到准确的参数，噪声严重时还容易造成错误解释。
所以，目前实际工作中的地震处理流程中，滤波成为了必不可少的一个步骤，它是进行后面处理解释工作的基础，它能否顺利并且高质量的进行将对后面的处理与解释造成非常大的影响。
这篇文章中，就针对提高信噪比的滤波技术这个环节，及其它在电脑C程序上的实现来展开讨论。
二、一维时间域滤波的设计与实现
根据软件工程的理论，像这样一个软件（程序）的生命周期包括以下几个步骤：
问题定义、可行性研究、需求分析、总体设计、详细设计、编码和单元测试、综合测试、软件维护这几个阶段。
前面已经对方法基本原理、思想进行了介绍，本章将主要按上面几个部分来进行介绍。
这次要实现的是时间域上的一维滤波。
之所以选择它来实现，其一，是因为它可以在很大程度上体现出数字滤波基本原理与思路。
其二，无论频率滤波还是时间域滤波都是对数据进行频率上的处理，本质思想相同。
其三，和其他滤波方法相比，这个更容易在VC++上实现。
（一）程序的总体设计
在结构上，程序如果要顺利运行出结果，就应该包括以下几个组成部分，即：
SEGY读取部分，傅立叶变换部分，滤波部分，输出部分。
下面逐个部分的给出程序组成部分
第一部分：
读取SEGY数据
SEGY数据是地震数据的通常标准形式。
不同的处理软件所接受的数据格式并不相同，所以一般都先存为SEGY文件，然后再转为可处理的文件格式。
这里就是把SEGY转为C语言可以接受的数据文件的步骤。
使用的时候，要注明输入文件与输出文件的名称。
SEGY文件中包含有许多信息，包括3200字节的卷头信息，400字节道头信息，最后才是数据。
前两项可以根据需要进行取舍，数据才是要读取的核心内容，所以在读取的时候要根据需要跳过一些字节。
最后处理完输出到指定的文件之内。
第二部分：
傅立叶正变换
将SEGY文件处理完之后，还需要对其进行傅立叶变换，变换为频谱，才能进行频率滤波。
不过单纯的傅立叶变换对计算机运行速度非常浪费，所以这里采用的是一种快速傅立叶算法。
第三部分：
滤波部分
如图2-1为未经滤波处理原始炮集剖面，从上图可以清晰看到有效信号中有噪音和面波的干扰。
而通过滤波所要去掉的就是图中的干扰。
图2-1原始地震数据图
由于屏幕限制，上面只截取了1到111道，每道中200到900多个样点的图形。
增益后的信号简单的说就是将地震波在地下传播中损失的能量进行补充处理后的信号。
具体过程比较复杂，这里不再过多解释，但是从图形可以更明显的看出噪音干扰和面波干扰。
利用OMEGA软件可以先划定出需要处理的区域，然后在此区域进行去面波处理。
可以得到精度更高的结果。
至于频率带宽的选取则不尽相同。
应当根据所收集数据的地区和以往经验来确定。
而且往往需要不断的调整以便得到最佳效果。
第四部分：
傅立叶逆变换
由于进行的是时间域的滤波变换，所以在滤波完成后可以不用再进行傅立叶逆变换。
但是如果是进行的频率域滤波，则需要再做傅立叶逆变换，然后才能读取出来。
第五部分：
数据转换为SEGY文件
处理完的数据一般最后都转变为SEGY的通用格式。
整个程序的结构如图2-2：
图2-2程序结构图
（二）程序编译
编译所使用的工具为创天中文VC++.操作系统及WINDOWSXP2。
程序大体结构是按照上面所说的结构来安排的。
在主程序中调用读取、傅立叶变换、确定主频、滤波等的子程序。
下面给出所使用参量的对应关系。
表2-1程序内参量对应关系
Dt
tl
L
Nj
LN
flag
T
NR
X
maxf
采样间隔
采样点数
接收道数
2的N次方
滤波因子长度
数据内包含的道数
正/负傅立叶变换
2的阶数
存放读出来的数据
所选道的主频
此外在各个模块中，zN代表噪音函数。
fft子函数用于进行傅立叶正/逆变换。
filter子函数用于进行滤波。
具体源码可以参见附属文件，具体在编译中遇见的问题也将在相关的地方有具体说明。
（三）程序测试
测试内容包括以下几个方面：
一个是在编译过程中利用编写工具查看是否存在编写错误等问题，再一个是要用实际数据来验证此程序能否按照编写的意愿进行工作，此外，还测试了本程序适用的范围。
在东方地球物理公司的实习中，我有幸得到了2个单炮数据，其中一个是另一个数据经增益处理后的结果。
这两个文件就成为测试程序的基础。
具体使用中由于增益文件的效果更明显，所以图片都是在增益文件的基础上做成的。
首先是编译问题。
经过反复的调试，目前程序可以正常进行。
但是，当出现以下情况时，可能出现错误。
1、带入傅立叶变换模块时的数组不是2的N次方；
2、在用不同数据测试时少改了参数；
3、输入输出文件时不是4个字节；
4、写入模式不正确；
5、未对多出来的数组内容赋值等等。
其次就是代入所用的数据来进行滤波处理。
在代入的过程中由于所用的文件内每道含有的样点数过多，所以一开始在程序执行中总是出错。
如下图2-3。
后来将2500点截取其中的前1000个点后问题得到解决。
通过不断的测试也发现每道样点的上限为2000。
此外在输入666道的时候不会出现计算溢出。
所以难以得到道数的上限是多少。
图2-3出错时候的表现
此外，在开始测试中由于效果不理想，里面还设计了噪音函数。
用于给输入的数据加入噪音，以便得到更明显的处理结果。
现在此函数仍然保留，但未被使用。
测试过程中，通过反复的调整频率参数来验证有效信号的范围。
确定出最佳效果。
（四）程序使用
第一步：
打开程序，如果打开成功会显示参数文件打开成功，并要求输入要用来分析主频的道号。
如图2-4
图2-4步骤一
第二步：
得到结果。
当显示如下图的时候说明运行成功。
文件被写入到了滤波波地震记录.dat。
如图2-5
图2-5步骤二
（五）小结
本章主要介绍了该程序的设计思路，实现方法，使用过程及其编译及编译中所遇到的困难。
整个过程运用到了地球物理方面的地震勘探理论，VC++编程技术，并按照软件工程的基本思想进行了程序的设计与测试，最终完成了结果。
接下来，本文就要对本程序的滤波效果进行分析。
这也是这个程序所要起到的目的。
三、滤波处理的效果分析
滤波的目的是滤去噪音，滤波的条件可以是频率，可以是振幅等等，以能否明显的区分出有用信号和无用信号作为标准。
所编辑的时间域滤波本质上仍然是一种频率滤波
（一）滤波效果
由于计算机的条件限制，无法把所有的图形一一完全列举，这里仅仅截取部分来做示例。
图3-1是原始地震数据的图形：
图3-1原始地震图象
可以看到图中有不少的杂乱信号，比如弧顶的偏下一带。
（图中仅仅截取显示的前80来道，每道的前700来点。
既然寻找到了主频，那么利用主频来改变频率范围就是接下来滤波要做的事情。
首先选择一个相对适中的范围。
通过编译的程序进行处理后就如图3-2：
图3-2第一次滤波的结果
设置的过滤范围参数为低截：
10，低通：
15，高通：
40，高截：
50（假设的主频为20，实际测量出来的主频平均值大概在20左右。
根据选择的道的不同会有变化）。
可以明显看出图形清晰了许多。
10以下，50以上的信号被过滤掉了。
被过滤掉的信号主要在中下部一带，左中部一带，信号的最外一层。
左边5，6道一带的噪音也变弱了不少。
为了作个对比，下面扩大频带范围作一下比较。
这次设置的过滤参数为低截：
5，低通：
10，高通过50，高截60。
图3-3第二次的滤波结果
图3-3明显比上一幅要多出了一些信号，成图的效果并不好，看上去很杂乱无章。
但是如果和前面一个滤波图和原始图做比较的话可以看到5以下60以上的一些信号还是被过滤掉了。
现在可以确定这个范围过大，应当缩小范围。
这回将范围缩小为低截：
10，低通：
15，高通过35，高截45，下图就是处理结果：
图3-4第三次滤波后的结果
由于有效的频带范围被降低了，所以信号也显得也更少了。
从图3-4可以看出有些信号的幅度被明显降低了。
此外主信号也开始更清晰了一些。
不过总体上来看不少信号并不明显，这可能是频率过低所导致的。
此外左边的明显噪音也并没有被过滤掉。
那个噪音的频率范围也在10到30之间。
这样就比较难以处理，因为这个范围内已经主要是有效信号了。
为了取得更好的效果，在上图的基础上再略微放大一些频率范围。
将参数设置为10，15，45，55。
图3-5第四次滤波的结果
图3-5可以看出信号略有加强，但是仍显的不太清晰。
不过经过和其他范围的对比可以看出这个范围左右还是比较适合的。
（二）软件分析
在目前世面上，用于进行地震数据处理的软件并不在少数。
例如OMEGA（美国），CGG（法国）GeoEast（中国），GRISYS（中国）等等，前两者都是由国际知名大公司所研制，后两者则是国人自主研发的地震资料处理软件。
我在东方公司研究院所实习期间所学习使用的便是OMEGA软件。
与其相比，我的程序与他们的不同体现在以下几个方面。
首先，我的程序仅仅是用来进行一维的时间域滤波。
而其他的成熟的商业软件往往还能进行动/静校正，叠加等处理。
其次，就运行平台而言，我的软件运行平台为WINDOWSXP2操作系统。
而其他的往往是在WINDOWSNT，UNIX，LINUX等工作站上运行。
因为在实际工作中只有工作站才能应付这样大规模的处理运算。
再次，在操作界面上。
现在的主流软件都已经实现了图形操作界面。
如同下面的OMEGA软件的图3-6一样。
但我这次编译的界面则相对非常简陋，是在DOS下的指令。
仅仅能够输入相关参数而已。
图3-6OMEGA制作处理任务时的操作界面
在处理效果上，它们也占有极大的优势，图3-7是一个OMEGA绘制出来的单炮数据。
图3-7OMEGA所做出来的图形
最后，从维护检修的角度来说。
他们这样的大型商业软件一旦出了问题会非常难以解决。
而我的程序对于任何熟悉C语言的使用者都可以轻松上手调试。
四、结论
本次毕业论文设计，我主要的工作就在了解滤波基本原理的基础上，运用软件工程的基本思想进行滤波软件的设计、编译和分析。
本次设计的软件可对去掉道头的地震SEGY数据进行读取和一维时间域滤波。
然后利用第三方软件将滤波的结果绘制出图形，用于分析滤波结果。
最后选择不同的滤波参数，对滤波效果进行了分析。
通过本次毕业设计，得到如下几点认识：
（1）通过本次毕业设计，进一步理解了滤波的基本原理，利用本次在C语言上编译出来的程序能够进行实际地震资料数据处理，能够将一些噪音信号过滤掉，绘成的图比原来更加清晰明了。
（2）本次毕业设计虽然只做了一个简单的一维时间域滤波，但是在整个程序设计过程中，尝试将软件工程思想融入其中，对软件工程各个环节在实际程序设计中进行了具体体现，通过本次实践，为将来进一步设计复杂的软件打下了基础。
（3）这次程序设计对我来说是一次非常有挑战的尝试，在编辑程序过程中收获到了许多原来书本上得不到的知识。
加深了对滤波处理计算的理解，对于公式有了更深刻的认识，也进一步熟悉了计算机编程技术。
了解了许多实用的编辑技巧，锻炼了自己的实际工作能力。
在编程上，已经可以熟练的使用写入、读取、打印、循环引用等语句；掌握了利用scanf命令来寻找程序错误的技巧；知道了噪音的添加方法，积累了进行应用软件测试经验；同时学会了软件开发的流程，为日后其它相关软件的开发奠定了基础。
这次设计为我将地球物理相关知识和计算机技术的结合积累了宝贵的经验。
（4）相比起预期而言，处理效果也有不理想的地方。
在理论上，一维滤波存在着自己的局限性。
对一个由许多不同频率成分的简谐波组成的地震脉冲波，经过频率滤波后，组成这个脉冲的简谐成分发生了变化，譬如某些频率成分被滤掉了，波剖面的形状要发生变化。
因此，单独的频率滤波存在不足。
目前比较常用的解决方法是使用二维视速度滤波技术、非线性滤波、最小平方滤波、小波变换以及更多维数的滤波等等。
参考文献
[1]李录明,罗省贤.信号与系统分析[M].成都：
电子科技大学出版社2002年8月
[2]王世一.数字信号处理[M].北京：
北京理工大学出版社1997年1月
[3]牟永光.地震勘探资料数字处理方法[M].北京：
石油工业出版社1986年3月
[4]程佩青.快速滤波与傅立叶变换[M].北京：
清华大学出版社1990年10月
[5]陆基孟.地震勘探原理[M].北京：
石油工业出版社1986年10月
[5]渥伊尔马滋（美）.地震资料分析[M].北京：
石油工业出版社2006年8月
[6]张文坡等.辽河油田地震资料精细处理[M].北京：
石油工业出版社2007年3月
[7]江玉乐,雷宛.地球物理数据处理教程[M].北京：
地质出版社2006年7月
[8]井上伸雄（日）.数字信号处理的应用[M].北京：
科学出版社1991年年5月
[9]张海藩.软件工程导轮[M].北京:
清华大学出版社2003年12月
[10]邓重一.滤波技术的发现状[D].衡阳:
湖南建材高等专科学校
升级会员 