地震数据处理课程设计(报告)

《地震资料数据处理》课程设计总结报告专业班级：姓名：学号：设计时间：指导老师：2011年5月30日目录一、设计内容………………………………………………………………（1）褶积滤波………………………………………………（2）快变滤波………………………………………………（3）褶积滤波与快变滤波的比较…………………………（4）设计高通滤波因子……………………………………（5）频谱分析………………………………………………（6）分析补零对振幅谱的影响……………………………（7）线性褶积与循环褶积…………………………………（8）最小平方反滤波………………………………………（9）零相位转换……………………………………………（10）最小相位转换…………………………………………（11）静校正…………………………………………………二、附录…………………………………………………………………………（1）附录1：相关程序……………………………………（2）附录2：相关图件……………………………………【附录1：有关程序】1.褶积滤波CCCCCCCCCCCCCCCCC 褶积滤波CCCCCCCCCCCCCCCCCPROGRAM MAINDIMENSION X(100),H1(-50:50),H2(-50:50),Y_LOW(200),Y_BAND(200)PARAMETER (PI=3.141592654)CCCCCCCC H1是低通滤波因子，H2为带通滤波因子CCCCCCREAL X,H1,H2,Y_LOW,Y_BANDREAL dt,F,F1,F2INTEGER Idt=0.002F=70.0F1=10.0F2=80.0OPEN(1,FILE='INPUT1.DA T',FORM='FORMATTED',STATUS='UNKNOWN') READ(1,*)(X(I),I=1,100)CCCCCCCCCCCCCCCCCC低通滤波器CCCCCCCCCCCCCCCCC DO 10 I=-50,50IF (I.EQ.0)THENH1(I)=2*F*PI/PIELSEH1(I)=SIN(2*PI*F*I*dt)/(PI*I*dt)END IF10 CONTINUECCCCCCCCCCCCCCCC输出低通滤波因子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输出滤波后的数据CCCCCCCCCCCCCCCC 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,50IF(I.EQ.0)THENH2(I)=140ELSEH2(I)=SIN(2*PI*F2*I*dt)/(PI*I*dt)-SIN(2*PI*F1*I*dt)/(PI*I*dt)END IF20 CONTINUECCCCCCCCCCCCCCC输出带通滤波因子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)ENDCCCCCCCCCCCCCCCCCCCCC褶积函数CCCCCCCCCCCCCCCCCCCC SUBROUTINE CON(A,B,C,I,J,K)DIMENSION A(I),B(J),C(K)DO 1 K1=1,K1 C(K1)=0.0DO 2 I1=1,IDO 2 I2=1,JII=I1+I2-12 C(II)=C(II)+A(I1)*B(I2)*0.002RETURNEND2.快变滤波CCCCCCCCCCCCCCC频率滤波CCCCCCCCCCCCCCCCCCCC PROGRAM MAINPARAMETER (PI=3.141592654)REAL H_LOW(1:200),H_BAND(1:200),Y_LOW(1:200),Y_BAND(1:200)REAL X(1:200)INTEGER I,NFFT,K,NPCOMPLEX C(1:200),TEMP(1:200)REAL DT,DF,F,F1,F2F=70.0F1=10.0F2=80.0DT=0.002OPEN(1,FILE='INPUT1.DA T',FORM='FORMATTED',STATUS='UNKNOWN') READ(1,*)(X(I), I=1,100)NP=100K=LOG(100.0)/LOG(2.0)IF(2**K.LT.100)THENK=K+1NFFT=2**KEND IFDF=1.0/(NFFT*DT)DO 10 I=101,128X(I)=0.010 CONTINUEDO 20 I=1,NFFT20 C(I)=CMPLX(X(I),0.0)CCCCCCCCCCC 向频率域转换CCCCCCCCCCCCCCCCCCALL FFT(NFFT,C,1.0)CCCCCCCCCCC 低通滤波因子的设计CCCCCCCCCCCDO 30 I=1,NFFT/2IF(I*DF.LE.F)THENH_LOW(I)=1.0ELSEH_LOW(I)=0.0END IF30 CONTINUECCCCCCCCC 使滤波器理想化（对称）CCCCCCCCCCDO 40 I=NFFT/2+1,NFFTF=I*DFH_LOW(I)=H_LOW(NFFT-I+1)40 CONTINUECCCCCCCCCCCCCCC 实现滤波CCCCCCCCCCCCCCCCCDO 50 I=1 , NFFT50 TEMP(I)=C(I)*H_LOW(I)CCCCCCCCCCCCCCC 向时间域变换CCCCCCCCCCCCCCALL FFT(NFFT,TEMP,-1.0)DO 60 I=1 , NFFT60 Y_LOW(I)=REAL(TEMP(I))OPEN(2,FILE='LOWPASS.DAT',FORM='FORMA TTED',STATUS='UNKNOWN') WRITE(2,*)(Y_LOW(I),I=1,NFFT)CLOSE(2)CCCCCCCCCCC 带通滤波因子的设计CCCCCCCCCCCDO 70 I=1,NFFT/2IF((I*DF.GE.F1).AND.(I*DF.LE.F2))THENH_BAND(I)=1.0ELSEH_BAND(I)=0.0END IF70 CONTINUECCCCCCCCC 使滤波器理想化（对称）CCCCCCCCCCDO 80 I=NFFT/2+1,NFFTF=I*DFH_LOW(I)=H_BAND(NFFT-I+1)80 CONTINUECCCCCCCCCCCCCCC 实现滤波CCCCCCCCCCCCCCCCCDO 90 I=1 , NFFT90 TEMP(I)=C(I)*H_BAND(I)CALL FFT(NFFT,TEMP,-1.0)DO 100 I=1 , NFFT100 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)ENDCCCCCCCCCCCCCCCCC子程序CCCCCCCCCCCCCCCCCCCCC LX 为输入数据的点数CCCCCCCCCCCCCCCCCCCCCC CX 为复型数组CCCCCCCCCCCCCCCCCCCCCCCCCCC SIGNI 为转换标志CCCCCCCCCCCCCCCCCCCCSUBROUTINE FFT(LX,CX,SIGNI)COMPLEX CX(LX),CARG,CEXP,CW,CTEMPJ=1SC=1.0/LXIF(SIGNI.EQ.1.0)SC=1.0SIG=-SIGNIDO 300 I=1,LXIF(I.GT.J)GO TO 100CTEMP=CX(J)*SCCX(J)=CX(I)*SCCX(I)=CTEMP100 M=LX/2200 IF(J.LE.M)GO TO 300J=J-MM=M/2IF(M.GE.1)GO TO 200300 J=J+ML=1400 ISTEP=2*LDO 500 M=1,LCARG=(0.0,1.0)*(3.14159265*SIG*(M-1))/LCW=CEXP(CARG)DO 500 I=M,LX,ISTEPCTEMP=CW*CX(I+L)CX(I+L)=CX(I)-CTEMP500 CX(I)=CX(I)+CTEMPL=ISTEPIF(L.LT.LX)GO TO 400RETURNEND3.褶积滤波与快变滤波的比较CCCCCCCCCCCCC 褶积滤波与递归滤波的比较CCCCCCCCCCCCCC CCCCCCCCCCCCCCCC 褶积滤波的设计CCCCCCCCCCCCCCCCCC PROGRAM MAINPARAMETER(DT=0.002)C H_NZ为非零相位滤波因子,H_Z为零相位滤波因子 CC Y_CNZ为非零相位滤波结果,Y_CZ为零相位滤波结果 CREAL 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,NREAL A(0:100),B(0:100)REAL Y_DNZ(1:100),Y_DZ(1:100)REAL X3(1:100),X4(1:100),X1(1:100),X2(1:100)REAL DFDATA 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,FILE='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 递归滤波的设计CCCCCCCCCCCCCCCC X 存放INPUT3里数据的数组 CC Y_DNZ 为非零相位递归滤波，Y_DZ 为零相位递归滤波CC X1,2,X,3,X4 为运算过程中的过度变量 CA(0)=1.0A(1)=4.0A(2)=6.0A(3)=4.0A(4)=1.0b(0)=0.0B(1)=-1.254B(2)=0.987B(3)=-0.341B(4)=0.0524CCCCCCCCCCCCCCC 对A 补零CCCCCCCCCCCCCCCCCC DO 40 i=5,4940 A(I)=0.0CCCCCCCCCCCCCCC 对B 补零CCCCCCCCCCCCCCCCCC DO 50 i=5,4950 B(I)=0.0CCCCCCC 从外部数据向X 导入数据CCCCCCCCCCCCOPEN(1,FILE='INPUT3.DA T',FORM='FORMATTED',STATUS='UNKNOWN') READ(1,*)(X(I),I=1,50)CCCCCCCCCCC 非零相位递归滤波的设计CCCCCCCCCDO 10 I=0,49X1(I)=0.0X2(I)=0.0DO 20 J=1,I20 X1(I)=X1(I)+A(J)*X(I-J)IF(I.EQ.0) THENX2(I)=0.0ELSEDO 30 K=1,I30 X2(i)=X2(i)+B(k)*Y_DNZ(I-K)END IFY_DNZ(I)=X1(I)-X2(I)10 CONTINUEOPEN(2,FILE='Y_DNZ.DAT',FORM='FORMATTED',STATUS='UNKNOWN') WRITE(2,*)(Y_DNZ(I),I=1,49)CLOSE(2)CCCCCCCCCC 零相位递归滤波的设计CCCCCCCCCCCCDO 60 I=49,0,-1X3(I)=0.0X4(I)=0.0DO 70 J=0,49-I70 X3(I)=X3(I)+A(J)*Y_DNZ(I+J)IF(I.EQ.49) THENX4(I)=0.0ELSEDO 80 K=2,49-I80 X4(I)=X4(I)+B(K)*Y_DZ(I+K)END IFY_DZ(I)=X3(I)-X4(I)60 CONTINUEOPEN(3,FILE='Y_DZ.DAT',FORM='FORMATTED',STATUS='UNKNOWN') WRITE(3,*)(Y_DZ(I),I=1,49)CLOSE(3)CLOSE(1)ENDCCCCCCCCCCCCCCCCCCC 褶积子程序CCCCCCCCCCCCCCCCCCC SUBROUTINE CON(A,B,C,I,J,K)DIMENSION A(I),B(J),C(K)DO 1 K1=1,K1 C(K1)=0.0DO 2 I1=1,IDO 2 I2=1,JII=I1+I2-12 C(II)=C(II)+A(I1)*B(I2)*0.004RETURNENDCCCCCCCCCCCCCCCCCC 自相关子程序CCCCCCCCCCCCCCCCCC SUBROUTINE AUTCOR(J,K,C,A)DIMENSION C(J),A(K)DO 7 N=1,KA(N)=0.0I=NDO 6 IN=I,JIL=IN-N+16 A(N)=A(N)+C(IN)*C(IL)7 CONTINUERETURNEND4.设计高通滤波因子CCCCCCCCCCC 高通滤波器的设计CCCCCCCCCCCPROGRAM MAINPARAMETER (N=101,DT=0.004,PI=3.141592654,F=30.0)REAL H(150),H2_R(128),H2_I(128),H_S(128)INTEGER K,NFFTCOMPLEX 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)THENK=K+1ENDIFNFFT=2**KDO 10,I=-50,50IF(I.EQ.0)THENH(I+51)=1.0/DT-2*FELSEH(I+51)=-SIN(2*PI*30.0*I*DT)/(PI*I*DT)END IF10 CONTINUEDO 20,I=102,12820 H(I)=0.0DO 30,I=1,12830 H2(I)=CMPLX(H(I),0.0)CALL FFT(128,H2,1.0)DO 40,I=1,128H2_R(I)=REAL(H2(I))H2_I(I)=AIMAG(H2(I))H_S(I)=(H2_R(I)**2+H2_I(I)**2)40 CONTINUEWRITE(2,*)(H_S(I),I=1,128)WRITE(1,*)(H(I),I=1,128)CLOSE(1)CLOSE(2)ENDCCCCCCCCCCCCC FFT 子程序CCCCCCCCCCCCCC SUBROUTINE FFT(LX,CX,SIGNI)COMPLEX CX(LX),CARG,CEXP,CW,CTEMPJ=1SC=1.0/LXIF(SIGNI.EQ.1.0)SC=1.0SIG=-SIGNIDO 30 I=1,LXIF(I.GT.J)GO TO 10CTEMP=CX(J)*SCCX(J)=CX(I)*SCCX(I)=CTEMP10 M=LX/220 IF(J.LE.M)GO TO 30J=J-MM=M/2IF(M.GE.1)GO TO 2030 J=J+ML=140 ISTEP=2*LDO 50 M=1,LCARG=(0.0,1.0)*(3.14159265*SIG*(M-1))/LCW=CEXP(CARG)DO 50 I=M,LX,ISTEPCTEMP=CW*CX(I+L)CX(I+L)=CX(I)-CTEMP50 CX(I)=CX(I)+CTEMPL=ISTEPIF(L.LT.LX)GO TO 40RETURNEND5.频谱分析CCCCCCCCCCCCC 频谱分析CCCCCCCCCCCCCCPROGRAM MAINPARAMETER (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 DFCOMPLEX XSIN_C(200),X_C(200),X_WA VE_C(200)CCCCCCCCCCCC 对SIN 的分析CCCCCCCCCCCCDO 100 I=1,128XSIN(I)=SIN(2*I*PI/64.0)100 CONTINUECCCCCCCCC 是数据转换成复数形式CCCCCCCCCCDO 200 I=1,128200 XSIN_C(I)=CMPLX(XSIN(I),0.0)CALL FFT(128,XSIN_C,1.0)DO 300 I=1,128300 X_SIN_R(I)=REAL(XSIN_C(I))DO 400 I=1,128400 X_SIN_I(I)=AIMAG(XSIN_C(I))OPEN(1,FILE='XSINSPRT.DAT',FORM='FORMATTED',STATUS='UNKNOWN') DO 1 I=1,1281 SINSPRT(I)=SQRT(X_SIN_R(I)**2+X_SIN_I(I)**2)WRITE(1,*)(SINSPRT(I), I=1,128)CLOSE(1)CCCCCCCCCCCCCC 对脉冲信号的分析CCCCCCCCCCCCCX(1)=1.0DO 800 I=2,128800 X(I)=0.0DO 900 I=1,128900 X_C(I)=CMPLX(X(I),0.0)CALL FFT(128,X_C,1.0)DO 1000 I=1,1281000 X_R(I)=REAL(X_C(I))DO 1100 I=1,1281100 X_I(I)=AIMAG(X_C(I))OPEN(5,FILE='XSPRT.DAT',FORM='FORMATTED',STATUS='UNKNOWN')DO 2 I=1,1282 XSPRT(I)=SQRT(X_R(I)**2+X_I(I)**2)WRITE(5,*)(XSPRT(I), I=1,128)CLOSE(5)CCCCCCCCCCCC 对地震波W A VE 的分析CCCCCCCCCCCOPEN(3,FILE='WA VE.DAT',FORM='FORMATTED',STATUS='UNKNOWN')READ(3,*)(WA VE(I), I=1,128)DO 500 I=1,128500 X_WA VE_C(I)=CMPLX(WA VE(I),0.0)CALL FFT(128,X_WA VE_C,1.0)DO 600 I=1,128600 X_WA VE_R(I)=REAL(X_WA VE_C(I))DO 700 I=1,128700 X_WA VE_I(I)=AIMAG(X_WA VE_C(I))OPEN(4,FILE='WA VESPRT.DAT',FORM='FORMATTED',STATUS='UNKNOWN') DO 3 I=1,1283 WA VESPRT(I)=SQRT(X_WA VE_R(I)**2+X_WA VE_I(I)**2)WRITE(4,*)(WA VESPRT(I), I=1,128)CLOSE(4)CLOSE(3)ENDCCCCCCCCCCCCCCC FFT 子程序CCCCCCCCCCCCCCCCCSUBROUTINE FFT(LX,CX,SIGNI)COMPLEX CX(LX),CARG,CEXP,CW,CTEMPJ=1SC=1.0/LXIF(SIGNI.EQ.1.0)SC=1.0SIG=-SIGNIDO 30 I=1,LXIF(I.GT.J)GO TO 10CTEMP=CX(J)*SCCX(J)=CX(I)*SCCX(I)=CTEMP10 M=LX/220 IF(J.LE.M)GO TO 30J=J-MM=M/2IF(M.GE.1)GO TO 2030 J=J+ML=140 ISTEP=2*LDO 50 M=1,LCARG=(0.0,1.0)*(3.14159265*SIG*(M-1))/LCW=CEXP(CARG)DO 50 I=M,LX,ISTEPCTEMP=CW*CX(I+L)CX(I+L)=CX(I)-CTEMP50 CX(I)=CX(I)+CTEMPL=ISTEPIF(L.LT.LX)GO TO 40RETURNEND6.分析补零对振幅谱的影响CCCCCCCCCCC 分析补零对振幅谱的影响CCCCCCCCCCCPROGRAM MAINCCCCCCCCCCC 不补零时的振幅普CCCCCCCCCCCCCCCCC CCCCCCCCCCC 对函数SIN 的分析CCCCCCCCCCCCCCCCPARAMETER (PI=3.1415)REAL XSIN(150),XWA VE(150)REAL XWA VE_R(150),XWA VE_I(150)REAL XSIN_Z(150),XWA VE_Z(150)REAL XSIN_R(150),XSIN_I(150)COMPLEX XSIN_60C(150),XWA VE_60C(150)COMPLEX XSIN_64C(150),XWA VE_64C(150)COMPLEX XSIN_128C(150),XWA VE_128C(150)REAL XSIN_60S(150),XWA VE_60ZS(150)REAL XSIN_64S(150),XWA VE_64ZS(150)REAL XSIN_128S(150),XWA VE_128ZS(150)OPEN(1,FILE='XSIN_60S.DA T',FORM='FORMATTED',STATUS='UNKNOWN') DO 1 I=1,601 XSIN(I)=SIN(2*PI*I/30.0)DO 11 I=1,6011 XSIN_60C(I)=CMPLX(XSIN(I),0.0)CALL DFT(60,XSIN_60C,1.0)DO 12 I=1,6012 XSIN_R(I)=REAL(XSIN_60C(I))DO 13 I=1,6013 XSIN_I(I)=AIMAG(XSIN_60C(I))DO 14 I=1,6014 XSIN_60S(I)=SQRT(XSIN_R(I)**2+XSIN_I(I)**2)WRITE(1,*)(XSIN_60S(I),I=1,60)CLOSE(1)CCCCCCCCCCCCC 对地震数据的分析CCCCCCCCCCCCCCOPEN(2,FILE='WA VE.DAT',FORM='FORMATTED',STATUS='UNKNOWN')READ(2,*)(XWA VE(I),I=1,60)CLOSE(2)DO 21 I=1,6021 XWA VE_60C(I)=CMPLX(XWA VE(I),0.0)CALL DFT(60,XWA VE_60C,1.0)DO 22 I=1,6022 XWA VE_R(I)=REAL(XWA VE_60C(I))DO 23 I=1,6023 XWA VE_I(I)=AIMAG(XWA VE_60C(I))DO 24 I=1,6024 XWA VE_60ZS(I)=SQRT(XWA VE_R(I)**2+XWA VE_I(I)**2)OPEN(3,FILE='XWA VE_60ZS.DA T',FORM='FORMATTED',STA TUS='UNKNOWN') WRITE(3,*)(XWA VE_60ZS(I),I=1,60)CLOSE(3)CCCCCCCCCCCCC 对补零数据的分析CCCCCCCCCCCCCCCC CCCCCCCCCCCC 对函数SIN 的分析CCCCCCCCCCCCCCCCOPEN(4,FILE='XSIN_64S.DA T',FORM='FORMATTED',STATUS='UNKNOWN')DO 31 I=61,6431 XSIN(I)=0.0DO 32 I=1,6432 XSIN_64C(I)=CMPLX(XSIN(I),0.0)CALL DFT(64,XSIN_64C,1.0)DO 33 I=1,6433 XSIN_R(I)=REAL(XSIN_64C(I))DO 34 I=1,6434 XSIN_I(I)=AIMAG(XSIN_64C(I))DO 35 I=1,6435 XSIN_64S(I)=SQRT(XSIN_R(I)**2+XSIN_I(I)**2)WRITE(4,*)(XSIN_64S(I),I=1,64)CLOSE(4)OPEN(5,FILE='XSIN_128S.DA T',FORM='FORMATTED',STATUS='UNKNOWN')DO 41 I=65,12841 XSIN(I)=0.0DO 42 I=1,12842 XSIN_128C(I)=CMPLX(XSIN(I),0.0)CALL DFT(128,XSIN_128C,1.0)DO 43 I=1,12843 XSIN_R(I)=REAL(XSIN_128C(I))DO 44 I=1,12844 XSIN_I(I)=AIMAG(XSIN_128C(I))DO 45 I=1,12845 XSIN_128S(I)=SQRT(XSIN_R(I)**2+XSIN_I(I)**2)WRITE(5,*)(XSIN_128S(I),I=1,128)CLOSE(5)CCCCCCCCCCCCC 对地震数据的分析CCCCCCCCCCCCCCOPEN(6,FILE='XWA VE_64ZS.DAT',FORM='FORMATTED',STATUS='UNKNOWN') DO 51 I=61,6451 XWA VE(I)=0.0DO 52 I=1,6452 XWA VE_64C(I)=CMPLX(XWA VE(I),0.0)CALL DFT(64,XWA VE_64C,1.0)DO 53 I=1,6453 XWA VE_R(I)=REAL(XWA VE_64C(I))DO 54 I=1,6454 XWA VE_I(I)=AIMAG(XWA VE_64C(I))DO 55 I=1,6455 XWA VE_64ZS(I)=SQRT(XWA VE_R(I)**2+XWA VE_I(I)**2)WRITE(6,*)(XWA VE_64ZS(I),I=1,64)CLOSE(6)OPEN(7,FILE='XWA VE_128ZS.DAT',FORM='FORMATTED',STATUS='UNKNOWN') DO 61 I=65,12861 XWA VE(I)=0.0DO 62 I=1,12862 XWA VE_128C(I)=CMPLX(XWA VE(I),0.0)CALL DFT(128,XWA VE_128C,1.0)DO 63 I=1,12863 XWA VE_R(I)=REAL(XWA VE_128C(I))DO 64 I=1,12864 XWA VE_I(I)=AIMAG(XWA VE_128C(I))DO 65 I=1,12865 XWA VE_128ZS(I)=SQRT(XWA VE_R(I)**2+XWA VE_I(I)**2)WRITE(7,*)(XWA VE_128ZS(I),I=1,128)CLOSE(7)ENDCCCCCCCCCCCCCCC DFT 子程序CCCCCCCCCCCCCCCCSUBROUTINE DFT(N,C,S)COMPLEX C(N),Y(512),WIF(S.EQ.1.0)Z=1.0IF(S.EQ.-1.0)Z=1.0/FLOAT(N)DO 20 I=1,NY(I)=(0.0,0.0)DO 20 J=1,NA=COS(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))B=SIN(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))*(-S)W=CMPLX(A,B)20 Y(I)=Y(I)+C(J)*WDO 30 I=1,N30 C(I)=Y(I)*ZRETURNEND7.线性褶积与循环褶积CCCCCCCCC 线性滤波与循环褶积比较CCCCCCCCCPROGRAM MAINCCCCCCCCCCCCC 线性滤波CCCCCCCCCCCCCCCCCCCCCCCCCCC H是低通滤波因CCCCCCCCCCCCPARAMETER (PI=3.141592654, DT=0.002,F=70.0)REAL X(100),X2(101),H(-50:50),Y_LINEAR(200),Y_CIRCLE(200)OPEN(1,FILE='INPUT1.DA T',FORM='FORMATTED',STATUS='UNKNOWN')READ(1,*)(X(I),I=1,100)CCCCCCCCCCCCCCC低通滤波器CCCCCCCCCCCCCDO 10 I=-50,50IF (I.EQ.0)THENH(I)=2*F*PI/PIELSEH(I)=SIN(2*PI*F*I*DT)/(PI*I*DT)END IF10 CONTINUECALL CON(X,H,Y_LINEAR,100,101,200)CCCCCCCCCCCC输出滤波后的数据CCCCCCCCCCCCCCCOPEN(2,FILE='Y_LINEAR.DA T',FORM='FORMATTED',STATUS='UNKNOWN')WRITE(2,*)(Y_LINEAR(I),I=1,200)CLOSE(2)CCCCCCCCCCCCC 圆周滤波的设计CCCCCCCCCCCCOPEN(3,FILE='Y_CIRCLE.DAT',FORM='FORMATTED',STATUS='UNKNOWN')DO 20 I=1,10020 X2(I)=X(I)X2(100)=0.0CALL CIRCON(X2,H,Y_CIRCLE,101)WRITE(3,*)(Y_CIRCLE(I),I=1,101)CLOSE(3)CCCCCCC 线性卷积输出点等于圆周卷积CCCCCCOPEN(4,FILE='Y_LINEARCUT.DAT',FORM='FORMATTED',STATUS='UNKNOWN') WRITE(4,*)(Y_LINEAR(I),I=1,101)CLOSE(4)ENDCCCCCCCCCCCCCC 线性卷积子程序CCCCCCCCCCCCCCSUBROUTINE CON(A,B,C,I,J,K)DIMENSION A(I),B(J),C(K)DO 1 K1=1,K1 C(K1)=0.0DO 2 I1=1,IDO 2 I2=1,JII=I1+I2-12 C(II)=C(II)+A(I1)*B(I2)*0.004RETURNENDCCCCCCCCCCCCCCC 圆周卷积子程序CCCCCCCCCCCCCCSUBROUTINE CIRCON(X,H,Y,N)REAL X(N),H(N),Y(N)DO I=1,NY(I)=0.0ENDDODO I=1,NDO J=1,NIF(I-J.GT.0) Y(I)=Y(I)+X(J)*H(I-J)ENDDOENDDOEND8.最小平方反滤波CCCCCCCCCC 两种反射系数序列的对比CCCCCCCCCCPROGRAM MAINREAL B(12),R(200),X_CONV(211),RXX(211)REAL A(211),A_CONV(270)DATA B/17.5,12.5,7.0,0.0,-3.5,-7.0,-3.0,-1.0,0.0,1.4,3.5,2.0/OPEN(1,FILE='INPUT8.DAT',FORM='FORMATTED',STATUS='UNKNOWN') READ(1,*)(R(I),I=1,200)CCCCCCCCCCCCCC 求合成地震记录CCCCCCCCCCCCCCCALL CON(R,B,X_CONV,200,12,211)CCCCCCCCCCCC 求地震子波的自相关CCCCCCCCCCCCCALL AUTCOR(211,211,X_CONV,RXX)CCCCCCCCCCCCCCC 求反滤波因子CCCCCCCCCCCCCCCCALL PEO(211,RXX,A)OPEN(2,FILE='A.DAT',FORM='FORMA TTED',STATUS='UNKNOWN')WRITE(2,*)(A(I),I=1,60)CCCCCCC 通过反滤波因子求反射系数序列CCCCCCCCALL CON(A,X_CONV,A_CONV,60,211,270)OPEN(3,FILE='A_CONV.DA T',FORM='FORMATTED',STATUS='UNKNOWN') WRITE(3,*)(A_CONV(I),I=1,270)CLOSE(3)CLOSE(2)CLOSE(1)ENDCCCCCCCCCCCCCCC 褶积子程序CCCCCCCCCCCCCCCC SUBROUTINE CON(A,B,C,I,J,K)DIMENSION A(I),B(J),C(K)DO 1 K1=1,K1 C(K1)=0.0DO 2 I1=1,IDO 2 I2=1,JII=I1+I2-12 C(II)=C(II)+A(I1)*B(I2)*0.004RETURNENDCCCCCCCCCCCCCCC 自相关子程序CCCCCCCCCCCCCCC SUBROUTINE AUTCOR(J,K,C,A)DIMENSION C(J),A(K)DO 7 N=1,KA(N)=0.0I=NDO 6 IN=I,JIL=IN-N+16 A(N)=A(N)+C(IN)*C(IL)7 CONTINUERETURNENDCCCCCCCCCCCCCC 求反滤波因子CCCCCCCCCCCCCCCC SUBROUTINE PEO(LR,R,A)C PEO IS THE AUXILIARY TOEPLITZ RECURSIONC IT GIVES THE PREDICTION-ERROR OPERA TOR A(I)DIMENSION R(LR),A(LR)V=R(1)A(1)=1.0IF(LR.EQ.1)RETURNDO 6 L=2,LRD=0.0L3=L-1DO 1 J=1,L3K=L-J+11 D=D+A(J)*R(K)C=D/VIF(L.EQ.2)GOTO 4L1=(L-2)/2L2=L1+1IF(L2.LT.2)GOTO 3DO 2 J=2,L2HOLD=A(J)K=L-J+1A(J)=A(J)-C*A(K)2 A(K)=A(K)-C*HOLDIF(2*L1.EQ.L-2)GOTO 43 LT3=L2+1A(LT3)=A(LT3)-C*A(LT3)4 A(L)=-C6 V=V-C*DRETURNEND9.零相位转换CCCCCCCCCCC 零相位的转换CCCCCCCCCCCCCPROGRAM MAINREAL B(32),B1(32),BC_S(32),B_CONV(32)REAL BC_R(32),BC_I(32)DATA B/0,1,3,5,7,9,10,8,6,4,2,0,-1,0,1,2,3,1.5,0,-1,0,2,1,0,0/COMPLEX B_C(32),BC_SZ(32)DO 10,I=1,32IF(I.GE.26)THENB_C(I)=(0.0,0.0)ELSEB_C(I)=CMPLX(B(I),0.0)END IF10 CONTINUECALL FFT(32,B_C,1.0)DO 20,I=1,32BC_R(I)=REAL(B_C(I))BC_I(I)=AIMAG(B_C(I))BC_S(I)=SQRT(BC_R(I)**2+BC_I(I)**2)20 BC_SZ(I)=CMPLX(BC_S(I),0.0)CALL FFT(32,BC_SZ,-1.0)DO 30,I=1,3230 B_CONV(I)=REAL(BC_SZ(I))OPEN (1,FILE='B_CONV.DA T',FORM='FORMATTED',STATUS='UNKNOWN') WRITE(1,*)(B_CONV(I),I=1,32)CLOSE(1)ENDCCCCCCCCCCCCCC FFT子程序CCCCCCCCCCCCCCSUBROUTINE FFT(LX,CX,SIGNI)COMPLEX CX(LX),CARG,CEXP,CW,CTEMPJ=1SC=1.0/LXIF(SIGNI.EQ.1.0)SC=1.0SIG=-SIGNIDO 30 I=1,LXIF(I.GT.J)GO TO 10CTEMP=CX(J)*SCCX(J)=CX(I)*SCCX(I)=CTEMP10 M=LX/220 IF(J.LE.M)GO TO 30J=J-MM=M/2IF(M.GE.1)GO TO 2030 J=J+ML=140 ISTEP=2*LDO 50 M=1,LCARG=(0.0,1.0)*(3.14159265*SIG*(M-1))/LCW=CEXP(CARG)DO 50 I=M,LX,ISTEPCTEMP=CW*CX(I+L)CX(I+L)=CX(I)-CTEMP50 CX(I)=CX(I)+CTEMPL=ISTEPIF(L.LT.LX)GO TO 40RETURNEND10.最小相位转换CCCCCCCCCCCC 最小相位转换CCCCCCCCCCCCPROGRAM MAIMREAL B(25),B0(75)DATA B/0,1,2,3,1.5,0,-1,0,1,3,5,7,9,10,8,6,4,2,0,-1,0,2,3,2,0/REAL RBB(25),Y0(25),A0(25)REAL SUMX,C(49),G(49)OPEN(1,FILE='B0.DA T',FORM='FORMATTED',STA TUS='UNKNOWN') OPEN(2,FILE='B.DAT',FORM='FORMATTED',STATUS='UNKNOWN') WRITE(2,*)(B(I),I=1,25)CALL AUTCOR(25,25,B,RBB)CALL PEO(25,RBB,Y0)SUMX=0.0DO 10,I=1,25CALL CON(B,Y0,C,25,25,49)10 SUMX=SUMX+(C(I))**2B0(1)=1.0/SQRT(SUMX)20 A0(I)=B0(1)*Y0(I)CALL CON(B,A0,G,25,25,49)CALL COAUTCOR(25,49,73,B,G,B0)WRITE(1,*)(B0(I),I=1,73)CLOSE(2)CLOSE(1)ENDCCCCCCCCCCCC 求反滤波因子CCCCCCCCCCCC SUBROUTINE PEO(LR,R,A)C PEO IS THE AUXILIARY TOEPLITZ RECURSIONC IT GIVES THE PREDICTION-ERROR OPERA TOR A(I)DIMENSION R(LR),A(LR)V=R(1)A(1)=1.0IF(LR.EQ.1)RETURNDO 6 L=2,LRD=0.0L3=L-1DO 1 J=1,L3K=L-J+11 D=D+A(J)*R(K)C=D/VIF(L.EQ.2)GOTO 4L1=(L-2)/2L2=L1+1IF(L2.LT.2)GOTO 3DO 2 J=2,L2HOLD=A(J)K=L-J+1A(J)=A(J)-C*A(K)2 A(K)=A(K)-C*HOLDIF(2*L1.EQ.L-2)GOTO 43 LT3=L2+1A(LT3)=A(LT3)-C*A(LT3)4 A(L)=-C6 V=V-C*DRETURNENDCCCCCCCCCCCCCC 自相关子程序CCCCCCCCCCCCCCC SUBROUTINE AUTCOR(J,K,C,A)DIMENSION C(J),A(K)A(N)=0.0I=NDO 6 IN=I,JIL=IN-N+16 A(N)=A(N)+C(IN)*C(IL)7 CONTINUERETURNENDCCCCCCCCCCCCCC 互相关子程序CCCCCCCCCCCCCCCSUBROUTINE COAUTCOR(J,L,K,A,B,C)DIMENSION A(J),B(L),C(K)DO 7 N=1,KC(N)=0.0IF (J.LE.L+1-N) THENMIN=JELSEMIN=L+1-NENDIFDO 6 IA=1,MINIB=IA+N-16 C(N)=C(N)+A(IA)*B(IB)7 CONTINUERETURNENDCCCCCCCCCCCCCCC 卷积子程序CCCCCCCCCCCCCCCCSUBROUTINE CON(A,B,C,I,J,K)DIMENSION A(I),B(J),C(K)DO 1 K1=1,K1 C(K1)=0.0DO 2 I1=1,IDO 2 I2=1,JII=I1+I2-12 C(II)=C(II)+A(I1)*B(I2)*0.004RETURNEND11.静校正CCCCCCCCCCCCCC 静校正程序设计CCCCCCCCCCCCCCCPROGRAM MAININTEGER I,J,NREAL X(10,100),X_SAMPLE(100),X_CT(0:19),X_COR(-19:19),XDT(10),TH(100) REAL MAXOPEN(1,FILE='SEIS.DA T',FORM='FORMATTED',STATUS='UNKNOWN') OPEN(2,FILE='SEIS11.DA T',FORM='FORMATTED',STATUS='UNKNOWN') DO 100 I=1,10READ(1,*)(X(I,J),J=1,100)100 CONTINUEPRINT*,'PLEASE INPUT SAMPLE LINE NUMBER !'READ*,NDO 200 J=1,100X_SAMPLE(J)=X(N,J)200 CONTINUEDO 300 I=1,10DO 400 J=1,100TH(J)=X(I,J)400 CONTINUECALL COAUTCOR(100,100,20,TH,X_SAMPLE,X_CT)DO J=0,19X_COR(J)=X_CT(J)ENDDOCALL COAUTCOR(100,100,20,X_SAMPLE,TH,X_CT)DO J=0,19X_COR(-J)=X_CT(J)ENDDOMAX=X_COR(-19)DO J=-18,19IF(MAX.LT.X_COR(J))THENMAX=X_COR(J)XDT(I)=JENDIFENDDOIF(XDT(I).GT.0)THENDO J=100-XDT(I),1,-1X(I,J+XDT(I))=X(I,J)ENDDODO J=1,XDT(I)-1X(I,J)=0ENDDOELSEIF(XDT(I).LT.0)THENDO J=1,100+XDT(I)X(I,J)=X(I,J-XDT(I))ENDDODO J=101+XDT(I),100X(I,J)=0ENDDOENDIF300 CONTINUEWRITE(2,*)((X(I,J),J=1,100),I=1,10)CLOSE(2)CLOSE(1)ENDCCCCCCCCCCCCC 互相关子程序CCCCCCCCCCCCC SUBROUTINE COAUTCOR(J,L,K,A,B,C)DIMENSION A(J),B(L),C(K)DO 7 N=1,KC(N)=0.0IF (J.LE.L+1-N) THENMIN=JELSEMIN=L+1-NENDIFDO 6 IA=1,MINIB=IA+N-16 C(N)=C(N)+A(IA)*B(IB)7 CONTINUERETURNEND【附录2：有关图件】。