毕业设计地震数据处理.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（
升级会员 