精选MATLAB的地震数据信号的分析 doc资料.docx
《精选MATLAB的地震数据信号的分析 doc资料.docx》由会员分享，可在线阅读，更多相关《精选MATLAB的地震数据信号的分析 doc资料.docx（34页珍藏版）》请在冰豆网上搜索。
精选MATLAB的地震数据信号的分析doc资料
福建电脑2021年第2期
基于MATLAB的地震数据的分析
赵静
（中州大学工程技术学院河南郑州450044
【摘要】:
为了提高震相分析的准确性,给出了一种基于MATLAB的FIR数字滤波器的优化设计方案,并将其用于地震波数据的分析研究中。
仿真结果表明,该方法可以反映出地震波的真实信息,达到了预期目的。
【关键词】:
MATALB;数字滤波器;地震波
1、引言
地震带给人类的损失是巨大的。
地震观测资料是否准确、可靠,是地震学家进行地震预测的基础[1]。
但是地震波信号变化的不平稳性、复杂性以及各种干扰,都会给地震波的分析和预测带来严重影响,甚至导致错误结果。
为了提高地震波分析的准确性,可先画出其频谱图,然后选择合适的滤波器滤除干扰信号,最后再对数据进行分析处理。
MATLAB软件具有强大的运算处理能力,很容易实现Fourier变换和各种数字滤波器的设计,在地震数据的分析处理中起着重要作用。
本文给出了快速Fourier变换和FIR数字滤波器的MATLAB实现方法,并对一个存在干扰的地震波实例进行仿真研究。
2、快速Fourier变换的MATLAB实现
为了获取信号序列的频谱特性,可以采用离散Fourier变换（DFT。
设f（n是一个长度为M的有限长序列,则f（n的M点离散傅里叶变换定义为:
（1
由于M较大时,（1式的计算量很大,因此可以将f（n分解成许多子序列,然后利用子序列的离散Fourier变换实现整个序列的离散Fourier变换,这种方法就是快速Fourier变换（FFT。
在MATLAB中对信号序列进行快速傅立叶变换的函数为:
F=fft（f,M（2其中,f为信号序列,F为f序列的快速Fourier变换,即f的频谱特征。
3、FIR数字滤波器的MATLAB实现
数字滤波器可保留数字信号中有用频段的数据、滤去无用频段的数据。
根据实现的网络结构不同,可分为无限脉冲响应（IIR滤波器和有限脉冲响应（FIR滤波器两种。
考虑到地震波数据的特点,本文选用FIR数字滤波器,其传递函数为[2]
（3
其中,h（n是滤波器的单位脉冲响应。
若h（n是实序列,并且满足h（n=h（N-n-1或h（n=-h（N-n-1,则不但可以获得逼近平直的幅频特性,还可获得严格的线性相位特性。
利用MATLAB对FIR数字滤波器进行设计的步骤为:
（1根据地震波的频谱图确定滤波器的技术指标;
（2利用函数[M,F0,A0,W]=remezord（f,a,dev,Fs[3]估算等波纹逼近法的参数:
最低滤波器阶数M、频率向量F0、幅度向量A0和加权向量W。
其中,f是归一化频率;a为滤波器在各个频段上的幅值;dev为波纹振幅;Fs为采样频率。
（3利用函数h=remez（M,F0,A0,W完成FIR数字滤波器的设计,并调用函数filter对输入信号进行滤波。
4、仿真实例
以辽宁省营口台的数字地震记录资料为例进行地震数据的分析。
原始地震数据频谱图如图1所示,可知:
地震信号的优势频率为0.25Hz,主要频段为0~1Hz;干扰的优势频率为12.5Hz,主要干扰频段为10~15Hz。
为了滤除干扰信号,最大限度的保（下转第3页
项目:
河南省教育厅自然科学研究计划项目（2021C510002图1原始地震图2FIR带阻滤波器的数据的频谱图频率特性曲线
12
2021年第2期福建电脑
（上接第12页
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
留其中的有用信号,设计FIR带阻滤波器,参数选定为:
通带上截止频率Fp1=7,阻带下截止频率Fs1=7.1,阻带上截止频率Fs2=18.9,通带下截止频率Fp2=19,通带波纹峰值dp=0.01,阻带波纹峰值ds=0.01。
则该滤波器的频率特性曲线如图2所示。
利用该滤波器对原始地震信号进行滤波,滤波前后的时域和频域波形图分别如图3、图4所示。
可知:
干扰信号被滤除了,地震波信号很好的显示出来。
5、结束语
地震数据中的干扰信号会影响震相分析的准确性。
本文采用目前非常流行的MATLAB软件,利用快速Fourier变换和最优滤波器的设计方法,对采集到的地震数据进行分析、处理,滤除干扰并最大限度地保留有用信号[4],提高了震相分析的准确度。
该方法可用于结构地震动力分析、地震台等领域,对地震的观测、分析、预报和研究有着重要的意义。
参考文献:
[1]宋建锁.滤波在地震分析中的应用[J].防灾技术高等专科学校学报,2020,8（1:
75-79
[2]胡广书.数字信号处理[M].北京:
清华大学出版社,2003
[3]唐向红,岳恒立,郑雪峰.MATLAB及在电子信息类课程中的应用[M].北京:
电子工业出版社,2020
[4]李敬,甘延锋,黄友明.数字地震记录中干扰波的排除[J].防灾技术高等专科学校学报,2004,6（3:
20-25
图3滤波前后地震波图4滤波前后地震波时域波形图的比较
频域波形图的比较
4、结论
本文针对不一致信息系统,采用3中转换算法:
最大分布约简,分布约简和分配约简,将不一致决策信息系统转换为一致决策信息系统,然后对其进行属性约简,有效的处理了现实数据集中,由于采集能力有限等原因造成的不一致问题。
通过对约简结果的比较和分析,本文算法对输入数据没有限制,既可以处理不一致信息系统,对完备信息系统同样适用,算法执行效率及约简结果都优于经典粗糙集模型算法。
因而,本算法具有更优的实用价值。
参考文献:
[1]PawlakZ.Roughsets[J].InternationalJournalofComputerandInformationSciences,1982,1（11:
341-356.
[2]韩祯祥,张琦,文福拴.粗糙集理论及其应用[J].信息与控制.1998,27（1:
37-44.
[3]申锦标,吕跃进.粗糙集的近似约简及其算法[J].计算机技术与发展.2021,19（12:
17-20.
[4]PawlakZ.SomeIssuesonRoughSets[J].TransonRoughSetsI,Berlin:
Springer-Verlag,2004,3100:
1-58.
[5]PawlakZ,Grzymala-BusseJ,SlowinskiRetal.RoughSets[J].CommunicationoftheACM.1995,38（1:
89-95.
[6]刘少辉,吴斌等.Rough集高效算法研究[J].计算机学报.2003,26（5:
524-529.
[7]HanJC,HuXH,LinTY.ANewComputationModelforRoughSetTheoryBasedonDatabaseSystems[A].DaWaK2003,LNCS2737[C].2003,Berlin:
Springer-VerlagHeidelberg,pp.381-390..
[8]GrecoS,InuiguchiM,SlowinskiR.Fuzzyroughsetsandmul-tiple-premisegradualdecisionrules[J].IntJofApproximateRea-soning.2020,41（2:
179-211.
[9]曾黄磷.粗集理论及其应用-关于数据推理的新方法[M].重庆:
重庆大学出版社,1996.
[11]KryszkiewiczM.ComparativeStudyofAlternativeTypesofKnowledgeReductioninInconsistentSystems[J].IntJofIntelligentSystems.2001,16（1:
105-120.
[12]常犁云,王国胤等.一种粗糙集理论的属性约简及规则提取方法[J].软件学报.1999,10（11:
1206-1221.
表1本文算法与文献[12]算法属性约简结果比较
表2本文算法与文献[12]
算法执行效率比较
3
2021年第24卷第5期
ElectronicSci.＆Tech./May.15,2021
协议·算法及仿真
收稿日期:
2021-11-25
作者简介:
马茹（1986－,女,硕士研究生。
研究方向:
扩频通信干扰。
GPS信号和其干扰的Matlab仿真
马
茹1,胡建伟1,杨立成
2（1.西安电子科技大学电子工程学院,陕西西安710071;2.解放军总装备指挥技术学院光电装备系,北京101416
摘
要
文中主要以Matlab为平台,对GPS信号的信号结构（数据信号,C/A码、GPS的多种干扰方式进行仿真,并分析了不同干扰的干信比与误码率关系。
关键词GPS;干扰方式;Matlab仿真
中图分类号
文献标识码
A
文章编号1007－7820（202105－017－06
TheMatlabSimulationoftheGPSSignalandtheInterference
MaRu1,HuJianwei1,YangLicheng2
（1.SchoolofElectronicEngineering,XidianUniversity,Xi'an710071,China;
2.DepartmentofPhotoelectronicEquipment,TheAcademyofEquipmentCommand＆Technology,Beijing101416,China
Abstract
ThispapersimulatestheparameteroftheGPSsignal（datasignal,C/Acodeandthedifferentin-terferencemethodsofGPSwithMatlabasemulationmode．TherelationbetweenJNRandtheBNRindifferentinter-ferencemethodsisanalyzed．
Keywords
GPS;theinterferencemethods;Matlabsimulation
GPS接收系统最大的弱点是易受干扰性,目前,我国的GPS在军用、民用领域应用范围日益扩大,在这方面的研究工作也应加速进行。
本文详细说明了利用Matlab产生GPS信号,并对多种干扰进行仿真分析的过程。
为了对各种信号的干扰效果进行准确有效的分析,搭建仿真平台进行仿真是必要的。
将产生GPS并加不同干扰分析其性能的系统模块设计,如图1所示。
图1GPS系统干扰误码率仿真模型
由于接收机同时收到多路信号,本文简单讨论单路信号。
设接收机所处理的信号的序号为1,则对C/A码接收机,收到GPS信号可表示为[1]
R（t=2槡PD1（tCA1（tcos（ω1t+φ+J（t+N（t
（1式中,P为接收功率;D1（t为数据信号;CA1（t为信号
C/A码;J（t和N（t分别是多址干扰信号和高斯噪声。
1GPS卫星信号
假设每个导航数据包含30位[2]
每位导航数据包
含20组C/A码（1023位,每组C/A码的周期为1ms,载频取L1=154F=157542MHz,每位包含1540个载波。
当每个载波按8个采样点计算时,每个导航数据的仿真数据为30ˑ20ˑ1023ˑ1540ˑ8,由于数据量太大,导致仿真速度太慢,从而无法仿真。
数据码是数值为的数据流。
本文在编程实现数据码时,以每个导航数据为基本发送单位,用Randn函数产生长度为30bit的数据源码d,如图2所示。
图2每个导航数据的源数据
协议·算法及仿真马茹,等:
GPS信号和其干扰的Matlab仿真
数据信号信息码速率R
b
=30,数据信号的功率谱
如图3所示。
可以看出,频谱宽度约等于2R
b
=60,基
本符合发送数据源码d（t的功率谱图的理论值。
图3数据信号的功率谱图
（1C/A码的产生。
C/A码是由m序列优选对组合形成的Gold码,
由两个10级反馈移位寄存器组合产生。
图4为C/A
码产生原理图,其中的历元就是时刻,即接受卫星信
号的时刻。
图4C/A码产生原理图
在1.023MHz钟脉冲的驱动下,两个移位寄存器
分别产生码长为N=210－1=1023,周期为Nt0=
1ms的m序列G
1
和G
2
。
它们的多项式分别为
G
1
=1+x3+x10
G
2
=1+x2+x3+x6+x8+x9+x
{10（2
必须指出,这时第二个移位寄存器的输出,不是
在该移位寄存器的最后一个存储单元,而是选择该移
位寄存器中某两个存储单元的输出进行模二相加后再
输出,由此可以得到一个与G
2
（t平移等价的m序列
G
2i
;再将其与G
1
（t进行模二相加,便可产生结构不
同的C/A码,亦称Gold码。
由于G
2
（t的码元共有
1023位,故G
2
（t可能有1023种平移等价序列,这
些不同结构的伪随机码,称为一族或一组C/A码。
用pn=Golda（10调用Golda函数,产生阶数为
10的Gold序列pn,序列值－1,1分别对应了“0”码
元和“1”码元。
取前100个码片,如图5所示。
图5C/A码pn
（2C/A码的性质验证。
1相关函数特性验证。
本程序生成的Gold序列的自相关图如图6所示。
图6C/A码的自相关特性
图6验证了关于Gold码的自相关特性的结论。
C/A码中的Gold码的n=10,取偶数,除了τ=0时
有尖锐的自相关峰值,其他时刻,相关函数值为
－1－2n+22=－65,－1和－1+2n+22=63。
只在图上标
了3个点的坐标,可以很明显地看出生成的C/A码
符合自相关特性。
2平衡Gold码特性验证。
平衡的Gold序列中的一个周期内“1”码元数比“0”
码元数仅多一个,即要求这里,序列中值为“1”的码应
该要比值为“－1”的码多一个,才符合平衡码的要求。
对C/A码求和,得出sum（pn=1,可以验证生成的
C/A码符合平衡特性。
具体操作如图7所示。
马茹,等:
GPS信号和其干扰的Matlab仿真
协议·算法及仿真
1.3GPS信号的频谱分析
导航数据用daohang_d=kron（d,ones（1,2*1023产生每字对应2组C/A码,用pn=kron（ones（1,2*30,pn来实现导航数据对应的扩频序列pn,即30个每字包含2个C/A码组的pn序列c（t;用kuopin_d=daohang_d.*pn对30位导航数据进行扩展频谱调制,即d（tc（t;Sin=sin（[
0:
2*pi/8:
2*pi*7/8]是采样点为8的正弦波,carr=kron（ones（1,5,Sin表示导航数据中每个码元为5个载波,用Send_data=kron（kuopin_d,carr来表示对扩频后的数据进行载波调制,即发送数据s（t。
每个导航数据的仿真数据为30ˑ2ˑ1023ˑ5ˑ8=2.4552MB。
发送数据局部图如图8所示。
图8
发送端数据局部
本文以一个导航数据为单位,伪噪声码速率Rc=30ˑ2ˑ1023=61.380kHz。
在PSK调制的情况下,发送的导航数据s（t的带宽等于伪随机码速率的两倍,即BRF=2Rc,而几乎与数字信号d（t的码速率无关
[3]
。
图9为扩频码功率谱图。
图9
扩频码信号的功率谱
2GPS干扰效果定义
当GPS接收机受到干扰时,才能计算干扰机的
干扰效果。
描述这一干扰效果的最常用方法是用进入接收机内部的干扰信号功率除以GPS信号功率,所得到的效果比来表示。
这个比率称为干信比,用J/S来表示,单位为dB
[4]
。
GPS接收机接收到的信号功率
GPS信号从发射机向GPS接收机直视远程发射。
它到达GPS接收机的信号功率为
S=PT+GT－32－20lgF－20lgDS+Gr
（3
GPS接收机接收到的干扰功率
根据信号传输特性,干扰信号的传输是单程的。
一般情况下,干扰目标无论是通信接收机还是GPS接收机,其干扰信号的特征都是相同的。
在测定干信比是,重要的是计算接入GPS接收机工作频道内的干扰信号功率。
在上述条件下,到达GPS接收机输入端的干扰信号功率为
J=PJ+GJ－32－20lgF－20lgDJ+GrJ
（4
干信比的计算
在GPS接收机带宽内,干信比在纵坐标上是线
性的。
这是在GPS接收机频带为理想频宽并调谐于
GPS信号频率的前提下得到的。
由于J和S都以“dB”来表示,它们的功率比与其“dB”比是大体相同的。
对于单程的GPS信号传输情况来说,干信比公式则可以用分贝计算得到,如式（5所示。
J－S=PJ+GJ－32－20lgF－20lgDJ+GrJ－[PT+GT－32－20lgF－20lgGS+Gr]=PJ－PT+GJ－GT－20lgDJ+20lgDS+GrJ－Gr
（5
3
各种干扰的仿真
带限平稳高斯随机过程干扰分析
高斯白噪声是最常用也最容易获得的压制式干扰信号,其随机性有利于造成信号幅度的随机变化,使
得信号码元极性被错误判断,随着干扰信号的功率增加,误码率也逐渐增大。
当误码率达到一定程度时,可切断通信链路。
假设干扰信号为
J（t=n（tcos（2πfjt+φ
（6
n（t为高斯白噪声,进入直扩系统和本地扩频
码、本地载波相乘后滤除高频分量,这时的噪声表达式为
J（t=n（tc（tcosφ'
（7
协议·算法及仿真
马茹,等:
GPS信号和其干扰的Matlab仿真
其概率分布函数为
P[J'（t＜n]=P[c（tn（tcosφ'＜n]=
12∫
n
－ɕ
cosφ'12槡
πσen2（t2σdn（t+
12∫
ɕ
－n
cosφ'
1
2槡
πσen2（t2σdn（t（8
其中n为判决门限,任意实数,因此,其概率密度函数为
f[J'（t]=
dP（J'（t＜ndn=12槡
πσeJ2（t2σ2（9
所以高斯信号通过直扩通信系统后,J'（t仍为高斯
分布。
本系统的解调采取DBPSK方法,可计算得到系统的总误码率为
Pe=P（1Pe1+P（0Pe0=1
2
e－r
（10
其中,P（1是发送信息码元为“1”的概率;P（0是
发送信息码元为“0”
的概率;r为干信比。
仿真可得到在不同干信比下,干扰信号为高斯白噪声时,对应
的误码率曲线如图10所示。
图10
高斯白噪声误码率分析
由图10可见,信噪比越大,误码率越小。
以信
单频信号时比较常见的一种信号,单音干扰因为频谱窄,功率集中,进入宽带接收机的通频带后很容易引起阻塞,从而破坏接收机正常的工作。
单音干扰将一个大功率连续波音频放置在频谱的某一位置。
该位置通常取决于目标,目的是优化该位置,使其在某种意义上有利于干扰机。
当干扰机置于Wss内时,认为接收机的任意带宽滤波器对于干扰信号都不起作用。
换句话说,它离滤波带边缘足够远,滤波就不能明显地使它衰减。
频谱图如图11所示。
假设单音干扰的干扰信号为
J（t=2P槡jcos（2πfjt+φ
（11
图11
单音干扰频谱示意图
其中,Pj为干扰信号功率;fj为干扰信号频率;φ为[0,2π]均匀分布的随机变量。
单频信号在接收端和本地扩频码、本地载波相乘后,滤除高频分量。
这时的噪声表达式为
J（t=2P槡j
2N
cosφ
（12
其中,N为扩频码的周期。
则单音干扰下GPS系统的
误码率为
Pe=P（1P
2P槡j2N
cosφ+2P槡s2
（
＜0+P（0P2P槡j2Ncosφ－2P
槡s2
≥（
0=12πarccosN
PsP槡
j
（13
其中,P（1是发送信息码元为“1”的概率;P（0是
发送信息码元为“0”的概率;Ps是GPS信号功率。
仿真得出单音干扰的干信比与无码率曲线图如图12
所示。
图12单音干扰误码率分析
多音干扰分析
13所示。
马茹,等:
GPS信号和其干扰的Matlab仿真协议·算法及仿真
图13多音干扰频谱示意图
在图13（a中,5个等幅等间隔分布的音频,以
f
为中心,两边对称;图13（b显示了更普遍的情
况,即3个非等幅不等间隔分布的音频,不以f
对称。
对于多音干扰,此时有多个不同的音频和不同幅
度的单音干扰。
通过上面对单音干扰的分析,可将多
音干扰表示为
J（t=Acos（2π（f
J
ʃk
1
f
J
ʃk
2
f
J
…,ʃk
n
f
J
t+φ
（14
多音干扰的误码率曲线如图14所示。
图14多音干扰频谱示意图
从图14中可以看出,信噪比越大,误码率越小。
信噪比为－20dB时误码率曲线不明显,以信噪比为
－10dB为例分析,干信比越大,即干扰信号功率越
大,系统误码率越高。
干信比为15dB以后,误码率
基本保持在最大不变。
说明干信比超过一定范围时,
其变化已经不影响误码率。
已调波表示为
J（t=U
cos[w
t+θ（t]（15
其中,θ（t=K
FM∫
t
u
n
（tdt=K
FM
η（t;KFM为调频参
数,调制噪声电压u
n
（t为正态分布。
因为η（t=∫t0un（tdt,相当于函数un（t通过积
分网络,其功率谱密度为
G
η
（Ω=
G
n
（Ω
Ω2
（16
式中,G
n
（Ω为函数u
n
（t的功率谱密度;Ω2为线性
系统的幅频特性。
线性调频信号的能量90%以上都集中在带宽B
内,而且是均匀分布的。
线性调频信号容易获得、谱
宽和中心频率都易控制且能量集中,适合用于对直扩
系统的干扰。
其干信比和误码率曲线如图15所示。
图15调频噪声干扰频谱示意图
从仿真图15看出,信噪比越大,误码率越小。
以信噪比为－10dB为例分析,当干信比＜10dB时,
误码率全为0。
说明调频干扰的干扰功率,达到同样
干扰效果时要比其他干扰高。
当干信比为1030时,
误码率成曲线上升趋势。
说明误码率随着干扰信号增
加而增大。
干信比大于一定量（28dB时,误码率已
经基本保持不变。
噪声调幅是以视频噪声对射频振荡信号进行调
制,使射频振荡的振幅随视频噪声作不规则变化。
调幅噪声干扰表达式为
J（t=[U
+u
n
（t]cosω
t=A（tcosω
t（20
式中,U
为射频载波振幅;u
n
（t为调幅噪声。
假设调制噪声u
n
升级会员 