毕业设计地震数据处理.docx
- 文档编号:25434198
- 上传时间:2023-06-08
- 格式:DOCX
- 页数:52
- 大小:310.35KB
毕业设计地震数据处理.docx
《毕业设计地震数据处理.docx》由会员分享,可在线阅读,更多相关《毕业设计地震数据处理.docx(52页珍藏版)》请在冰豆网上搜索。
毕业设计地震数据处理
《地震资料数据处理》课程设计
总结报告
专业班级:
姓名:
学号:
设计时间:
指导老师:
一、设计内容………………………………………………………………
(1)褶积滤波………………………………………………
(2)快变滤波………………………………………………
(3)褶积滤波与快变滤波的比较…………………………
(4)设计高通滤波因子……………………………………
(5)频谱分析………………………………………………
(6)分析补零对振幅谱的影响……………………………
(7)线性褶积与循环褶积…………………………………
(8)最小平方反滤波………………………………………
(9)零相位转换……………………………………………
(10)最小相位转换…………………………………………
(11)静校正…………………………………………………
二、附录…………………………………………………………………………
(1)附录1:
相关程序……………………………………
(2)附录2:
相关图件……………………………………
【附录1:
有关程序】
1.褶积滤波
CCCCCCCCCCCCCCCCC褶积滤波CCCCCCCCCCCCCCCCC
PROGRAMMAIN
DIMENSIONX(100),H1(-50:
50),H2(-50:
50),Y_LOW(200),Y_BAND(200)
PARAMETER(PI=3.141592654)
CCCCCCCCH1是低通滤波因子,H2为带通滤波因子CCCCCC
REALX,H1,H2,Y_LOW,Y_BAND
REALdt,F,F1,F2
INTEGERI
dt=0.002
F=70.0
F1=10.0
F2=80.0
OPEN(1,FILE='INPUT1.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
READ(1,*)(X(I),I=1,100)
CCCCCCCCCCCCCCCCCC低通滤波器CCCCCCCCCCCCCCCCC
DO10I=-50,50
IF(I.EQ.0)THEN
H1(I)=2*F*PI/PI
ELSE
H1(I)=SIN(2*PI*F*I*dt)/(PI*I*dt)
ENDIF
10CONTINUE
CCCCCCCCCCCCCCCC输出低通滤波因子CCCCCCCCCCCCCCCC
OPEN(2,FILE='H1_LOW.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(2,*)(H1(I),I=-50,50)
CLOSE
(2)
CALLCON(X,H1,Y_LOW,100,101,200)
CCCCCCCCCCCCCCCC输出滤波后的数据CCCCCCCCCCCCCCCC
OPEN(3,FILE='Y_LOW.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(3,*)(Y_LOW(I),I=51,150)
CLOSE(3)
CCCCCCCCCCCCCCCCCC带通滤波器CCCCCCCCCCCCCCCCCCCC
DO20I=-50,50
IF(I.EQ.0)THEN
H2(I)=140
ELSE
H2(I)=SIN(2*PI*F2*I*dt)/(PI*I*dt)-SIN(2*PI*F1*I*dt)/(PI*I*dt)
ENDIF
20CONTINUE
CCCCCCCCCCCCCCC输出带通滤波因子CCCCCCCCCCCCCCCCC
OPEN(4,FILE='H2_BAND.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(4,*)(H2(I),I=-50,50)
CLOSE(4)
CALLCON(X,H2,Y_BAND,100,101,200)
CCCCCCCCCCCCCCCC输出滤波后的数据CCCCCCCCCCCCCCCCC
OPEN(5,FILE='Y_BAND.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(5,*)(Y_BAND(I),I=51,150)
CLOSE(5)
END
CCCCCCCCCCCCCCCCCCCCC褶积函数CCCCCCCCCCCCCCCCCCCC
SUBROUTINECON(A,B,C,I,J,K)
DIMENSIONA(I),B(J),C(K)
DO1K1=1,K
1C(K1)=0.0
DO2I1=1,I
DO2I2=1,J
II=I1+I2-1
2C(II)=C(II)+A(I1)*B(I2)*0.002
RETURN
END
2.快变滤波
CCCCCCCCCCCCCCC频率滤波CCCCCCCCCCCCCCCCCCCC
PROGRAMMAIN
PARAMETER(PI=3.141592654)
REALH_LOW(1:
200),H_BAND(1:
200),Y_LOW(1:
200),Y_BAND(1:
200)
REALX(1:
200)
INTEGERI,NFFT,K,NP
COMPLEXC(1:
200),TEMP(1:
200)
REALDT,DF,F,F1,F2
F=70.0
F1=10.0
F2=80.0
DT=0.002
OPEN(1,FILE='INPUT1.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
READ(1,*)(X(I),I=1,100)
NP=100
K=LOG(100.0)/LOG(2.0)
IF(2**K.LT.100)THEN
K=K+1
NFFT=2**K
ENDIF
DF=1.0/(NFFT*DT)
DO10I=101,128
X(I)=0.0
10CONTINUE
CCCCCCCCCCC数据变换成复数形式CCCCCCCCCCC
DO20I=1,NFFT
20C(I)=CMPLX(X(I),0.0)
CCCCCCCCCCC向频率域转换CCCCCCCCCCCCCCCCC
CALLFFT(NFFT,C,1.0)
CCCCCCCCCCC低通滤波因子的设计CCCCCCCCCCC
DO30I=1,NFFT/2
IF(I*DF.LE.F)THEN
H_LOW(I)=1.0
ELSE
H_LOW(I)=0.0
ENDIF
30CONTINUE
CCCCCCCCC使滤波器理想化(对称)CCCCCCCCCC
DO40I=NFFT/2+1,NFFT
F=I*DF
H_LOW(I)=H_LOW(NFFT-I+1)
40CONTINUE
CCCCCCCCCCCCCCC实现滤波CCCCCCCCCCCCCCCCC
DO50I=1,NFFT
50TEMP(I)=C(I)*H_LOW(I)
CCCCCCCCCCCCCCC向时间域变换CCCCCCCCCCCCC
CALLFFT(NFFT,TEMP,-1.0)
DO60I=1,NFFT
60Y_LOW(I)=REAL(TEMP(I))
OPEN(2,FILE='LOWPASS.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(2,*)(Y_LOW(I),I=1,NFFT)
CLOSE
(2)
CCCCCCCCCCC带通滤波因子的设计CCCCCCCCCCC
DO70I=1,NFFT/2
IF((I*DF.GE.F1).AND.(I*DF.LE.F2))THEN
H_BAND(I)=1.0
ELSE
H_BAND(I)=0.0
ENDIF
70CONTINUE
CCCCCCCCC使滤波器理想化(对称)CCCCCCCCCC
DO80I=NFFT/2+1,NFFT
F=I*DF
H_LOW(I)=H_BAND(NFFT-I+1)
80CONTINUE
CCCCCCCCCCCCCCC实现滤波CCCCCCCCCCCCCCCCC
DO90I=1,NFFT
90TEMP(I)=C(I)*H_BAND(I)
CCCCCCCCCCCCCCC向时间域变换CCCCCCCCCCCCC
CALLFFT(NFFT,TEMP,-1.0)
DO100I=1,NFFT
100Y_BAND(I)=REAL(TEMP(I))
OPEN(3,FILE='BANDPASS.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(3,*)(Y_BAND(I),I=1,NFFT)
CLOSE(3)
CLOSE
(1)
END
CCCCCCCCCCCCCCCCC子程序CCCCCCCCCCCCCCCC
CCCCCLX为输入数据的点数CCCCCCCCCCCCCCCCC
CCCCCCX为复型数组CCCCCCCCCCCCCCCCCCCCCC
CCCCCSIGNI为转换标志CCCCCCCCCCCCCCCCCCCC
SUBROUTINEFFT(LX,CX,SIGNI)
COMPLEXCX(LX),CARG,CEXP,CW,CTEMP
J=1
SC=1.0/LX
IF(SIGNI.EQ.1.0)SC=1.0
SIG=-SIGNI
DO300I=1,LX
IF(I.GT.J)GOTO100
CTEMP=CX(J)*SC
CX(J)=CX(I)*SC
CX(I)=CTEMP
100M=LX/2
200IF(J.LE.M)GOTO300
J=J-M
M=M/2
IF(M.GE.1)GOTO200
300J=J+M
L=1
400ISTEP=2*L
DO500M=1,L
CARG=(0.0,1.0)*(3.14159265*SIG*(M-1))/L
CW=CEXP(CARG)
DO500I=M,LX,ISTEP
CTEMP=CW*CX(I+L)
CX(I+L)=CX(I)-CTEMP
500CX(I)=CX(I)+CTEMP
L=ISTEP
IF(L.LT.LX)GOTO400
RETURN
END
3.褶积滤波与快变滤波的比较
CCCCCCCCCCCCC褶积滤波与递归滤波的比较CCCCCCCCCCCCCC
CCCCCCCCCCCCCCCC褶积滤波的设计CCCCCCCCCCCCCCCCCC
PROGRAMMAIN
PARAMETER(DT=0.002)
CH_NZ为非零相位滤波因子,H_Z为零相位滤波因子C
CY_CNZ为非零相位滤波结果,Y_CZ为零相位滤波结果C
REALY_CNZ(1:
100),Y_CZ(1:
200)
REALH_Z(1:
20),H_NZ(1:
20)
REALX(1:
50)
INTEGERI,J,K,N
REALA(0:
100),B(0:
100)
REALY_DNZ(1:
100),Y_DZ(1:
100)
REALX3(1:
100),X4(1:
100),X1(1:
100),X2(1:
100)
REALDF
DATAH_NZ/1.0,5.254,11.6,13.7,8.47,0.775,-3.325,-2.764,-0.364,
$1.099,0.97,0.15,-0.37,-0.344,-0.06,-0.126,0.122,0.025,-0.042,
$-0.043/
OPEN(1,FILE='INPUT3.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
READ(1,*)(X(I),I=1,50)
CLOSE
(1)
CALLCON(X,H_NZ,Y_CNZ,50,20,69)
OPEN(2,FILE='Y_CNZ.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(2,*)(Y_CNZ(i),I=1,69)
CLOSE
(2)
CALLAUTCOR(20,20,H_NZ,H_HZ)
CALLcon(X,H_CZ,Y_CZ,50,20,69)
OPEN(3,FILE='Y_CZ.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
write(3,*)(Y_CZ(i),I=1,69)
CLOSE(3)
CLOSE
(2)
CCCCCCCCCCCCCCC递归滤波的设计CCCCCCCCCCCCCCC
CX存放INPUT3里数据的数组C
CY_DNZ为非零相位递归滤波,Y_DZ为零相位递归滤波C
CX1,2,X,3,X4为运算过程中的过度变量C
A(0)=1.0
A
(1)=4.0
A
(2)=6.0
A(3)=4.0
A(4)=1.0
b(0)=0.0
B
(1)=-1.254
B
(2)=0.987
B(3)=-0.341
B(4)=0.0524
CCCCCCCCCCCCCCC对A补零CCCCCCCCCCCCCCCCCC
DO40i=5,49
40A(I)=0.0
CCCCCCCCCCCCCCC对B补零CCCCCCCCCCCCCCCCCC
DO50i=5,49
50B(I)=0.0
CCCCCCC从外部数据向X导入数据CCCCCCCCCCCC
OPEN(1,FILE='INPUT3.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
READ(1,*)(X(I),I=1,50)
CCCCCCCCCCC非零相位递归滤波的设计CCCCCCCCC
DO10I=0,49
X1(I)=0.0
X2(I)=0.0
DO20J=1,I
20X1(I)=X1(I)+A(J)*X(I-J)
IF(I.EQ.0)THEN
X2(I)=0.0
ELSE
DO30K=1,I
30X2(i)=X2(i)+B(k)*Y_DNZ(I-K)
ENDIF
Y_DNZ(I)=X1(I)-X2(I)
10CONTINUE
OPEN(2,FILE='Y_DNZ.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(2,*)(Y_DNZ(I),I=1,49)
CLOSE
(2)
CCCCCCCCCC零相位递归滤波的设计CCCCCCCCCCCC
DO60I=49,0,-1
X3(I)=0.0
X4(I)=0.0
DO70J=0,49-I
70X3(I)=X3(I)+A(J)*Y_DNZ(I+J)
IF(I.EQ.49)THEN
X4(I)=0.0
ELSE
DO80K=2,49-I
80X4(I)=X4(I)+B(K)*Y_DZ(I+K)
ENDIF
Y_DZ(I)=X3(I)-X4(I)
60CONTINUE
OPEN(3,FILE='Y_DZ.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(3,*)(Y_DZ(I),I=1,49)
CLOSE(3)
CLOSE
(1)
END
CCCCCCCCCCCCCCCCCCC褶积子程序CCCCCCCCCCCCCCCCCCC
SUBROUTINECON(A,B,C,I,J,K)
DIMENSIONA(I),B(J),C(K)
DO1K1=1,K
1C(K1)=0.0
DO2I1=1,I
DO2I2=1,J
II=I1+I2-1
2C(II)=C(II)+A(I1)*B(I2)*0.004
RETURN
END
CCCCCCCCCCCCCCCCCC自相关子程序CCCCCCCCCCCCCCCCCC
SUBROUTINEAUTCOR(J,K,C,A)
DIMENSIONC(J),A(K)
DO7N=1,K
A(N)=0.0
I=N
DO6IN=I,J
IL=IN-N+1
6A(N)=A(N)+C(IN)*C(IL)
7CONTINUE
RETURN
END
4.设计高通滤波因子
CCCCCCCCCCC高通滤波器的设计CCCCCCCCCCC
PROGRAMMAIN
PARAMETER(N=101,DT=0.004,PI=3.141592654,F=30.0)
REALH(150),H2_R(128),H2_I(128),H_S(128)
INTEGERK,NFFT
COMPLEXH2(128)
OPEN(1,FILE='H1_HIGHPASS.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
OPEN(2,FILE='H2_HIGHPASS.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
K=LOG(101.0)/LOG(2.0)
IF(2**K.LT.101)THEN
K=K+1
ENDIF
NFFT=2**K
DO10,I=-50,50
IF(I.EQ.0)THEN
H(I+51)=1.0/DT-2*F
ELSE
H(I+51)=-SIN(2*PI*30.0*I*DT)/(PI*I*DT)
ENDIF
10CONTINUE
DO20,I=102,128
20H(I)=0.0
DO30,I=1,128
30H2(I)=CMPLX(H(I),0.0)
CALLFFT(128,H2,1.0)
DO40,I=1,128
H2_R(I)=REAL(H2(I))
H2_I(I)=AIMAG(H2(I))
H_S(I)=(H2_R(I)**2+H2_I(I)**2)
40CONTINUE
WRITE(2,*)(H_S(I),I=1,128)
WRITE(1,*)(H(I),I=1,128)
CLOSE
(1)
CLOSE
(2)
END
CCCCCCCCCCCCCFFT子程序CCCCCCCCCCCCCC
SUBROUTINEFFT(LX,CX,SIGNI)
COMPLEXCX(LX),CARG,CEXP,CW,CTEMP
J=1
SC=1.0/LX
IF(SIGNI.EQ.1.0)SC=1.0
SIG=-SIGNI
DO30I=1,LX
IF(I.GT.J)GOTO10
CTEMP=CX(J)*SC
CX(J)=CX(I)*SC
CX(I)=CTEMP
10M=LX/2
20IF(J.LE.M)GOTO30
J=J-M
M=M/2
IF(M.GE.1)GOTO20
30J=J+M
L=1
40ISTEP=2*L
DO50M=1,L
CARG=(0.0,1.0)*(3.14159265*SIG*(M-1))/L
CW=CEXP(CARG)
DO50I=M,LX,ISTEP
CTEMP=CW*CX(I+L)
CX(I+L)=CX(I)-CTEMP
50CX(I)=CX(I)+CTEMP
L=ISTEP
IF(L.LT.LX)GOTO40
RETURN
END
5.频谱分析
CCCCCCCCCCCCC频谱分析CCCCCCCCCCCCCC
PROGRAMMAIN
PARAMETER(PI=3.141592654,DT=0.004)
REALXSIN(200),X(200),WAVE(200)
REALX_SIN_R(200),X_SIN_I(200)
REALX_R(200),X_I(200)
REALX_WAVE_R(200),X_WAVE_I(200)
REALSINSPRT(200),XSPRT(200),WAVESPRT(200)
REALDF
COMPLEXXSIN_C(200),X_C(200),X_WAVE_C(200)
CCCCCCCCCCCC对SIN的分析CCCCCCCCCCCC
DO100I=1,128
XSIN(I)=SIN(2*I*PI/64.0)
100CONTINUE
CCCCCCCCC是数据转换成复数形式CCCCCCCCCC
DO200I=1,128
200XSIN_C(I)=CMPLX(XSIN(I),0.0)
CALLFFT(128,XSIN_C,1.0)
DO300I=1,128
300X_SIN_R(I)=REAL(XSIN_C(I))
DO400I=1,128
400X_SIN_I(I)=AIMAG(XSIN_C(I))
OPEN(1,FILE='XSINSPRT.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
DO1I=1,128
1SINSPRT(I)=SQRT(X_SIN_R(I)**2+X_SIN_I(I)**2)
WRITE(1,*)(SINSPRT(I),I=1,128)
CLOSE
(1)
CCCCCCCCCCCCCC对脉冲信号的分析CCCCCCCCCCCCC
X
(1)=1.0
DO800I=2,128
800X(I)=0.0
DO900I=1,128
900X_C(I)=CMPLX(X(I),0.0)
CALLFFT(128,X_C,1.0)
DO1000I=1,128
1000X_R(I)=REAL(X_C(I))
DO1100I=1,128
1100X_I(I)=AIMAG(X_C(I))
OPEN(5,FILE='XSPRT.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
DO2I=1,128
2XSPRT(I)=SQRT(X_R(I)**2+X_I(I)**2)
WRITE(5,*)(XSPRT(I),I=1,128)
CLOSE(5)
CCCCCCCCCCCC对地震波WAVE的分析CCCCCCCCCCC
OPEN(3,FILE='WAVE.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
READ(3,*)(WAVE(I),I=1,128)
DO500I=1,128
500X_WAVE_C(I)=CMPLX(WAVE(I),0.0)
CALLFFT(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 毕业设计 地震 数据处理