zhangyiss/DSurfTomo
1
0
Fork 0
mirror of https://github.com/HongjianFang/DSurfTomo.git synced 2025-08-03 01:34:13 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

add model uncertainty estimation

Browse Source 
add model uncertainty estimation with random models, also indent all the
code to make them a little bit easy to read
bug fix about shift when output the velocity model
This commit is contained in:
Hongjian Fang 2016-08-07 20:05:49 +02:00
parent 851eb3418f
commit 4f2f7cbce5
7 changed files with 3316 additions and 3207 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
example/DSurfTomo.in
View File

@ -5,7 +5,7 @@ surfdataTB.dat c: data file
18 18 9 c: nx ny nz (grid number in lat lon and depth direction) 18 18 9 c: nx ny nz (grid number in lat lon and depth direction)
25.2 121.35 c: goxd gozd (upper left point,[lat,lon]) 25.2 121.35 c: goxd gozd (upper left point,[lat,lon])
0.015 0.017 c: dvxd dvzd (grid interval in lat and lon direction) 0.015 0.017 c: dvxd dvzd (grid interval in lat and lon direction)
449 c: nsrc*maxf 25i c: nsrc*maxf
4.0 0.0 c: weight damp 4.0 0.0 c: weight damp
3 c: nsublayer (numbers of sublayers for each grid interval:grid --> layer) 3 c: nsublayer (numbers of sublayers for each grid interval:grid --> layer)
0.5 2.8 c: minimum velocity, maximum velocity 0.5 2.8 c: minimum velocity, maximum velocity
@ -19,3 +19,5 @@ surfdataTB.dat c: data file
0 c: synthetic flag(0:real data,1:synthetic) 0 c: synthetic flag(0:real data,1:synthetic)
0.02 c: noiselevel 0.02 c: noiselevel
2.5 c: threshold 2.5 c: threshold
1 c: modest (1: estimate model variation, 0: no estimation)
30 c: number of random models

2422
src/CalSurfG.f90
View File

File diff suppressed because it is too large Load Diff

BIN
src/DSurfTomo
View File

Binary file not shown.

16
src/aprod.f90
View File

@ -4,11 +4,11 @@
!c Version for use with sparse matrix specified by !c Version for use with sparse matrix specified by
!c output of subroutine sparse for use with LSQR !c output of subroutine sparse for use with LSQR
subroutine aprod(mode, m, n, x, y, leniw, lenrw, iw, rw) subroutine aprod(mode, m, n, x, y, leniw, lenrw, iw, rw)
implicit none implicit none
!c Parameters: !c Parameters:
integer mode ! ==1: Compute y = y + a*x integer mode ! ==1: Compute y = y + a*x
! y is altered without changing x ! y is altered without changing x
! ==2: Compute x = x + a(transpose)*y ! ==2: Compute x = x + a(transpose)*y
@ -23,38 +23,38 @@
! iw[iw[1]+2:2*iw[1]+1] Column indices ! iw[iw[1]+2:2*iw[1]+1] Column indices
real rw(lenrw) ! [1..iw[1]] Non-zero elements of a real rw(lenrw) ! [1..iw[1]] Non-zero elements of a
!c Local variables: !c Local variables:
integer i1 integer i1
integer j1 integer j1
integer k integer k
integer kk integer kk
!c set the ranges the indices in vector iw !c set the ranges the indices in vector iw
kk=iw(1) kk=iw(1)
i1=1 i1=1
j1=kk+1 j1=kk+1
!c main iteration loop !c main iteration loop
do k = 1,kk do k = 1,kk
if (mode.eq.1) then if (mode.eq.1) then
!c compute y = y + a*x !c compute y = y + a*x
y(iw(i1+k)) = y(iw(i1+k)) + rw(k)*x(iw(j1+k)) y(iw(i1+k)) = y(iw(i1+k)) + rw(k)*x(iw(j1+k))
else else
!c compute x = x + a(transpose)*y !c compute x = x + a(transpose)*y
x(iw(j1+k)) = x(iw(j1+k)) + rw(k)*y(iw(i1+k)) x(iw(j1+k)) = x(iw(j1+k)) + rw(k)*y(iw(i1+k))
endif endif
enddo enddo
! 100 continue ! 100 continue
return return
end end

46
src/delsph.f90
View File

@ -1,28 +1,28 @@
subroutine delsph(flat1,flon1,flat2,flon2,del) subroutine delsph(flat1,flon1,flat2,flon2,del)
implicit none implicit none
real,parameter:: R=6371.0 real,parameter:: R=6371.0
REAL,parameter:: pi=3.1415926535898 REAL,parameter:: pi=3.1415926535898
real flat1,flat2 real flat1,flat2
real flon1,flon2 real flon1,flon2
real del real del
real dlat real dlat
real dlon real dlon
real lat1 real lat1
real lat2 real lat2
real a real a
real c real c
!dlat=(flat2-flat1)*pi/180 !dlat=(flat2-flat1)*pi/180
!dlon=(flon2-flon1)*pi/180 !dlon=(flon2-flon1)*pi/180
!lat1=flat1*pi/180 !lat1=flat1*pi/180
!lat2=flat2*pi/180 !lat2=flat2*pi/180
dlat=flat2-flat1 dlat=flat2-flat1
dlon=flon2-flon1 dlon=flon2-flon1
lat1=pi/2-flat1 lat1=pi/2-flat1
lat2=pi/2-flat2 lat2=pi/2-flat2
a=sin(dlat/2)*sin(dlat/2)+sin(dlon/2)*sin(dlon/2)*cos(lat1)*cos(lat2) a=sin(dlat/2)*sin(dlat/2)+sin(dlon/2)*sin(dlon/2)*cos(lat1)*cos(lat2)
c=2*atan2(sqrt(a),sqrt(1-a)) c=2*atan2(sqrt(a),sqrt(1-a))
del=R*c del=R*c
end subroutine end subroutine

6
src/gaussian.f90
View File

@ -1,6 +1,6 @@
real function gaussian() real function gaussian()
implicit none implicit none
! real rd ! real rd
real x1,x2,w,y1 real x1,x2,w,y1
real y2 real y2
@ -28,4 +28,4 @@
use_last=1 use_last=1
endif endif
gaussian=y1 gaussian=y1
end function end function

185
src/main.f90
View File

@ -1,28 +1,28 @@
! CODE FOR SURFACE WAVE TOMOGRAPHY USING DISPERSION MEASUREMENTS ! CODE FOR SURFACE WAVE TOMOGRAPHY USING DISPERSION MEASUREMENTS
! VERSION: ! VERSION:
! 1.0 ! 1.0
! AUTHOR: ! AUTHOR:
! HONGJIAN FANG. fanghj@mail.ustc.edu.cn ! HONGJIAN FANG. fanghj@mail.ustc.edu.cn
! PURPOSE: ! PURPOSE:
! DIRECTLY INVERT SURFACE WAVE DISPERSION MEASUREMENTS FOR 3-D ! DIRECTLY INVERT SURFACE WAVE DISPERSION MEASUREMENTS FOR 3-D
! STUCTURE WITHOUT THE INTERMEDIATE STEP OF CONSTUCTION THE PHASE ! STUCTURE WITHOUT THE INTERMEDIATE STEP OF CONSTUCTION THE PHASE
! OR GROUP VELOCITY MAPS. ! OR GROUP VELOCITY MAPS.
! REFERENCE: ! REFERENCE:
! Fang, H., Yao, H., Zhang, H., Huang, Y. C., & van der Hilst, R. D. ! Fang, H., Yao, H., Zhang, H., Huang, Y. C., & van der Hilst, R. D.
! (2015). Direct inversion of surface wave dispersion for ! (2015). Direct inversion of surface wave dispersion for
! three-dimensional shallow crustal structure based on ray tracing: ! three-dimensional shallow crustal structure based on ray tracing:
! methodology and application. Geophysical Journal International, ! methodology and application. Geophysical Journal International,
! 201(3), 1251-1263. ! 201(3), 1251-1263.
! HISTORY: ! HISTORY:
! 2015/01/31 START TO REORGONIZE THE MESSY CODE ! 2015/01/31 START TO REORGONIZE THE MESSY CODE
! !
program SurfTomo program SurfTomo
use lsmrModule, only:lsmr use lsmrModule, only:lsmr
use lsmrblasInterface, only : dnrm2 use lsmrblasInterface, only : dnrm2
implicit none implicit none
! VARIABLE DEFINE ! VARIABLE DEFINE
character inputfile*80 character inputfile*80
character logfile*100 character logfile*100
@ -101,24 +101,37 @@
real maxnorm real maxnorm
real threshold,threshold0 real threshold,threshold0
! OPEN FILES FIRST TO OUTPUT THE PROCESS ! FOR MODEL VARIATION
!------------------------------------------------
integer idx
integer counte
real stdvs
integer numrand
real,allocatable,dimension(:,:)::modstat
real,allocatable,dimension(:)::modsig
real gaussian
external gaussian
integer modest
counte=0
! OPEN FILES FIRST TO OUTPUT THE PROCESS
open(34,file='IterVel.out') open(34,file='IterVel.out')
nout=36 nout=36
open(nout,file='lsmr.txt') open(nout,file='lsmr.txt')
! OUTPUT PROGRAM INFOMATION ! OUTPUT PROGRAM INFOMATION
write(*,*) write(*,*)
write(*,*),' S U R F T O M O' write(*,*),' S U R F T O M O'
write(*,*),'PLEASE contact Hongjain Fang & write(*,*),'PLEASE contact Hongjain Fang &
(fanghj@mail.ustc.edu.cn) if you find any bug' (fanghj@mail.ustc.edu.cn) if you find any bug'
write(*,*) write(*,*)
! READ INPUT FILE ! READ INPUT FILE
if (iargc() < 1) then if (iargc() < 1) then
write(*,*) 'input file [SurfTomo.in(default)]:' write(*,*) 'input file [DSurfTomo.in(default)]:'
read(*,'(a)') inputfile read(*,'(a)') inputfile
if (len_trim(inputfile) <=1 ) then if (len_trim(inputfile) <=1 ) then
inputfile = 'SurfTomo.in' inputfile = 'DSurfTomo.in'
else else
inputfile = inputfile(1:len_trim(inputfile)) inputfile = inputfile(1:len_trim(inputfile))
endif endif
@ -205,11 +218,13 @@
read(10,*)ifsyn read(10,*)ifsyn
read(10,*)noiselevel read(10,*)noiselevel
read(10,*) threshold0 read(10,*) threshold0
read(10,*) modest
read(10,*) numrand
close(10) close(10)
nrc=nsrc nrc=nsrc
kmax=kmaxRc+kmaxRg+kmaxLc+kmaxLg kmax=kmaxRc+kmaxRg+kmaxLc+kmaxLg
! READ MEASUREMENTS ! READ MEASUREMENTS
open(unit=87,file=datafile,status='old') open(unit=87,file=datafile,status='old')
allocate(scxf(nsrc,kmax),sczf(nsrc,kmax),& allocate(scxf(nsrc,kmax),sczf(nsrc,kmax),&
rcxf(nrc,nsrc,kmax),rczf(nrc,nsrc,kmax),stat=checkstat) rcxf(nrc,nsrc,kmax),rczf(nrc,nsrc,kmax),stat=checkstat)
@ -286,6 +301,10 @@
allocate(norm(maxvp), stat=checkstat) allocate(norm(maxvp), stat=checkstat)
allocate(vsf(nx,ny,nz), stat=checkstat) allocate(vsf(nx,ny,nz), stat=checkstat)
allocate(vsftrue(nx,ny,nz), stat=checkstat) allocate(vsftrue(nx,ny,nz), stat=checkstat)
! FOR MODEL VARIATION
!------------------------------------------------
allocate(modstat(numrand,maxvp))
allocate(modsig(maxvp))
allocate(rw(maxnar), stat=checkstat) allocate(rw(maxnar), stat=checkstat)
if(checkstat > 0)then if(checkstat > 0)then
@ -302,7 +321,7 @@
allocate(cbst(dall+maxvp),dsyn(dall),datweight(dall),wt(dall+maxvp),dtres(dall+maxvp),& allocate(cbst(dall+maxvp),dsyn(dall),datweight(dall),wt(dall+maxvp),dtres(dall+maxvp),&
stat=checkstat) stat=checkstat)
! MEASUREMENTS STATISTICS AND READ INITIAL MODEL ! MEASUREMENTS STATISTICS AND READ INITIAL MODEL
write(*,'(a,i7)') 'Number of all measurements',dall write(*,'(a,i7)') 'Number of all measurements',dall
open(10,file='MOD',status='old') open(10,file='MOD',status='old')
@ -318,7 +337,7 @@
! CHECKERBOARD TEST ! CHECKERBOARD TEST
if (ifsyn == 1) then if (ifsyn == 1) then
write(*,*) 'Checkerboard Resolution Test Begin' write(*,*) 'Checkerboard Resolution Test Begin'
vsftrue = vsf vsftrue = vsf
@ -340,14 +359,14 @@
! ITERATE UNTILL CONVERGE ! ITERATE UNTILL CONVERGE
writepath = 0 writepath = 0
do iter = 1,maxiter do iter = 1,maxiter
iw = 0 iw = 0
rw = 0.0 rw = 0.0
col = 0 col = 0
! COMPUTE SENSITIVITY MATRIX ! COMPUTE SENSITIVITY MATRIX
if (iter == maxiter) then if (iter == maxiter) then
writepath = 1 writepath = 1
open(40,file='raypath.out') open(40,file='raypath.out')
@ -390,7 +409,7 @@
write(66,*)'Maximum and Average DWS values:',maxnorm,averdws write(66,*)'Maximum and Average DWS values:',maxnorm,averdws
write(66,*)'Threshold is:',threshold write(66,*)'Threshold is:',threshold
! WRITE OUT RESIDUAL FOR THE FIRST AND LAST ITERATION ! WRITE OUT RESIDUAL FOR THE FIRST AND LAST ITERATION
if(iter.eq.1) then if(iter.eq.1) then
open(88,file='residualFirst.dat') open(88,file='residualFirst.dat')
do i=1,dall do i=1,dall
@ -409,7 +428,7 @@
endif endif
! ADDING REGULARIZATION TERM ! ADDING REGULARIZATION TERM
weight=dnrm2(dall,cbst,1)**2/dall*weight0 weight=dnrm2(dall,cbst,1)**2/dall*weight0
nar_tmp=nar nar_tmp=nar
nars=0 nars=0
@ -465,7 +484,7 @@
enddo enddo
if (nar > maxnar) stop 'increase sparsity fraction(spfra)' if (nar > maxnar) stop 'increase sparsity fraction(spfra)'
! CALLING IRLS TO SOLVE THE PROBLEM ! CALLING IRLS TO SOLVE THE PROBLEM
leniw = 2*nar+1 leniw = 2*nar+1
lenrw = nar lenrw = nar
@ -487,6 +506,10 @@
if(istop==3) print*,'istop = 3, large condition number' if(istop==3) print*,'istop = 3, large condition number'
do i =1,dall
cbst(i)=cbst(i)/datweight(i)
enddo
mean = sum(cbst(1:dall))/dall mean = sum(cbst(1:dall))/dall
std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2) std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
write(*,'(i2,a)'),iter,'th iteration...' write(*,'(i2,a)'),iter,'th iteration...'
@ -520,7 +543,7 @@
enddo enddo
enddo enddo
enddo enddo
write(34,*)',OUTPUT S VELOCITY AT ITERATION',iter write(34,*)'OUTPUT S VELOCITY AT ITERATION',iter
do k=1,nz do k=1,nz
do j=1,ny do j=1,ny
write(34,'(100f7.3)') (vsf(i,j,k),i=1,nx) write(34,'(100f7.3)') (vsf(i,j,k),i=1,nx)
@ -535,7 +558,7 @@
enddo !end iteration enddo !end iteration
! OUTPUT THE VELOCITY MODEL ! OUTPUT THE VELOCITY MODEL
write(*,*),'Program finishes successfully' write(*,*),'Program finishes successfully'
write(66,*),'Program finishes successfully' write(66,*),'Program finishes successfully'
@ -547,8 +570,8 @@
do k=1,nz-1 do k=1,nz-1
do j=1,ny-2 do j=1,ny-2
do i=1,nx-2 do i=1,nx-2
write(65,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsftrue(i,j,k) write(65,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsftrue(i+1,j+1,k)
write(63,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i,j,k) write(63,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i+1,j+1,k)
enddo enddo
enddo enddo
enddo enddo
@ -568,7 +591,7 @@
do k=1,nz-1 do k=1,nz-1
do j=1,ny-2 do j=1,ny-2
do i=1,nx-2 do i=1,nx-2
write(64,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i,j,k) write(64,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i+1,j+1,k)
enddo enddo
enddo enddo
enddo enddo
@ -583,6 +606,73 @@
close(40) close(40)
close(nout) !close lsmr.txt close(nout) !close lsmr.txt
close(66) !close surf_tomo.log close(66) !close surf_tomo.log
! USE RANDOM MODEL TO OBTAIN THE MODEL VARIATION
!modest = 1
if (modest ==1) then
write(*,*) 'model variation estimation begin...'
do iter = 1,numrand
call init_random_seed()
vsftrue=vsf
DO K=1,NZ-1
DO J=2,NY-1
DO I=2,NX-1
idx = (k-1)*(ny-2)*(nx-2)+(j-2)*(nx-2)+i-1
dv(idx) = 0.1/EXP(2*NORM(idx)/maxnorm)*gaussian()
VSFTRUE(I,J,K) = VSF(I,J,K)+dv(idx)
ENDDO
ENDDO
ENDDO
write(*,*),'maximum and minimum velocity variation',maxval(dv),minval(dv)
call synthetic(nx,ny,nz,maxvp,vsftrue,dsyn,&
goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
scxf,sczf,rcxf,rczf,nrc1,nsrc1,kmax,&
nsrc,nrc,0.0)
do i = 1,dall
cbst(i) = obst(i) - dsyn(i)
enddo
write(*,*), dnrm2(dall,cbst,1)/sqrt(real(dall)), 1.05*std_devs
if (dnrm2(dall,cbst,1)/sqrt(real(dall)) < 1.05*std_devs) then
counte = counte + 1
modstat(counte,:) = dv
endif
enddo ! iteration for random models
write(*,*),'number of of models satisfy requirements',counte
modsig = 1.0
if (counte>0) then
do i=1,maxvp
!statis
!mean = sum(cbst(1:dall))/dall
!std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
mean = sum(modstat(1:counte,i))/counte
stdvs = sqrt(sum(modstat(1:counte,i)**2)/counte-mean**2)
modsig(i) = stdvs
enddo
endif
write(*,*),'write model variation to "model_variation.dat"'
open (64,file='model_variation.dat')
do k=1,nz-1
do j=1,ny-2
do i=1,nx-2
idx = (k-1)*(ny-2)*(nx-2)+(j-1)*(nx-2)+i
write(64,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),modsig(idx)
enddo
enddo
enddo
close(64)
write(*,*) 'finishing model variation estimation'
endif
deallocate(obst) deallocate(obst)
deallocate(dsyn) deallocate(dsyn)
deallocate(dist) deallocate(dist)
@ -611,4 +701,21 @@
deallocate(tLg) deallocate(tLg)
endif endif
end program end program
!-----------------------------------------------------------------------
! Generate seed for random number generator of fortran
! Note: only need to be called once inside one program
!-----------------------------------------------------------------------
subroutine init_random_seed()
integer :: i,n,clock
integer,dimension(:),allocatable :: seed
call random_seed(size=n)
allocate(seed(n))
call system_clock(count=clock)
seed=clock+37*(/(i-1,i=1,n)/)
call random_seed(PUT=seed)
deallocate(seed)
end subroutine