inverse.cc

二维射线追踪地震层析成像

/*
 * inverse.cc
 * 
 * Jun Korenaga, MIT/WHOI
 * January 1999
 */

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <geom.h>
#include <util.h>
#include "inverse.h"
#include "sparse_rect.h"
#include "lsqr.h"

//
// TomographicInversion2d
//
TomographicInversion2d::TomographicInversion2d(SlownessMesh2d& m, const char* datafn,
					       int xorder, int zorder, double crit_len,
					       int nintp, double cg_tol, double br_tol)

    : smesh(m), graph(smesh,xorder,zorder),
      betasp(1,0,nintp), bend(smesh,betasp,cg_tol,br_tol),
      nrefl(0), nnoded(0), do_full_refl(false),
      nnodev(smesh.numNodes()), nx(smesh.Nx()), nz(smesh.Nz()),
      itermax_LSQR(nnodev*100), LSQR_ATOL(1e-3), 
      smooth_velocity(false), logscale_vel(false), do_filter(false),
      wsv_min(0.0), wsv_max(0.0), dwsv(1.0), weight_s_v(0.0),
      smooth_depth(false), logscale_dep(false), wsd_min(0.0),
      wsd_max(0.0), dwsd(1.0), weight_s_d(0.0), corr_dep_p(0),
      damp_velocity(false), damp_depth(false),
      target_dv(0.0), target_dd(0.0),
      damping_is_fixed(false), weight_d_v(0.0), weight_d_d(0.0),
      do_squeezing(false),
      robust(false), crit_chi(1e3), refl_weight(1.0), target_chisq(0.0),
      printLog(false), verbose_level(-1),
      printFinal(false), printTransient(false), out_mask(false), out_level(0),
      gravity(false), out_grav_dws(false), ngravdata(0)
{
    if (crit_len>0) graph.refineIfLong(crit_len);

    // there can be two interfaces at most (seafloor and moho).
    start_i.reserve(2); end_i.reserve(2);
    interp.reserve(2);
    bathyp = new Interface2d(smesh);
    
    read_file(datafn);

    const int max_np = int(2*sqrt(float(nnodev)));
    path.reserve(max_np);
    pp.reserve(max_np+4);
    
    A.resize(ndata); 
    Rv_h.resize(nnodev); Rv_v.resize(nnodev);
    Tv.resize(nnodev);
    data_vec.resize(ndata);
    tmp_data.resize(ndata);
    total_data_vec.reserve(2*ndata+3*nnodev+2*max_np);
    r_dt_vec.resize(ndata);
    int idata=1;
    for (int i=1; i<=obs_dt.size(); i++){
	for (int j=1; j<=obs_dt(i).size(); j++){
	    r_dt_vec(idata) = 1.0/obs_dt(i)(j);
	    idata++;
	}
    }
    modelv.resize(nnodev);
    mvscale.resize(nnodev);
    dmodel_total.resize(nnodev);
    nnode_total = nnodev;
    path_length.resize(ndata);

    Q.resize(nintp);

    nodev_hit.resize(nnodev); 
    tmp_node.resize(nnodev);
    tmp_nodev.resize(nnodev);
}

void TomographicInversion2d::read_file(const char* datafn)
{
    ifstream in(datafn);
    if (!in){
	cerr << "TomographicInversion2d::cannot open " << datafn << "\n";
	exit(1);
    }

    int iline=0;

    int nsrc;
    in >> nsrc; iline++;
    if (nsrc<=0) error("TomographicInversion2d::invalid nsrc");
    src.resize(nsrc); rcv.resize(nsrc);
    raytype.resize(nsrc); obs_ttime.resize(nsrc); obs_dt.resize(nsrc);
    res_ttime.resize(nsrc);
    
    int isrc=0;
    while(in){
	char flag;
	double x, y;
	int nrcv;

	in >> flag >> x >> y >> nrcv; iline++;
	if (flag!='s'){
	    cerr << "TomographicInversion2d::bad input (s) at l."
		 << iline << '\n';
	    exit(1);
	}
	isrc++;
	src(isrc).set(x,y);

	rcv(isrc).resize(nrcv); raytype(isrc).resize(nrcv);
	obs_ttime(isrc).resize(nrcv); obs_dt(isrc).resize(nrcv);
	res_ttime(isrc).resize(nrcv);
	for (int ircv=1; ircv<=nrcv; ircv++){
	    int n;
	    double t_val, dt_val;
	    in >> flag >> x >> y >> n >> t_val >> dt_val; iline++;
	    if (flag!='r'){
		cerr << "TomographicInversion2d::bad input (r) at l."
		     << iline << '\n';
		exit(1);
	    }
	    if (dt_val < 1e-6){
		cerr << "TomographicInversion2d::bad input (zero obs_dt) at l."
		     << iline << '\n';
		exit(1);
	    }
	    rcv(isrc)(ircv).set(x,y);
	    raytype(isrc)(ircv) = n;
	    obs_ttime(isrc)(ircv) = t_val;
	    obs_dt(isrc)(ircv) = dt_val;
	}
	if (isrc==nsrc) break;
    }
    if (isrc != nsrc) error("TomographicInversion2d::mismatch in nsrc");

    ndata = 0;
    for (int isrc=1; isrc<=nsrc; isrc++){
	ndata += rcv(isrc).size();
    }
}

void TomographicInversion2d::doRobust(double a)
{
    if (a>0){
	robust = true;
	crit_chi = a;
    }
}

void TomographicInversion2d::removeOutliers()
{
    Array1d<double> Adm(ndata_valid);
    SparseRectangular sparseA(A,tmp_data,tmp_node,nnode_total);
    sparseA.Ax(dmodel_total,Adm);
    Array1d<int> tmp_data2(tmp_data.size());

    tmp_data2 = 0;
    int ndata_valid2=0;
    for (int i=1; i<=tmp_data.size(); i++){
	int j=tmp_data(i);
	if (j>0){
	    double res=Adm(j)-data_vec(i);
	    if (abs(res)<=crit_chi){
		tmp_data2(i) = ++ndata_valid2;
	    }
	}
    }

    // redefine data node vector and related stats
    ndata_valid = ndata_valid2;
    tmp_data = tmp_data2;
    rms_tres[0]=rms_tres[1]=0;
    init_chi[0]=init_chi[1]=0;
    ndata_in[0]=ndata_in[1]=0;
    int idata=1;
    for (int isrc=1; isrc<=src.size(); isrc++){
	for (int ircv=1; ircv<=rcv(isrc).size(); ircv++){
	    if (tmp_data(idata)>0){
		double res=res_ttime(isrc)(ircv);
		int icode = raytype(isrc)(ircv);
		double res2=res*r_dt_vec(idata);
		double res22=res2*res2;
		rms_tres[icode] += res*res;
		init_chi[icode] += res22;
		++ndata_in[icode];
	    }
	    idata++;
	}
    }
    rms_tres_total = sqrt((rms_tres[0]+rms_tres[1])/ndata_valid);
    init_chi_total = (init_chi[0]+init_chi[1])/ndata_valid;
    for (int i=0; i<=1; i++){
	rms_tres[i] = ndata_in[i]>0 ? sqrt(rms_tres[i]/ndata_in[i]) : 0.0;
	init_chi[i] = ndata_in[i]>0 ? init_chi[i]/ndata_in[i] : 0.0;
    }
}

void TomographicInversion2d::setLSQR_TOL(double a)
{
    if (a<=0){
	cerr << "TomographicInversion2d::setLSQR_TOL - invalid TOL (ignored)\n";
	return;
    }
    LSQR_ATOL = a;
}

void TomographicInversion2d::setLogfile(const char* fn)
{
    printLog = true;
    log_os_p = new ofstream(fn);
    if (!(*log_os_p)){
	cerr << "TomographicInversion2d::setLogfile - can't open "
	     << fn << '\n';
	exit(1);
    }
}

void TomographicInversion2d::setVerbose(int i){ verbose_level=i; }

void TomographicInversion2d::outStepwise(const char* fn, int i)
{
    printTransient = true;
    out_root = fn;
    out_level = i;
}

void TomographicInversion2d::outFinal(const char* fn, int i)
{
    printFinal = true;
    out_root = fn;
    out_level = i;
}

void TomographicInversion2d::targetChisq(double c)
{
    target_chisq = c;
}

void TomographicInversion2d::outMask(const char* fn)
{
    out_mask = true;
    vmesh_os_p = new ofstream(fn);
    if (!(*vmesh_os_p)){
	cerr << "TomographicInversion2d::outMask - can't open "
	     << fn << '\n';
	exit(1);
    }
}

void TomographicInversion2d::addonGravity(double _z, const Array1d<double>& _x,
					  const Array1d<double>& _g,
					  AddonGravityInversion2d* p, double w,
					  const char *dwsfn)
{
    out_grav_dws = true;
    grav_dws_osp = new ofstream(dwsfn);
    if (!(*grav_dws_osp)){
	cerr << "TomographicInversion2d::addonGravity - can't open "
	     << dwsfn << '\n';
	exit(1);
    }
    addonGravity(_z,_x,_g,p,w);
}

void TomographicInversion2d::addonGravity(double _z, const Array1d<double>& _x,
					  const Array1d<double>& _g,
					  AddonGravityInversion2d* p, double w)
{
    gravity = true;
    grav_z0 = _z;
    ngravdata = _x.size();
    grav_x.resize(ngravdata); grav_x = _x;
    obs_grav.resize(ngravdata); obs_grav = _g;
    res_grav.resize(ngravdata);
    B.resize(ngravdata);
    tmp_gravdata.resize(ngravdata);
    for (int i=1; i<=ngravdata; i++) tmp_gravdata(i) = i;
    ginv = p;
    weight_grav = w;
    if (w<0){
	error("TomographicInversion2d::addonGravity - negative weight detected.");
    }
}

void TomographicInversion2d::solve(int niter)
{
    typedef map<int,double>::iterator mapIterator;
    typedef map<int,double>::const_iterator mapBrowser;
    
    const int bend_nfac=1;
    if (printLog){
	*log_os_p << "# strategy " << jumping << " " << robust << " " << crit_chi << '\n';
	*log_os_p << "# ray_trace " << graph.xOrder() << " " << graph.zOrder() << " "
		  << graph.critLength() << " "
		  << betasp.numIntp() << " " << bend.tolerance() << '\n';
	*log_os_p << "# smooth_vel " << smooth_velocity << " "
		  << wsv_min << " " << wsv_max << " " << dwsv << " "
		  << do_filter << '\n';
	*log_os_p << "# smooth_dep " << smooth_depth << " "
		  << wsd_min << " " << wsd_max << " " << dwsd << '\n';
	if (damping_is_fixed){
	    *log_os_p << "# fixed_damping " << weight_d_v << " " << weight_d_d << '\n';
	}else{
	    *log_os_p << "# damp_vel " << damp_velocity << " " << target_dv << '\n';
	    *log_os_p << "# damp_dep " << damp_depth << " " << target_dd << '\n';
	}
	*log_os_p << "# ndata " << ndata << '\n';
	if (gravity){
	    *log_os_p << "# grav_data " << obs_grav.size() << " " << weight_grav << '\n';
	}
	*log_os_p << "# nnodes " << nnodev << " " << nnoded << " " << refl_weight << '\n';
	*log_os_p << "# LSQR " << LSQR_ATOL << endl;
    }

    // construct index mappers
    for (int i=1; i<=nnodev; i++){
	tmp_node(i) = i;
	tmp_nodev(i) = i;
    }
    if (nrefl>0){
	for (int i=1; i<=nnoded; i++){
	    int i_total = i+nnodev;
	    tmp_node(i+nnodev) = i_total;
	    tmp_nodedc(i) = i_total;
	    tmp_nodedr(i) = i;
	}
    }

    Array1d<int> kstart;
    Array1d<double> filtered_model, dx_vec, dxr_vec, dxn_vec;
    if (do_filter){
	filtered_model.resize(dmodel_total.size());
	dx_vec.resize(dmodel_total.size());
	dxr_vec.resize(dmodel_total.size());
	dxn_vec.resize(dmodel_total.size());
	kstart.resize(nx);
	for (int i=1; i<=nx; i++){
	    Point2d p = smesh.nodePos(smesh.nodeIndex(i,1));
	    double boundz = uboundp->z(p.x());
	    kstart(i) = 1;
	    for (int k=1; k<=nz; k++){
		if (smesh.nodePos(smesh.nodeIndex(i,k)).y()>boundz) break;
		kstart(i) = k;
	    }
	}
    }
    
    // pre-allocate temporary arrays
    Array1d<double> tmp_modelv(modelv.size());
    Array1d<double> tmp_modeld(modeld.size());
    if (jumping) dmodel_total_sum.resize(dmodel_total.size());

    bool isFinal=false;
    int iter=1;
    while(iter<=niter){
	if (verbose_level>=0) cerr << "TomographicInversion2d::iter="
				   << iter << "(" << niter << ")\n";
	if (iter==niter) isFinal=true;
	double graph_time=0.0, bend_time=0.0;
    
	smesh.get(modelv);
	if (nrefl==1) reflp->get(modeld);
	if (iter==1 || !jumping){ // set model scaling vectors
	    mvscale = modelv;
	    if (nrefl==1) mdscale = modeld;
	}

	reset_kernel(); nodev_hit=0;
	int idata=1;
	for (int isrc=1; isrc<=src.size(); isrc++){
	    if (verbose_level>=0){
		cerr << "isrc=" << isrc << " nrec=" << rcv(isrc).size()
		     << " ";
	    }
	    ofstream *tres_os_p, *ray_os_p;
	    if (printTransient || (printFinal && isFinal)){
		if (out_level >= 1){
		    char transfn[MaxStr];
		    sprintf(transfn, "%s.tres.%d.%d", out_root, iter, isrc);
		    tres_os_p = new ofstream(transfn);
		}
		if (out_level >= 2){
		    char transfn[MaxStr];
		    sprintf(transfn, "%s.ray.%d.%d", out_root, iter, isrc);
		    ray_os_p = new ofstream(transfn);
		}
	    }

	    // limit range first
	    double xmin=smesh.xmax();
	    double xmax=smesh.xmin();
	    for (int ircv=1; ircv<=rcv(isrc).size(); ircv++){
		double x = rcv(isrc)(ircv).x();
		if (x < xmin) xmin = x;
		if (x > xmax) xmax = x;
	    }
	    double srcx=src(isrc).x();
	    if (srcx < xmin) xmin = srcx;
	    if (srcx > xmax) xmax = srcx;
	    graph.limitRange(xmin,xmax);
	    if (verbose_level>=1){
		cerr << "xrange=(" << xmin << "," << xmax << ") ";
	    }

	    //
	    if (do_full_refl){
		double pwater = 1.0/1.5;
		tmp_modelv = modelv;
		Array1d<int> irefl(nx);
		for (int i=1; i<=nx; i++){
		    Point2d p = smesh.nodePos(smesh.nodeIndex(i,1));
		    double reflz = reflp->z(p.x());
		    irefl(i) = nz;
		    for (int k=1; k<=nz; k++){
			if (smesh.nodePos(smesh.nodeIndex(i,k)).y()>reflz){
			    irefl(i) = k;
			    break;
			}
		    }
		}
		for (int i=1; i<=nx; i++){
		    for (int k=2; k<=nz; k++){
			if (k==irefl(i)){
			    // this is to make the velocity at the node just below the reflector to
			    // be the same sa the velocity just above. If I use pwater for this node,
			    // there would be a unwanted low velocity gradient surrounding the reflector.
			    tmp_modelv(smesh.nodeIndex(i,k)) = tmp_modelv(smesh.nodeIndex(i,k-1));
			}else if (k>irefl(i)){
			    tmp_modelv(smesh.nodeIndex(i,k)) = pwater;
			}
		    }
		}
	    }
	    
	    clock_t start_t=clock();
	    graph.solve(src(isrc));
	    bool is_refl_solved = false;
	    clock_t end_t=clock();
	    graph_time += end_t-start_t;

	    start_t = clock();
	    double graph_refl_time=0;
	    for (int ircv=1; ircv<=rcv(isrc).size(); ircv++){
		Point2d r = rcv(isrc)(ircv);
		double orig_t, final_t;
		int iterbend;
		int icode = raytype(isrc)(ircv);
		if (icode == 0){ // refraction
		    if (smesh.inWater(r)){
			if (verbose_level>=0) cerr << "*";
			int i0, i1;
			graph.pickPathThruWater(r,path,i0,i1);
			start_i.resize(1); end_i.resize(1); interp.resize(1);
			start_i(1) = i0; end_i(1) = i1; interp(1) = bathyp;
			iterbend=bend.refine(path,orig_t,final_t,
					     start_i, end_i, interp);
			if (verbose_level>=1) cerr << "(" << iterbend << ")";
		    }else{
			if (verbose_level>=0) cerr << ".";
			graph.pickPath(rcv(isrc)(ircv),path);
			iterbend=bend.refine(path,orig_t,final_t,bend_nfac);
			if (verbose_level>=1) cerr << "(" << iterbend << ")";
		    }
		    if (iterbend<0){
			cerr << "TomographicInversion2d::iteration = " << iter << '\n';
			cerr << "TomographicInversion2d::too many iterations required\n";
			cerr << "TomographicInversion2d::for bending refinement at (s,r)="
			     << isrc << "," << ircv << '\n';
			exit(1);
		    }
		    add_kernel(idata,path);
		}else if (icode == 1){ // reflection
		    if (nrefl==0){
			error("TomographicInversion2d:: reflector not specified.");
		    }
		    if (do_full_refl){
			// temporarily replace sub-reflector velocity field
			// with water velocity
			smesh.set(tmp_modelv);
		    }
		    if (!is_refl_solved){
			clock_t start_t=clock();
			graph.solve_refl(src(isrc), *reflp);
			clock_t end_t=clock();
			graph_refl_time = end_t-start_t;
			graph_time += graph_refl_time;
			is_refl_solved = true;
		    }
		    int ir0, ir1;
		    if (smesh.inWater(r)){
			if (verbose_level>=0) cerr << "#";
			int i0, i1;
			graph.pickReflPathThruWater(r,path,i0,i1,ir0,ir1);
			start_i.resize(2); end_i.resize(2); interp.resize(2);
			start_i(1) = i0; end_i(1) = i1; interp(1) = bathyp;
			start_i(2) = ir0; end_i(2) = ir1; interp(2) = reflp;