《地震数据处理课程设计-浙江大学.docx》由会员分享，可在线阅读，更多相关《地震数据处理课程设计-浙江大学.docx（21页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、地震数据处理课程设计-浙江大学 地震资料数据处理课程设计 总结报告 专业班级： 姓名： 学号： 设计时间： 指导老师： 目录 一、设计内容 （1）褶积滤波 （2）快变滤波 （3）褶积滤波与快变滤波的比较 （4）设计高通滤波因子 （5）频谱分析 （6）分析补零对振幅谱的影响 （7）线性褶积与循环褶积 （8）最小平方反滤波 （9）零相位转换 （10）最小相位转换 （11）静校正 二、附录 （1）附录1：相关程序 （2）附录2：相关图件 1.褶积滤波 CCCCCCCCCCCCCCCCC 褶积滤波CCCCCCCCCCCCCCCCC PROGRAM MAIN DIMENSION X(100),H1(-5
2、0:50),H2(-50:50),Y_LOW(200),Y_BAND(200) PARAMETER (PI=3.141592654) CCCCCCCC H1是低通滤波因子，H2为带通滤波因子CCCCCC REAL X,H1,H2,Y_LOW,Y_BAND REAL dt,F,F1,F2 INTEGER I dt=0.002 F=70.0 F1=10.0 F2=80.0 OPEN(1,FILE=INPUT1.DA T,FORM=FORMATTED,STATUS=UNKNOWN) READ(1,*)(X(I),I=1,100) CCCCCCCCCCCCCCCCCC低通滤波器CCCCCCCCCCCC
3、CCCCC DO 10 I=-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) END IF 10 CONTINUE CCCCCCCCCCCCCCCC输出低通滤波因子CCCCCCCCCCCCCCCC OPEN(2,FILE=H1_LOW.DAT,FORM=FORMATTED,STATUS=UNKNOWN) WRITE(2,*)(H1(I),I=-50,50) CLOSE(2) CALL CON(X,H1,Y_LOW,100,101,200) CCCCCCCCCCCCCCCC输出滤波后的数据CC
4、CCCCCCCCCCCCCC OPEN(3,FILE=Y_LOW.DA T,FORM=FORMATTED,STATUS=UNKNOWN) WRITE(3,*)(Y_LOW(I),I=51,150) CLOSE(3) CCCCCCCCCCCCCCCCCC带通滤波器CCCCCCCCCCCCCCCCCCCC DO 20 I=-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) END IF 20 CONTINUE CCCCCCCCCCCCCCC输出带通滤波因
5、子CCCCCCCCCCCCCCCCC OPEN(4,FILE=H2_BAND.DAT,FORM=FORMA TTED,STATUS=UNKNOWN) WRITE(4,*)(H2(I),I=-50,50) CLOSE(4) CALL CON(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 CCCCCCCCCCCCCC
6、CCCCCCC褶积函数CCCCCCCCCCCCCCCCCCCC SUBROUTINE CON(A,B,C,I,J,K) DIMENSION A(I),B(J),C(K) DO 1 K1=1,K 1 C(K1)=0.0 DO 2 I1=1,I DO 2 I2=1,J II=I1+I2-1 2 C(II)=C(II)+A(I1)*B(I2)*0.002 RETURN END 2.快变滤波 CCCCCCCCCCCCCCC频率滤波CCCCCCCCCCCCCCCCCCCC PROGRAM MAIN PARAMETER (PI=3.141592654) REAL H_LOW(1:200),H_BAND(1
7、:200),Y_LOW(1:200),Y_BAND(1:200) REAL X(1:200) INTEGER I,NFFT,K,NP COMPLEX C(1:200),TEMP(1:200) REAL DT,DF,F,F1,F2 F=70.0 F1=10.0 F2=80.0 DT=0.002 OPEN(1,FILE=INPUT1.DA T,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 END IF DF=
8、1.0/(NFFT*DT) DO 10 I=101,128 X(I)=0.0 10 CONTINUE DO 20 I=1,NFFT 20 C(I)=CMPLX(X(I),0.0) CCCCCCCCCCC 向频率域转换CCCCCCCCCCCCCCCCC CALL FFT(NFFT,C,1.0) CCCCCCCCCCC 低通滤波因子的设计CCCCCCCCCCC DO 30 I=1,NFFT/2 IF(I*DF.LE.F)THEN H_LOW(I)=1.0 ELSE H_LOW(I)=0.0 END IF 30 CONTINUE CCCCCCCCC 使滤波器理想化（对称）CCCCCCCCCC DO
9、40 I=NFFT/2+1,NFFT F=I*DF H_LOW(I)=H_LOW(NFFT-I+1) 40 CONTINUE CCCCCCCCCCCCCCC 实现滤波CCCCCCCCCCCCCCCCC DO 50 I=1 , NFFT 50 TEMP(I)=C(I)*H_LOW(I) CCCCCCCCCCCCCCC 向时间域变换CCCCCCCCCCCCC CALL FFT(NFFT,TEMP,-1.0) DO 60 I=1 , NFFT 60 Y_LOW(I)=REAL(TEMP(I) OPEN(2,FILE=LOWPASS.DAT,FORM=FORMA TTED,STATUS=UNKNOWN
10、) WRITE(2,*)(Y_LOW(I),I=1,NFFT) CLOSE(2) CCCCCCCCCCC 带通滤波因子的设计CCCCCCCCCCC DO 70 I=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 END IF 70 CONTINUE CCCCCCCCC 使滤波器理想化（对称）CCCCCCCCCC DO 80 I=NFFT/2+1,NFFT F=I*DF H_LOW(I)=H_BAND(NFFT-I+1) 80 CONTINUE CCCCCCCCCCCCCCC 实现滤波CC
11、CCCCCCCCCCCCCCC DO 90 I=1 , NFFT 90 TEMP(I)=C(I)*H_BAND(I) CALL FFT(NFFT,TEMP,-1.0) DO 100 I=1 , NFFT 100 Y_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 CCCCC LX 为输入数据的点数CCCCCCCC
12、CCCCCCCCC CCCCC CX 为复型数组CCCCCCCCCCCCCCCCCCCCCC CCCCC SIGNI 为转换标志CCCCCCCCCCCCCCCCCCCC SUBROUTINE FFT(LX,CX,SIGNI) COMPLEX CX(LX),CARG,CEXP,CW,CTEMP J=1 SC=1.0/LX IF(SIGNI.EQ.1.0)SC=1.0 SIG=-SIGNI DO 300 I=1,LX IF(I.GT.J)GO TO 100 CTEMP=CX(J)*SC CX(J)=CX(I)*SC CX(I)=CTEMP 100 M=LX/2 200 IF(J.LE.M)GO T
13、O 300 J=J-M M=M/2 IF(M.GE.1)GO TO 200 300 J=J+M L=1 400 ISTEP=2*L DO 500 M=1,L CARG=(0.0,1.0)*(3.14159265*SIG*(M-1)/L CW=CEXP(CARG) DO 500 I=M,LX,ISTEP CTEMP=CW*CX(I+L) CX(I+L)=CX(I)-CTEMP 500 CX(I)=CX(I)+CTEMP L=ISTEP IF(L.LT.LX)GO TO 400 RETURN END 3.褶积滤波与快变滤波的比较 CCCCCCCCCCCCC 褶积滤波与递归滤波的比较CCCCCCCC
14、CCCCCC CCCCCCCCCCCCCCCC 褶积滤波的设计CCCCCCCCCCCCCCCCCC PROGRAM MAIN PARAMETER(DT=0.002) C H_NZ为非零相位滤波因子,H_Z为零相位滤波因子 C C Y_CNZ为非零相位滤波结果,Y_CZ为零相位滤波结果 C REAL Y_CNZ(1:100),Y_CZ(1:200) REAL H_Z(1:20),H_NZ(1:20) REAL X(1:50) INTEGER I,J,K,N REAL A(0:100),B(0:100) REAL Y_DNZ(1:100),Y_DZ(1:100) REAL X3(1:100),X4
15、(1:100),X1(1:100),X2(1:100) REAL DF DATA H_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) CALL CON(X,H_NZ,Y_CNZ,50,20,69) OPEN(2,FI
16、LE=Y_CNZ.DAT,FORM=FORMATTED,STATUS=UNKNOWN) WRITE(2,*)(Y_CNZ(i),I=1,69) CLOSE(2) CALL AUTCOR(20,20,H_NZ,H_HZ) CALL con(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 C X 存放INPUT3里数据的数组 C C
17、Y_DNZ 为非零相位递归滤波，Y_DZ 为零相位递归滤波C C X1,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 DO 40 i=5,49 40 A(I)=0.0 CCCCCCCCCCCCCCC 对B 补零CCCCCCCCCCCCCCCCCC DO 50 i=5,49 50 B(I)=0.0 CCCCCCC 从外
18、部数据向X 导入数据CCCCCCCCCCCC OPEN(1,FILE=INPUT3.DA T,FORM=FORMATTED,STATUS=UNKNOWN) READ(1,*)(X(I),I=1,50) CCCCCCCCCCC 非零相位递归滤波的设计CCCCCCCCC DO 10 I=0,49 X1(I)=0.0 X2(I)=0.0 DO 20 J=1,I 20 X1(I)=X1(I)+A(J)*X(I-J) IF(I.EQ.0) THEN X2(I)=0.0 ELSE DO 30 K=1,I 30 X2(i)=X2(i)+B(k)*Y_DNZ(I-K) END IF Y_DNZ(I)=X1(I
19、)-X2(I) 10 CONTINUE OPEN(2,FILE=Y_DNZ.DAT,FORM=FORMATTED,STATUS=UNKNOWN) WRITE(2,*)(Y_DNZ(I),I=1,49) CLOSE(2) CCCCCCCCCC 零相位递归滤波的设计CCCCCCCCCCCC DO 60 I=49,0,-1 X3(I)=0.0 X4(I)=0.0 DO 70 J=0,49-I 70 X3(I)=X3(I)+A(J)*Y_DNZ(I+J) IF(I.EQ.49) THEN X4(I)=0.0 ELSE DO 80 K=2,49-I 80 X4(I)=X4(I)+B(K)*Y_DZ(I+
20、K) END IF Y_DZ(I)=X3(I)-X4(I) 60 CONTINUE 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 SUBROUTINE CON(A,B,C,I,J,K) DIMENSION A(I),B(J),C(K) DO 1 K1=1,K 1 C(K1)=0.0 DO 2 I1=1,I DO 2 I2=1,J II=I1+I2-1
21、2 C(II)=C(II)+A(I1)*B(I2)*0.004 RETURN END CCCCCCCCCCCCCCCCCC 自相关子程序CCCCCCCCCCCCCCCCCC SUBROUTINE AUTCOR(J,K,C,A) DIMENSION C(J),A(K) DO 7 N=1,K A(N)=0.0 I=N DO 6 IN=I,J IL=IN-N+1 6 A(N)=A(N)+C(IN)*C(IL) 7 CONTINUE RETURN END 4.设计高通滤波因子 CCCCCCCCCCC 高通滤波器的设计CCCCCCCCCCC PROGRAM MAIN PARAMETER (N=101,D
22、T=0.004,PI=3.141592654,F=30.0) REAL H(150),H2_R(128),H2_I(128),H_S(128) INTEGER K,NFFT COMPLEX H2(128) OPEN(1,FILE=H1_HIGHPASS.DAT,FORM=FORMATTED,STATUS=UNKNOWN) OPEN(2,FILE=H2_HIGHPASS.DA T,FORM=FORMATTED,STATUS=UNKNOWN) K=LOG(101.0)/LOG(2.0) IF(2*K.LT.101)THEN K=K+1 ENDIF NFFT=2*K DO 10,I=-50,50 I
23、F(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) END IF 10 CONTINUE DO 20,I=102,128 20 H(I)=0.0 DO 30,I=1,128 30 H2(I)=CMPLX(H(I),0.0) CALL FFT(128,H2,1.0) DO 40,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) 40 CONTINUE WRITE(2,*)(H_S(I),I=1,
24、128) WRITE(1,*)(H(I),I=1,128) CLOSE(1) CLOSE(2) END CCCCCCCCCCCCC FFT 子程序CCCCCCCCCCCCCC SUBROUTINE FFT(LX,CX,SIGNI) COMPLEX CX(LX),CARG,CEXP,CW,CTEMP J=1 SC=1.0/LX IF(SIGNI.EQ.1.0)SC=1.0 SIG=-SIGNI DO 30 I=1,LX IF(I.GT.J)GO TO 10 CTEMP=CX(J)*SC CX(J)=CX(I)*SC CX(I)=CTEMP 10 M=LX/2 20 IF(J.LE.M)GO TO
25、 30 J=J-M M=M/2 IF(M.GE.1)GO TO 20 30 J=J+M L=1 40 ISTEP=2*L DO 50 M=1,L CARG=(0.0,1.0)*(3.14159265*SIG*(M-1)/L CW=CEXP(CARG) DO 50 I=M,LX,ISTEP CTEMP=CW*CX(I+L) CX(I+L)=CX(I)-CTEMP 50 CX(I)=CX(I)+CTEMP L=ISTEP IF(L.LT.LX)GO TO 40 RETURN END 5.频谱分析 CCCCCCCCCCCCC 频谱分析CCCCCCCCCCCCCC PROGRAM MAIN PARAM
26、ETER (PI=3.141592654,DT=0.004) REAL XSIN(200),X(200),WA VE(200) REAL X_SIN_R(200),X_SIN_I(200) REAL X_R(200),X_I(200) REAL X_WA VE_R(200),X_WA VE_I(200) REAL SINSPRT(200),XSPRT(200),WA VESPRT(200) REAL DF COMPLEX XSIN_C(200),X_C(200),X_WA VE_C(200) CCCCCCCCCCCC 对SIN 的分析CCCCCCCCCCCC DO 100 I=1,128 XS
27、IN(I)=SIN(2*I*PI/64.0) 100 CONTINUE CCCCCCCCC 是数据转换成复数形式CCCCCCCCCC DO 200 I=1,128 200 XSIN_C(I)=CMPLX(XSIN(I),0.0) CALL FFT(128,XSIN_C,1.0) DO 300 I=1,128 300 X_SIN_R(I)=REAL(XSIN_C(I) DO 400 I=1,128 400 X_SIN_I(I)=AIMAG(XSIN_C(I) OPEN(1,FILE=XSINSPRT.DAT,FORM=FORMATTED,STATUS=UNKNOWN) DO 1 I=1,128 1 SINSPRT(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 DO 800 I=2,128 800 X(I)=0.0 DO 900 I=1,128 900 X_C(I)=CMPLX(X(I),0.0) CALL FFT(128,X_C,1.0) DO 1000 I=1,128 1000 X_R(I)=REAL(X_C(I)