地球科学学院
地球物理反演理论与方法课程作业
课程名称:地球物理反演理论与方法
指导老师:
学生姓名:
学号:
1 电阻率测深一维最小二乘最优化反演
本次反演的理论模型层数为三层,反演的曲线为H型。反演出的曲线与理论曲线画Grapher成图后肉眼无法分别,故图1.1只画了初始模型和理论模型曲线。
表 1.1 反演初始模型设置
第一层第二层第三层厚度15 35
电阻率200 60 1000
表 1.2 反演理论模型设置
第一层第二层第三层厚度10 20
电阻率250 40 600
表1.3 反演结果
第一层第二层第三层厚度10.00000 20.00001
电阻率250.00000 40.00002 599.99999
图1.1 初始模型及理论结果正演曲线图1.2 反演结果
2 电阻率测深一维模拟退火反演反演
根据模拟退火的特点,初始模模型为随机设置,理论模型同上。
模拟退火的容易编程实现,前人的研究也取得了一定的成果。但是模拟退火反演结果的好坏与初始温度及降火方案有很大的联系,由于时间的关系这里没有作过多的研究,因此反演的结果一般。
表 2.1 反演理论模型设置
第一层第二层第三层厚度10 20
电阻率250 40 600
表2.2 反演结果
第一层第二层第三层厚度7.78051 42.53932
电阻率252.79179 80.50919 600.43488
图2.1 正演曲线及反演结果
! ------------------------------------------------------------------------------------------------------ ! ================电阻率测深一维最小二乘最优化反演================ ! --------------------------------------HAIBIN 2012.06.28---------------------------------------- ! -----------------------------------------------主程序----------------------------------------------- PROGRAM MIAN
IMPLICIT NONE
INTEGER :: N0,M,N,I,J,TT,NN,COUNT,GN
REAL(8) :: TM,ERROR,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)
CHARACTER(20) :: SOUNDING_DATA_NAME,THEORETICAL_MODEL_NAME REAL(8),DIMENSION(:),ALLOCATABLE :: B,W
REAL(8),DIMENSION(:,:),ALLOCATABLE :: A,V
OPEN(10,FILE='ORIGINAL_MODEL.TXT')
OPEN(20,FILE='INVERSION_RESULT_DATA.TXT')
OPEN(30,FILE='RESULT COMPARSION.TXT')
OPEN(110,FILE='ORIGINAL_FORWARD_DA TA.TXT')
TT=1000
ERROR=1.0D-5
WRITE(*,*) '请选择反演模式!(输入0为模型反演,输入1为实测数据反演。)' READ(*,*) DC ! 小于等于0,模型反演,大于0实测数据反演
! ------------读入初始模型参数--------------
READ(10,*) N0 ! 模型层数
DO I=1,N0-1
READ(10,*) H0(I),P0(I)
END DO
READ(10,*) P0(N0)
WRITE(30,5) (H0(I),I=1,N0-1),(P0(I),I=1,N0)
5 FORMA T('初始模型:',/(10X,2F15.5/(10X,3F15.5)))
! -------------------------------------------
IF(DC) 10,10,20
20 WRITE(*,*)'请输入野外数据文件名!'
READ(*,*) SOUNDING_DA TA_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 FORMA T('理论模型:',/(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 CALL RS(M,N,XK,H0,P0,AP0)
CALL PLOT(M,N,XK,H0,P0,'ORIGINAL_PLOT.TXT')
DO I=1,M
WRITE(110,*) XK(I),AP0(I)
END DO
NN=N*2-1
ALLOCATE(A(M,NN))
ALLOCATE(B(M))
ALLOCATE(W(NN))
ALLOCATE(V(NN,NN))
CALL BI(AP0,AP,TM,M,B)
DO WHILE (TM.GT.ERROR)
COUNT=COUNT+1
CALL COUQ(N,H0,P0,REX0)
CALL AIJ(M,NN,XK,REX0,AP0,A)
CALL SVDCMP(A,M,NN,W,V)
CALL WIGGINSMETHOD(NN,W)
CALL SVBKSB(A,W,V,M,NN,B,REX)
CALL DM(REX0,REX,NN)
CALL DEUQ(N,H0,P0,REX0)
CALL RS(M,N,XK,H0,P0,AP0)
CALL BI(AP0,AP,TM,M,B)
IF(COUNT.GE.TT) EXIT
END DO
PRINT*,'均方误差:',TM
CALL PLOT(M,N,XK,H0,P0,'INVERSION_PLOT.TXT')
WRITE(30,15) (H0(I),I=1,N0-1),(P0(I),I=1,N0)
15 FORMA T('反演结果:',/(10X,2F15.5/(10X,3F15.5)))
DO I=1,M
WRITE(20,*) XK(I), AP0(I)
END DO
CLOSE(20)
WRITE(*,*)'共反演',COUNT,'次'
END
! ================================================================ ! 画出模型曲线或反演曲线
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 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
! ================================================================ ! 为了防止模型参量修改过量所做的规定。
! XX0----输入初始参量,返回修改后的参量。
! XX ----线性方程组的解,既修改量。
SUBROUTINE DM(XX0,XX,N)
IMPLICIT NONE
INTEGER :: I, N
REAL(8) :: XX0(N), XX(N)
DO I=1,N
IF(XX(I).GT.0.5) THEN
XX(I)=0.5
ELSEIF(XX(I).LT.-0.5) THEN
XX(I)=-0.5
END IF
END DO
DO I=1,N
XX0(I)=(1+XX(I))*XX0(I)
END DO
RETURN
END SUBROUTINE DM
! ================================================================ 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)
DA TA 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
! ================================================================ ! 求理论值对各参量的偏导数,形成雅克比矩阵A。
! 模型参数的无量纲化处理。
! 见《地球物理数据处理基础》107页
SUBROUTINE AIJ(M,N,R,RE,AP0,A)
IMPLICIT NONE
INTEGER :: M, N,I,J,NP
REAL(8) :: AP(M),H(N),P(N),RE(N)
REAL(8) :: AP0(M),A(M,N),R(M)
NP=(N+1)/2
DO I=1,N
RE(I)=1.1*RE(I)
CALL DEUQ(NP,H,P,RE)
CALL RS(M,NP,R,H,P,AP)
CALL COUQ(NP,H,P,RE)
RE(I)=RE(I)/1.1
DO J=1,M
A(J,I)=10.*(AP(J)/AP0(J)-1.0)
END DO
END DO
RETURN
END SUBROUTINE AIJ
! ================================================================
! 形成方程右端.
! 见《地球物理数据处理基础》107页.
SUBROUTINE BI(DZ,GM,TM,M,B)
IMPLICIT NONE
INTEGER :: K, M
REAL(8) :: B(M), GM(M),DZ(M)
REAL(8) :: TM, AA
TM=0
DO K=1,M
AA=(GM(K)-DZ(K))**2
TM=TM+AA
IF(DZ(K).EQ.0) THEN
B(K)=GM(K)
ELSE
B(K)=GM(K)/DZ(K)-1.0
END IF
END DO
TM=SQRT(TM/M)
RETURN
END SUBROUTINE BI
! ================================================================ ! -------------------用奇异值分解求解线性方程组-------------------
! ================================================================ ! N----整型变量,输入参数,A的列数。
! M----整型变量,输入参数,A的行数。
! A----实型数组,输入、输出参数,输入时按列存放实矩阵A,输出矩阵U。
! W----实型数组,输入参数,输出奇异值。
! V----实型数组,输入参数,输出变换矩阵V。
SUBROUTINE SVDCMP(A,M,N,W,V)
IMPLICIT NONE
INTEGER :: M, N, NMAX
REAL(8) :: A(M,N),W(N),V(N,N)
PARAMETER(NMAX=6000)
INTEGER :: I,L,K,NM,J,ITS
REAL(8) :: G,S,SCALE,F,H,C,X,Y,Z,RV1(NMAX),ANORM,QQ
REAL(8),EXTERNAL:: JUDSQ
IF(M G=0.0 SCALE=0.0 ANORM=0.0 DO I=1,N L=I+1 RV1(I)=SCALE*G G=0.0 S=0.0 SCALE=0.0 IF(I<=M) THEN DO K=I,M SCALE=SCALE+ABS(A(K,I)) END DO IF(SCALE/=0.0) THEN DO K=I,M A(K,I)=A(K,I)/SCALE S=S+A(K,I)*A(K,I) END DO F=A(I,I) G=-SIGN(SQRT(ABS(S)),F) H=F*G-S DO J=L,N S=0.0 DO K=I,M S=S+A(K,I)*A(K,J) END DO F=S/H DO K=I,M A(K,J)=A(K,J)+F*A(K,I) END DO END DO DO K=I,M A(K,I)=SCALE*A(K,I) END DO ENDIF ENDIF W(I)=SCALE*G G=0.0 S=0.0 SCALE=0.0 IF ((I<=M).AND.(I/=N)) THEN DO K=L,N SCALE=SCALE+ABS(A(I,K)) END DO IF (SCALE/=0.0) THEN DO K=L,N A(I,K)=A(I,K)/SCALE S=S+A(I,K)*A(I,K) END DO F=A(I,L) G=-SIGN(SQRT(ABS(S)),F) H=F*G-S A(I,L)=F-G DO K=L,N RV1(K)=A(I,K)/H END DO DO J=L,M S=0.0 DO K=L,N S=S+A(J,K)*A(I,K) END DO DO K=L,N A(J,K)=A(J,K)+S*RV1(K) END DO END DO DO K=L,N A(I,K)=SCALE*A(I,K) END DO ENDIF ENDIF ANORM=MAX(ANORM,(ABS(W(I))+ABS(RV1(I)))) END DO DO I=N,1,-1 IF (I IF (G/=0.0) THEN DO J=L,N V(J,I)=(A(I,J)/A(I,L))/G END DO S=0.0 DO K=L,N S=S+A(I,K)*V(K,J) END DO DO K=L,N V(K,J)=V(K,J)+S*V(K,I) END DO END DO ENDIF DO J=L,N V(I,J)=0.0 V(J,I)=0.0 END DO ENDIF V(I,I)=1.0 G=RV1(I) L=I END DO DO I=N,1,-1 L=I+1 G=W(I) DO J=L,N A(I,J)=0.0 END DO IF (G/=0.0) THEN G=1.0/G DO J=L,N S=0.0 DO K=L,M S=S+A(K,I)*A(K,J) END DO F=(S/A(I,I))*G DO K=I,M A(K,J)=A(K,J)+F*A(K,I) END DO END DO DO J=I,M A(J,I)=A(J,I)*G END DO ELSE DO J=I,M A(J,I)=0.0 END DO ENDIF A(I,I)=A(I,I)+1.0 END DO DO K=N,1,-1 DO ITS=1,30 DO L=K,1,-1 NM=L-1 IF((ABS(RV1(L))+ANORM)==ANORM) GOTO 10 IF((ABS(W(NM))+ANORM)==ANORM) GOTO 20 END DO 20 C=0.0 S=1.0 DO I=L,K F=S*RV1(I) RV1(I)=C*RV1(I) IF((ABS(F)+ANORM).EQ.ANORM) GOTO 10 G=W(I) H=JUDSQ(F,G) W(I)=H H=1.0/H C=(G*H) S=-(F*H) DO J=1,M Y=A(J,NM) Z=A(J,I) A(J,NM)=(Y*C)+(Z*S) A(J,I)=-(Y*S)+(Z*C) END DO END DO 10 Z=W(K) IF(L==K) THEN IF(Z<0.0) THEN W(K)=-Z DO J=1,N V(J,K)=-V(J,K) END DO ENDIF EXIT ENDIF IF(ITS==10000) PAUSE 'NO CONVERGENCE IN 400 ITERATIONS' X=W(L) NM=K-1 Y=W(NM) G=RV1(NM) H=RV1(K) F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y) QQ=1. G=JUDSQ(F,QQ) F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X C=1.0 S=1.0 DO J=L,NM I=J+1 G=RV1(I) Y=W(I) H=S*G G=G*C Z=JUDSQ(F,H) RV1(J)=Z C=F/Z S=H/Z F=(X*C)+(G*S) G=-(X*S)+(G*C) H=Y*S Y=Y*C DO NM=1,N X=V(NM,J) Z=V(NM,I) V(NM,J)=(X*C)+(Z*S) V(NM,I)=-(X*S)+(Z*C) END DO Z=JUDSQ(F,H) W(J)=Z IF(Z/=0.0) THEN Z=1.0/Z C=F*Z S=H*Z ENDIF F=(C*G)+(S*Y) X=-(S*G)+(C*Y) DO NM=1,M Y=A(NM,J) Z=A(NM,I) A(NM,J)=(Y*C)+(Z*S) A(NM,I)=-(Y*S)+(Z*C) END DO END DO RV1(L)=0.0 RV1(K)=F W(K)=X END DO END DO RETURN END SUBROUTINE SVDCMP ! ================================================================ ! U,W,V输入参数,由子程序SVDCMP所输出。 ! N----整型变量,输入参数,A的列数。 ! M----整型变量,输入参数,A的行数。 ! B----实型数组,输入参数,输入方程组右端的向量。 ! X----实型数组,输入参数,输出方程组的最小二乘解。 SUBROUTINE SVBKSB(A,W,V,M,N,B,X) IMPLICIT NONE INTEGER :: M, N, NMAX REAL(8) :: A(M,N),W(N),V(N,N),B(M),X(N) PARAMETER (NMAX=10000) REAL(8) :: TNMAX(NMAX),S INTEGER :: I, J, JJ DO J=1,N S=0.0 IF (W(J)/=0.0) THEN DO I=1,M S=S+A(I,J)*B(I) END DO S=S/W(J) END IF TNMAX(J)=S END DO DO J=1,N S=0.0 DO JJ=1,N S=S+V(J,JJ)*TNMAX(JJ) END DO X(J)=S END DO RETURN END SUBROUTINE SVBKSB ! ================================================================ ! 判别根号里面小于等于零的情况 FUNCTION JUDSQ(A,B) IMPLICIT NONE REAL(8) :: A,B,JUDSQ REAL(8) :: ABSA,ABSB ABSA=ABS(A) ABSB=ABS(B) IF(ABSA.GT.ABSB)THEN JUDSQ=ABSA*SQRT(1.+(ABSB/ABSA)**2) ElSE IF(ABSB.EQ.0.)THEN JUDSQ=0. ElSE JUDSQ=ABSB*SQRT(1.+(ABSA/ABSB)**2) ENDIF ENDIF RETURN END FUNCTION JUDSQ ! ================================================================ ! 优化奇异值分解法的解,得出较为精确的广义逆解。 ! 见《地球物理数据处理基础》108页 SUBROUTINE WIGGINSMETHOD(N,W) IMPLICIT NONE INTEGER :: K, N REAL(8) :: WMAX, W(N),WMIN WMAX=0.0 DO K=1,N IF(W(K)>WMAX) WMAX=W(K) !找到最大的奇异值。 END DO WMIN=WMAX*(1.0E-03) !定义最小值。 DO K=1,N IF(W(K) END DO RETURN END SUBROUTINE WIGGINSMETHOD ! ================================================================ ! 调用三次样条插值 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