电阻率测深一维模拟退火反演反演

操作方法

• 01

  设置理论模型。

  

• 02

  模拟退火程序。 PROGRAM MIAN IMPLICIT NONE INTEGER  :: N0,M,N,I,J,TT,NN,COUNT,GN REAL(8)  :: TM,DX REAL(8)  :: H0(20),H(20),P0(20),P(20) REAL(8)  :: REX(50),REX0(50),CX REAL(8)  :: DC,XK(500),AP0(500),AP(500) REAL(8)  :: XX(1000),YY(1000),GX(1000) REAL(8)  :: MO,T,T0,RE0(50),E,EE,MT,Y REAL(8)  :: AAP0(500),DE,U INTEGER  :: IN,ACCEPTPOINTS CHARACTER(20) :: SOUNDING_DATA_NAME,THEORETICAL_MODEL_NAME REAL(8),DIMENSION(:),ALLOCATABLE   :: B,W REAL(8),DIMENSION(:,:),ALLOCATABLE :: A,V OPEN(20,FILE='INVERSION_RESULT_DATA.TXT') OPEN(30,FILE='RESULT COMPARSION.TXT') T0=500. TT=10000 WRITE(*,*) '请选择反演模式!(输入0为模型反演，输入1为实测数据反演。)' READ(*,*) DC ! 小于等于0，模型反演，大于0实测数据反演 ! ------------读入初始模型参数-------------- WRITE(*,*) '请输入初始模型层数' READ(*,*) N0  ! 模型层数 CALL RANDOM_SEED() DO I=1,N0-1 CALL RANDOM_NUMBER(MO) H0(I)=100.*MO P0(I)=1000.*MO END DO CALL RANDOM_SEED() CALL RANDOM_NUMBER(MO) P0(N0)=1000.*MO ! ------------------------------------------- IF(DC) 10,10,20 20  WRITE(*,*)'请输入野外数据文件名！' READ(*,*) SOUNDING_DATA_NAME OPEN(40,FILE=SOUNDING_DATA_NAME) M=0 N=N0 DO READ(40,*,END=2) GX(M+1),AP(M+1) M=M+1 END DO 2   CONTINUE ! 当采集的数据量不够的时候进行样条插值 CX=GX(M) GN=M IF(M.LT.N*2-1) THEN M=INT(N*1.03) DX=CX/(M-1)     ! 插值后的X间距。 CALL SP(GX,AP,GN,CX,M,XX,YY) DO I=1,M AP(I)=YY(I) END DO DO I=1,M-1 XK(I)=DX*I END DO PRINT*,'插值后的采样点数',M ELSE DO I=1,M XK(I)=GX(I) END DO END IF GOTO 101 ! 如果是理论模型反演则读入理论模型 10  WRITE(*,*)'请输理论模型反演参数文件名！' READ(*,*) THEORETICAL_MODEL_NAME OPEN(50,FILE=THEORETICAL_MODEL_NAME) WRITE(*,*)'请输入采样数据个数' READ(*,*) M READ(50,*) N  ! 模型层数 DO I=1,N-1 READ(50,*) H(I),P(I) END DO READ(50,*) P(N) DO I=1,M XK(I)=10.**(I/6.) END DO WRITE(30,25) (H(I),I=1,N-1),(P(I),I=1,N) 25  FORMAT('理论模型:',/(10X,2F15.5/(10X,3F15.5))) CALL RS(M,N,XK,H,P,AP) CALL PLOT(M,N,XK,H,P,'THEORETICAL_PLOT.TXT') OPEN(115,FILE='THEORETICAL_FORWARD_DATA.TXT') DO I=1,M WRITE(115,*) XK(I),AP(I) END DO 101 NN=N*2-1 CALL RS(M,N,XK,H0,P0,AP0) IN=0 T=T0 CALL COUQ(N,H0,P0,RE0) ! ------------------------------------------------------------- ! 模拟退火部分 call random_seed () DO WHILE(IN.LT.TT) CALL DEUQ(N,H0,P0,RE0) CALL RS(M,N,XK,H0,P0,AP0) CALL INTE(M,AP,AP0,E) DO I=1,NN REX0(I)=RE0(I) END DO DO I=1,NN call random_number(MT) call RAND(MT,T,Y) REX0(I)=RE0(I)+Y*126. IF(REX0(I).LT.0.) REX0(I)=0. IF(REX0(I).GT.10000.) REX0(I)=10000. END DO CALL DEUQ(N,H0,P0,REX0) CALL RS(M,N,XK,H0,P0,AAP0) CALL INTE(M,AP,AAP0,EE) DE=EE-E call random_number (U) IN=IN+1 T=T*0.95 IF(DE.LE.0.)         THEN              ! 表达移动后得到更优解,则总是接受移动 DO I=1,NN RE0(I)=REX0(I) END DO AcceptPoints=AcceptPoints+1 ELSEIF(exp(-DE/T)>U) THEN DO I=1,NN RE0(I)=REX0(I) END DO AcceptPoints=AcceptPoints+1 END IF END DO PRINT*,IN,EE,E PRINT*,'AcceptPoints:',AcceptPoints WRITE(30,15) (H0(I),I=1,N-1),(P0(I),I=1,N) 15  FORMAT('理论模型:',/(10X,2F15.5/(10X,3F15.5))) CLOSE(30) CALL DEUQ(N,H0,P0,RE0) CALL PLOT(M,N,XK,H0,P0,'INVERSION_PLOT.TXT') CALL RS(M,N,XK,H0,P0,AP0) DO I=1,M WRITE(20,*) XK(I), AP0(I) END DO CLOSE(20) END ! =================================================================== ! 求目标函数 SUBROUTINE INTE(M,GM,DZ,E) INTEGER :: I,M REAL(8) :: E,GM(M),DZ(M) E=0. DO I=1,M E=E+(GM(I)-DZ(I))**2 END DO END ! =================================================================== !产生随机数及作用于符号函数SGN SUBROUTINE RAND(U,T,Y) IMPLICIT NONE REAL(8) :: U,SGNN,T,Y IF((U-0.5).LT.0) THEN SGNN=-1. ELSEIF((U-0.5).EQ.0) THEN SGNN=0. ELSE SGNN=1. END IF Y=T*SGNN*((1.+1./T)**ABS(2.*U-1.)-1.) END ! ================================================================ ! 把所有的未知量放在一个数组里面。 SUBROUTINE COUQ(N,H,P,RE) IMPLICIT NONE INTEGER :: I,N REAL(8) :: H(N),P(N),RE(N*2-1) DO I=1,N RE(I)=P(I) END DO DO I=1,N-1 RE(N+I)=H(I) END DO RETURN END SUBROUTINE COUQ ! ================================================================ ! 把总数组里面的值分别输出到各个未知变量或初始变量中。 SUBROUTINE DEUQ(N,H,P,RE) IMPLICIT NONE INTEGER :: I,N REAL(8) :: H(N),P(N),RE(N*2-1) DO I=1,N P(I)=RE(I) END DO DO I=1,N-1 H(I)=RE(N+I) END DO RETURN END SUBROUTINE DEUQ ! ================================================================ ! ================================================================ SUBROUTINE RS(M,N,R,H,P,PS) IMPLICIT NONE INTEGER:: K,M,N,I,J REAL(8):: C(20),PS(M),P(N),R(M),H(N),T(N),MM(20) DATA C/0.003042,-0.001198,0.01284,0.0235,0.08688,0.2374,0.6194, & 1.1817,0.4248,-3.4507, 2.7044,-1.1324,0.393,-0.1436,0.05812, & -0.02521,0.01125,-0.004978,0.002072,-0.000318/ T(N)=P(N) DO I=1,M PS(I)=0. DO K=1,20 MM(K)=(0.11396*10**(K/6.))/(R(I)) DO J=N-1,1,-1 T(J)=P(J)*(P(J)*TANH(MM(K)*H(J))+T(J+1))/(P(J)+ & T(J+1)*TANH(MM(K)*H(J))) END DO PS(I)=PS(I)+T(1)*C(K) END DO END DO RETURN END SUBROUTINE ! ================================================================ ! 画出模型曲线或反演曲线 SUBROUTINE PLOT(M,N,XK,H,P,PLOTNAME) IMPLICIT NONE INTEGER       :: M,N,I REAL(8)       :: DEP(N+1),H(N),P(N),XK(M) CHARACTER(20) :: PLOTNAME OPEN(10,FILE=PLOTNAME) H(N)=XK(M) DEP(1)=0. DO I=2,N+1 DEP(I)=DEP(I-1)+H(I-1) END DO WRITE(10,*) 0.1,P(1) WRITE(10,*) H(1),P(1) DO I=2,N WRITE(10,*) DEP(I),P(I) WRITE(10,*) DEP(I+1),P(I) END DO END SUBROUTINE PLOT ! ================================================================ ! 调用三次样条插值 SUBROUTINE SP(X,Y,N,CX,M,XX,YY) IMPLICIT NONE INTEGER :: I,J,N,M REAL(8) :: X(N),Y(N),YY(M),XX(M),D1,DN REAL(8) :: DX,CX D1=0.           ! 左端点的一阶导数 DN=0.           ! 右端点的一阶导数 DX=CX/(M-1) DO J=1,M XX(J)=DX*(J-1) END DO CALL SPL1(X,Y,N,D1,DN,M,XX,YY) END SUBROUTINE SP ! ====================================================== ! 三次样条插值子程序 SUBROUTINE SPL1(X,Y,N,DY1,DYN,M,XX,S) IMPLICIT NONE INTEGER :: I,J,N,M REAL(8) :: X(N),Y(N),XX(M),S(M),H(N),DY(N) REAL(8) :: DY1,DYN,H0,H1,BETA,ALPHA DY(1)=0.0 H(1)=DY1 H0=X(2)-X(1) DO J=2,N-1 H1=X(J+1)-X(J) ALPHA=H0/(H0+H1) BETA=(1.0-ALPHA)*(Y(J)-Y(J-1))/H0 BETA=3.0*(BETA+ALPHA*(Y(J+1)-Y(J))/H1) DY(J)=-ALPHA/(2.0+(1.0-ALPHA)*DY(J-1)) H(J)=(BETA-(1.0-ALPHA)*H(J-1)) H(J)=H(J)/(2.0+(1.0-ALPHA)*DY(J-1)) H0=H1 END DO DY(N)=DYN DO J=N-1,1,-1 DY(J)=DY(J)*DY(J+1)+H(J) END DO DO J=1,N-1 H(J)=X(J+1)-X(J) END DO H1=H(N-1)*H(N-1) DO J=1,M IF (XX(J).GE.X(N)) THEN I=N-1 ELSE I=1 60          IF (XX(J).GT.X(I+1)) THEN I=I+1 GOTO 60 END IF END IF H1=(X(I+1)-XX(J))/H(I) S(J)=(3.0*H1*H1-2.0*H1*H1*H1)*Y(I) S(J)=S(J)+H(I)*(H1*H1-H1*H1*H1)*DY(I) H1=(XX(J)-X(I))/H(I) S(J)=S(J)+(3.0*H1*H1-2.0*H1*H1*H1)*Y(I+1) S(J)=S(J)-H(I)*(H1*H1-H1*H1*H1)*DY(I+1) END DO RETURN END SUBROUTINE SPL1 ! ==============================================================

• 03

  反演结果对比。