radar2d4.f90

基于matlab的探地雷达信号处理开源程序

program radar2d4
!
!********************************************************************
!*                SPLIT-STEP 2D GPR MODELLING                       *	
!*         After Grandjean, G.and Durand, H., 1999,                 *
!*         Computers and Geosciences25 141-149.                     *
!*         Based on the C program rada2d4 by G. Grandjean           *
!*         Coded by : Andreas Tzanis, Department of Geophysics,     *
!*                    University of Athens, [e-mail]     *
!********************************************************************
!
implicit none
real(8),parameter   :: pi  = 3.141592653589793d0
real(8),parameter   :: mu0 =4.0d-7*pi
real(8),parameter   :: eps0=8.8592d-12
integer(4),parameter :: mhz = 1000000
! - - - local declarations - - -
integer(2) :: status
integer    :: npfft, nargs, alloc_err, istat, m2
integer    :: argc,fd,fdo,fdq,fdm
integer    :: n,ikx,ix,it,iz, iw, iw1, iw2, iwc, i, j, k
integer    :: nx2,nt2,nt,nz,ntr,f1,f2,fc,fc2 
integer    :: tet1,tet2,percent,dix
real(4)    :: kz,ctr,radc,difc,tet,dx,dz,nq,wr,expq,source 
real(4)    :: band,cref,sigma,kx,den,rshift,dt,den1,den2,bomega2
real(4)    :: aomega,bomega,aomegaxz,bomegaxz,vomegaxz
complex(4) :: shift,phase,denkc,denc,denkz,zdif,komegaxz,shift1,shift2,cref1,cref2
real(4),allocatable :: per(:,:),mu(:,:),v(:,:),q(:,:),vmean(:)
real(4),allocatable :: qm(:),omega(:),waven(:), g(:)
real(4),allocatable :: rcw2(:), icw2(:)
complex*16,allocatable :: cp(:,:),rc(:,:),komega(:)
character (len=256)  :: permit,muname,qname,outfname
character (len=256) :: argv
! - - - i/o channels - - -!
data fd,fdo,fdq,fdm/20, 21, 22, 23/

! - - - begin - - -
argc = NARGS()       ! #args from command line
!**** COMMAND LINE ARGUMENTS
if (argc == 13) then
   call getarg(int2(1),argv,status)
   permit  = adjustl(argv(1:status))
   call getarg(int2(2),argv,status)
   muname  = adjustl(argv(1:status))
   call getarg(int2(3),argv,status)
   qname = adjustl(argv(1:status))
   call getarg(int2(4),argv,status)
   outfname  = adjustl(argv(1:status))
   call getarg(int2(5),argv,status)
   read(argv(1:status),*) ntr
   call getarg(int2(6),argv,status)
   read (argv(1:status),*) nt
   call getarg(int2(7),argv,status)
   read (argv(1:status),*) dx
   call getarg(int2(8),argv,status)
   read (argv(1:status),*) nz
   call getarg(int2(9),argv,status)
   read (argv(1:status),*) dz
   call getarg(int2(10),argv,status)
   read (argv(1:status),*) dt
   call getarg(int2(11),argv,status)
   read (argv(1:status),*) fc
   call getarg(int2(12),argv,status)
   read (argv,*) band
else
  write (*,'(a)') " relative permittivity file name ? "
  read (*,'(a)') permit
  write (*,'(a)') " magnetic permeability file name ? "
  read (*,'(a)') muname
  write (*,'(a)') " quality factor file name ? "
  read (*,'(a)') qname
  write (*,'(a)') " output file ? "
  read (*,'(a)') outfname
  write (*,'(a)') " number of traces ?"
  read (*,*) ntr
  write (*,'(a)') " number of time samples ?"
  read (*,*) nt
  write (*,'(a)') " trace spacing (m)  ?"
  read (*,*) dx
  write (*,'(a)') " number of samples in z  ?"
  read (*,*) nz
  write (*,'(a)') " sampling interval along z (m)  ?"
  read (*,*) dz
  write (*,'(a)') " sampling interval in time (ns)  ?"
  read (*,*) dt
  write (*,'(a)') " antenna frequency (mhz) ?"
  read (*,*) fc
  write (*,'(a)') " number of bandwidths ?"
  read (*,*) band
end if

! Determine sizes of FFTs and work arrays
nx2  = npfft(int(ntr*1.1),2*ntr)
nt2  = npfft(int(nt*1.1),2*nt)
! determine frequency range 
tet1 = 60
tet2 = 80
wr   = 2.*pi*fc*mhz
dt   = dt/1.0e9
f1   = int(fc*mhz*0.25)
f2   = int(fc*mhz*2)
fc2  = int(fc*mhz)
iw1  = int(nt2*dt*f1)
iw2  = int(nt2*dt*f2)
iwc  = int(nt2*dt*fc2)
if (band < 0.0 .or. band > 0.0 ) then
   iw1 = 1
   iw2 = nt2/2+1
end if
write (*,'(a)') " > Modeling: parameters input OK."

!Allocate all complex arrays
 allocate(cp(nt2,nx2),rc(nz,nx2),komega(nt2),stat=alloc_err)
!Allocate real 2-D arrays  
 allocate(per(nz,ntr),mu(nz,ntr),v(nz,ntr),q(nz,ntr),stat=alloc_err)
!Allocate real 1-D arrays
 allocate(vmean(nz),waven(nx2),omega(nt2),stat=alloc_err)
 m2 = max(nx2,nt2)
 allocate(rcw2(m2),icw2(m2),stat=alloc_err)
 allocate(qm(nz),g(nt2), stat=alloc_err)
  WRITE (*,'(a)') "> Modeling: memory allocation OK."

! - - - import data - - -
open(fd,file=permit(1:len_trim(permit)),action='read',form='binary')
open(fdq,file=qname(1:len_trim(qname)),action='read',form='binary')
open(fdm,file=muname(1:len_trim(muname)),action='read',form='binary')
open(fdo,file=outfname(1:len_trim(outfname)),action='write')
do ix = 1,ntr
   read(fd) (vmean(i),i=1,nz)
   do iz = 1,nz
     per(iz,ix) = vmean(iz)
   end do
   read(fdm) (vmean(i),i=1,nz)
   do iz = 1,nz
     mu(iz,ix) = vmean(iz)
   end do
   read(fdq) (vmean(i),i=1,nz)
   do iz = 1,nz
     q(iz,ix) = vmean(iz)
   end do
end do
close (fd)
close (fdm)
close (fdq)
!**** compute velocity
do ix = 1,ntr
   do iz = 1,nz
     v(iz,ix) = 1.0/(sqrt(mu0*mu(iz,ix)*eps0*per(iz,ix))*  &
         cos((pi/4)*(1.0-((2/pi)*atan(q(iz,ix))))))
   end do
end do
write (*,'(a)') "> modeling: model input ok."
!**** compute mean velocity and quality factor
call mean_velocity(vmean, v, ntr, nz)
call mean_q(qm, q, ntr, nz)
!**** initialize work arrays
do ix = 1,nx2
  do iz = 1,nz
    rc(iz,ix) = cmplx(0.0,0.0)
  end do
  do iz = 1,nt2
    cp(iz,ix) = cmplx(0.0,0.0)
  end do
end do
!**** reflection coefficient
DO ix = 1,ntr
   DO iz = 2,nz
      cref1 = csqrt( cmplx(mu(iz,ix), 0.0) / &
              cmplx( per(iz,ix)*sin(atan(Q(iz,ix))), &
	          cos(per(iz,ix)*atan(Q(iz,ix))) ))
      cref2 = csqrt( cmplx(mu(iz-1,ix),0.0) / &
              cmplx( per(iz-1,ix)*sin(atan(Q(iz-1,ix))), &
              cos(per(iz-1,ix)*atan(Q(iz-1,ix))) ))
      rc(iz,ix) = (cref1 - cref2)/(cref1 + cref2)
   ENDDO
enddo
!**** Transform X -> kx
do iz = 1,nz
   do ix = 1,ntr
	  rcw2(ix)   = real(rc(iz,ix))
	  icw2(ix)   = imag(rc(iz,ix))
   end do
   do ix = ntr+1,nx2
	  rcw2(ix)   = 0.0
	  icw2(ix)   = 0.0
   end do
   call fft(rcw2,icw2,nx2,nx2,nx2,-1)
   do ix = 1,nx2
      rc(iz,ix) = cmplx(rcw2(ix), icw2(ix))
   end do
end do
write (*,'(a)') "> Modeling: FFT along X OK."
!**** compute angular frequenciws and wavenumbers
do ikx = 1,nx2
   kx = ikx*2.0*pi/nx2
   if (kx > pi) kx = 2.0*pi-kx
   kx = kx/dx
   waven(ikx) = kx
end do
do iw = iw1,iw2
   omega(iw) = (iw*2.0*pi/nt2)/dt
end do
if (band > 0.0) then
   sigma = (iw2-iw1)*band
   do iw = iw1,iw2
       g(iw) = exp(-1.0*((iw-iwc)/sigma)*((iw-iwc)/sigma))
   end do
end if
write (*,'(a)') "> Modeling: processing wavenumbers OK."
!*******************************
!*     LOOP OVER DEPTHS
!*******************************
dix = 0
percent = 0
write(*,'(''> modeling z='',f4.2,''m, ('',i3,'' percent)'')')  nz*dz,percent
do iz = nz,1,-1
   percent = int(((10*(nz-iz)*dz)/(nz*dz)))
   percent = percent*10
   if (percent > dix) then
      write (*,'(''> Modeling Z='',f4.2,''m, ('',i3,'' percent)'')')  iz*dz,percent
      dix = dix+ 10
   end if
   nq = (2.0/pi)*atan(qm(iz))
   expq = (1. - nq)/2.0
!*******************************
!*    LOOP OVER WAVENUMBERS
!*******************************
   do ikx = 1,nx2
      den = (waven(ikx)*waven(ikx))
      denkc = cmplx(den,0.0)
!*******************************
!*    LOOP OVER FREQUENCIES
!*******************************
      do iw = iw1,iw2
         ctr = (omega(iw)/wr)**expq
         bomega = omega(iw)/( vmean(iz) * ctr )
         aomega = bomega*dtan((1.0d0 - nq)*(pi/4.0d0))
         komega(iw) = dcmplx(bomega,aomega)
!*******************************
!*  INCLUDE RADIATION PATTERN
!*******************************
         bomega2 = (omega(iw)/vmean(1))*(omega(iw)/vmean(1))
         difc = bomega2-den
         if (difc > 0.0) then
            kz = sqrt(difc)
            tet = atan(waven(ikx)/kz)
            radc = source(tet1,tet2,tet)
          else
            radc = 1.0
          end if
          denc = komega(iw)*komega(iw)
          zdif = denc - denkc
          denkz = csqrt(zdif)
          den1 = real(denkz)
          den2 = imag(denkz)
          if ( den2 > 0.0) then
             rshift = exp((-1.0)*dz*den2)
          else
             rshift = exp((1.0)*dz*den2)
          end if
          shift1 = cexp(cmplx(0.0, dz*den1))
          shift2 = shift1*rshift
          cp(iw,ikx) = cp(iw,ikx)*shift2
          if (band > 0.0) cp(iw,ikx) = cp(iw,ikx)*g(iw)
          cp(iw,ikx) = cp(iw,ikx) + rc(iz,ikx)
          cp(iw,ikx) = cp(iw,ikx)*radc
      end do
   end do
!**** Inverse transform Kx -> X
   do iw = iw1,iw2
      do ikx = 1,nx2
  	     rcw2(ikx)   = real(cp(iw,ikx)*(1.0/nx2))
	     icw2(ikx)   = imag(cp(iw,ikx)*(1.0/nx2))
      end do
      call fft(rcw2,icw2,nx2,nx2,nx2,1)
      do ikx = 1,nx2
         cp(iw,ikx) = cmplx(rcw2(ikx), icw2(ikx))
     end do
   end do
!**************************************************
!*    CORRECTION FOR LATERAL VELOVITY VARIATIONS
!**************************************************
   do ix = 1,ntr
      nq = (2.0/pi)*atan(q(iz,ix))
      expq = (1. - nq)/2.0
      do iw = iw1,iw2
        ctr = (omega(iw)/wr)**expq
        vomegaxz = (v(iz,ix) * ctr)/2.0
        bomegaxz = omega(iw)/vomegaxz
        aomegaxz = bomegaxz*tan((1.-nq)*(pi/4.))
        komegaxz = cmplx(bomegaxz,aomegaxz)
        zdif = komega(iw) - komegaxz
        phase = cmplx(0.0,-1.0)*zdif
        shift = cexp(phase*dz)
        cp(iw,ix) = cp(iw,ix)*shift
      end do
   end do
!**** Transform X -> Kx
   do iw = iw1,iw2
      do ikx = 1,nx2
  	     rcw2(ikx)   = real(cp(iw,ikx))
	     icw2(ikx)   = imag(cp(iw,ikx))
      end do
      call fft(rcw2,icw2,nx2,nx2,nx2,-1)
      do ikx = 1,nx2
         cp(iw,ikx) = cmplx(rcw2(ikx), icw2(ikx))
      end do
   end do
end do
!*************************
!* END LOOP OVER DEPTHS
!*************************
write(*,'(a)') "> Modeling: direct FFT along X OK."
!**** inverse transform kx -> x
do iw = iw1,iw2
   do ikx = 1,nx2
      rcw2(ikx)   = real(cp(iw,ikx)*(1.0/nx2))
      icw2(ikx)   = imag(cp(iw,ikx)*(1.0/nx2))
   end do
   call fft(rcw2,icw2,nx2,nx2,nx2,1)
   do ikx = 1,nx2
      cp(iw,ikx) = cmplx(rcw2(ikx), icw2(ikx))
   end do
end do
write(*,'(a)') " > Modeling: inverse FFT along X OK."
!**** DONE ...
do ix = 1,ntr
   do iw = 1,iw1
  	  rcw2(iw)   = 0.
	  rcw2(nt2-iw+1) = 0.
	  icw2(iw)   = 0.
	  icw2(nt2-iw+1) = 0.
   end do
   do iw = iw1,iw2
     rcw2(iw)   = real(cp(iw,ix))
     rcw2(nt2-iw+1) = real(cp(iw,ix))
     icw2(iw)   = imag(cp(iw,ix))
	 icw2(nt2-iw+1) = -imag(cp(iw,ix))
   end do
   do iw = iw2+1,nt2
       rcw2(iw)   = 0.
       icw2(iw)   = 0.
   end do
   call fft(rcw2,icw2,nt2,nt2,nt2,-1)
   write(fdo,'(<nt>(f12.8,1x))' ) (rcw2(it),it=1,nt)
end do
close (fdo)
write (*,'(a)') "> Modeling: END."
end
! --------------------------------------------------
subroutine mean_velocity(vmean, v, ntr, nz)
!**** computes mean velocity
implicit none
integer :: nz,ntr,i,j
real(4) :: vmean(nz), v(nz,ntr)
! - - - begin - - -