[转载]利用fk批量合成理论地震记录----一事件对多台站

原文地址：利用fk批量合成理论地震记录----一事件对多台站作者：5Geo

在做台阵数据处理过程中需要验证某个台阵处理方法如beamforming，frequency-wavenumber,vespa-stack等方法，本来想用一个简单的雷克子波或者正弦波去做，但为了更贴近实际，用fk合成理论地震记录势必要信服点，以下来那个个子程序用于调用Lupei Zhu老师写的fk程序来得到合成地震记录。

fk合成地震记录需要我们给的地震到台站的方位角以及震中距，这个必须得我们自己算，所以以下fortran代码用于产生一个地震到多个台站的理论地震记录，需要台站的经纬度以及事件的经纬度信息。

代码见

http://pan.baidu.com/s/1dFsMXvr

如果合成的记录偏短，想对其进行补零操作，应用sac中的 命令

sac << eof

cuterr fillz

cut 0 e

r  *

w over

q

eof

这是将所有的到头b改成0，并且长度为0-e,其实e也可以改成比e更大的数，大于e的话则后面记录补零

这个代码主要产生震中距和方位角文件，后续代码用python调用fk

program main

implicit none

integer,parameter::nst=1000

character(len= 130) line

character(len = 60) staname(nst)

real    sta_lat(nst),sta_lon(nst),sta_elev(nst)

            real    sta_lat1,sta_lon1,sta_elev1

real    xout,yout,zout,theta

real    sta_x(nst),sta_y(nst),sta_z(nst)

real    dist

            real    az

real    del

real    elat,elon

integer     i,nsta

        open(10,file='syn_input.inc')

        read(10,*) nsta

        read(10,*) elat,elon

        if(nsta.gt.nst) stop 'nsta is too large'

        do i =1 ,nsta

read(10,'(a)') line 

read(line,*) staname(i),sta_lat(i),sta_lon(i),sta_elev(i)

 enddo

        close(10)   

! USED FOR syn in fk calculating the waveform

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcCCCCCCCCCCCCCCC  used for syn        

 open(11,file='azi_dist.info')

        write(*,*) "staname----sta_lat-------sta_lon------elat-------elon---------dist------az"

        do i= 1 , nsta

sta_lat1 =sta_lat(i)

sta_lon1= sta_lon(i)

sta_elev1=sta_elev(i)

        call delaz(elat,elon,sta_lat1,sta_lon1,del,dist,az)

        write(11,*) trim(staname(i)),sta_lat1,sta_lon1,elat,elon,dist,az

        write(*,*) trim(staname(i)),sta_lat1,sta_lon1,elat,elon,dist,az

        enddo

        close(11)

        end

!c Compute distance and azimuth on a sphere

subroutine delaz(alat, alon, blat, blon, del, dist, az)

!c computes distance and azimuth from a to b

!c a and b are in decimal degrees and n-e coordinates

!c del -- delta in degrees

!c dist -- distance in km

!c az -- azimuth from a to b clockwise from north in degrees

!c Original author:  Bruce Julian

!c new values for radius used (from Edi Kissling)

implicit none

!c Parameters:

real alat, alon ! Coordinates of first point

real blat, blon ! Coordinates of second point

real del ! Central angle (degrees)

real dist ! Distance (km)

real az ! Azimuth from a to b (degrees)

real xtop

real xden

!c Local variables:

real*8 acol, bcol

real*8 azr

real*8 blatr, blonr

real*8 colat

real*8 cosdel

real*8 delr

real*8 flat

real*8 geoa

real*8 geob

real*8 rad

real*8 radius

real*8      alatr, alonr

real*8      diflon

real*8      pi2

real*8      tana, tanb

!c Built-in functions:  Declarations not needed

real*8      dtan

real*8 datan

real*8 dsin

real*8 dcos

real*8 dacos

!c doubleprecision top,den

! character rcsid*150

! save rcsid

data pi2/1.570796d0/

data rad/1.745329d-02/

data flat/.993231d0/

!c-----convert to radians

alatr=alat*rad

alonr=alon*rad

blatr=blat*rad

blonr=blon*rad

!c-----convert latitudes to geocentric colatitudes

tana=flat*dtan(alatr)

geoa=datan(tana)

acol=pi2-geoa

tanb=flat*dtan(blatr)

geob=datan(tanb)

bcol=pi2-geob

!c-----calculate delta

diflon=blonr-alonr

cosdel=dsin(acol)*dsin(bcol)*dcos(diflon)+dcos(acol)*  &

      dcos(bcol)

delr=dacos(cosdel)

!c-----calculate azimuth from a to b

!c***** Note the use of single precision xtop and xden instead

!c of the double precision top and den in the original

!c program.

!c***** Note also the call to atan2 instead of datan2.

!c Both of these changes were made so that dyn.load

!c would work in Splus.  For some reason, the ld command

!c ld -r -d didn't find _d_atan2

!c WLE 10/16/91

xtop = dsin(diflon)

xden=(dsin(acol)/dtan(bcol))-dcos(diflon)*dcos(acol)

azr=atan2(xtop,xden)

!c----- convert to degrees

del=delr/rad

az=azr/rad

if(az.lt.0.0) az=360.+az

!c-----compute distance in kilometers

colat=pi2-(alatr+blatr)/2.d0

radius=6378.163d0*   &

      (1.d0+3.35278d-3*((1.d0/3.d0)-(dcos(colat)**2)))

dist=delr*radius

return

!c  ***** end of subroutine delaz *****

endsubroutine

#python代码批量合成地震记录

import os

import string

# read the array position

#array_file = open('array.station.info','r')

array_file = open('azi_dist.info','r')

line=array_file.readlines()

stname=[]

temp=[]

stlat=[]

stlon=[]

stelev=[]

total=[]

az = []

dist = []

for data in line:

    temp =  data.split()

    total.append(temp)

    stname.append(temp[0])

    stlat.append(temp[1])

    stlon.append(temp[2])

    dist.append(temp[5])

    az.append(temp[6])

i=0

while i

    dist[i]=str(dist[i])

    i=i+1

os.environ['dist']=str(" ".join(dist))

os.system('echo ${dist}')

os.environ['depth']="1"   # 修改震源深度 km

os.system('echo ${depth}')

os.system('./fk/fk.pl -Mhk/${depth}/k -N1024/0.01 ${dist}')  # 根据hk模型合成格林函数，存储在hk文件夹中，hk为模型文件在fk文件夹下

for azimuth in az:

    index=az.index(azimuth)

    os.environ['dist']=str(dist[index])

    name=stname[index]

    os.environ['stn'] =str(name[0:3])

    os.environ['azimuth']=str(azimuth)

    os.system('./fk/syn   -M5/335/80/-70 -D0.5 -A${azimuth} -OWave/${stn}.Z -Ghk_${depth}/${dist}.grn.0')

#默认合成5级地震，修改命令行请见fk文件夹中的syn.c -OWave表示输出文件夹

fk合成地震记录的用法见fk下的fk.pl以及syn.c代码打印的信息。这里调用fk.pl需要将fk.pl脚本的绝对路径或者相对路径给上，hk速度模型必须在当前合成地震记录的文件夹下，不是在fk本身文件夹下，否则会报错"can not open hk",需要事先建立波形输出的文件夹 Wave。

笔记比较杂乱,有事请留言！