inverse.cc

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

			int 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.pickReflPath(r,path,ir0,ir1);
			start_i.resize(1); end_i.resize(1); interp.resize(1);
			start_i(1) = ir0; end_i(1) = ir1; interp(1) = reflp;
			int iterbend=bend.refine(path,orig_t,final_t,
						 start_i, end_i, interp);
			if (verbose_level>=1) cerr << "(" << iterbend << ")";
		    }
		    if (do_full_refl){
			// revert to the original smesh
			smesh.set(modelv);
		    }
		    add_kernel_refl(idata,path,ir0,ir1);
		}else{
		    error("TomographicInversion2d:: illegal raycode detected.");
		}
		res_ttime(isrc)(ircv)=obs_ttime(isrc)(ircv)-final_t;
		idata++;

		if (printTransient || (printFinal && isFinal)){
		    if (out_level>=1){
			*tres_os_p << rcv(isrc)(ircv).x() << " "
				   << res_ttime(isrc)(ircv) << '\n';
		    }
		    if (out_level>=2){
			*ray_os_p << ">\n";
			printCurve(*ray_os_p, path, betasp);
		    }
		}
	    }
	    if (printTransient || (printFinal && isFinal)){
		if (out_level>=1){ tres_os_p->flush(); delete tres_os_p;}
		if (out_level>=2){ ray_os_p->flush(); delete ray_os_p;}
	    }
		    
	    if (verbose_level>=0) cerr << '\n';
	    end_t = clock();
	    bend_time += end_t-start_t-graph_refl_time;
	}
	graph_time /= CLOCKS_PER_SEC;
	bend_time /= CLOCKS_PER_SEC;

	// construct data vector
	idata=1;
	rms_tres[0]=rms_tres[1]=0;
	init_chi[0]=init_chi[1]=0;
	ndata_in[0]=ndata_in[1]=0;
	ndata_valid=0;
	tmp_data=0;
	for (int isrc=1; isrc<=src.size(); isrc++){
	    for (int ircv=1; ircv<=rcv(isrc).size(); ircv++){
		double res=res_ttime(isrc)(ircv);
		data_vec(idata) = res;
		double res2=res*r_dt_vec(idata);
		double res22=res2*res2;
		int icode = raytype(isrc)(ircv);
		rms_tres[icode] += res*res;
		init_chi[icode] += res22;
		++ndata_in[icode];
		tmp_data(idata) = ++ndata_valid;
		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;
	}
	if (init_chi_total<target_chisq) isFinal=true;

	// rescale kernel A and data vector
	// note: averaging matrices are scaled upon their construction.
	for (int i=1; i<=ndata; i++){
	    double data_wt=r_dt_vec(i);
	    data_vec(i) *= data_wt;
	    for (mapIterator p=A(i).begin(); p!=A(i).end(); p++){
		int j = p->first;
		double m;
		if (j<=nnodev){
		    m = mvscale(j);
		}else{
		    m = mdscale(j-nnodev)*refl_weight;
		}
		p->second *= m*data_wt;
	    }
	}

	// construct total kernel matrix and data vecter, and solve Ax=b
	if (smooth_velocity) calc_averaging_matrix();
	if (nrefl>0 && smooth_depth) calc_refl_averaging_matrix();
	if (damp_velocity) calc_damping_matrix();
	if (nrefl>0 && damp_depth) calc_refl_damping_matrix();

	// (optional) joint inversion with gravity anomalies
	if (gravity){
	    if (verbose_level>=0) cerr << "calculating residual gravity anomalies...";
	    ginv->calcGravity(grav_z0, grav_x, res_grav);
//	    for (int i=1; i<=res_grav.size(); i++){
//		cerr << grav_x(i) << " " << res_grav(i) << '\n';
//	    }
	    if (verbose_level>=0) cerr << "done.\n";

	    if (verbose_level>=0) cerr << "calculating gravity kernel...";
	    rms_grav = 0.0;
	    for (int i=1; i<=ngravdata; i++){
		if (verbose_level>=1) cerr << i << " ";
		// residual gravity anomalies
		double rg = obs_grav(i)-res_grav(i);
		res_grav(i) = rg*weight_grav;
		rms_grav += rg*rg;

		// gravity kernels
		ginv->calcGravityKernel(grav_z0, grav_x(i), B(i));
		// model normalization
		for (mapIterator p=B(i).begin(); p!=B(i).end(); p++){
		    int j = p->first;
		    double m;
		    if (j<=nnodev){
			m = mvscale(j);
		    }else{
			m = mdscale(j-nnodev)*refl_weight;
		    }
		    p->second *= m*weight_grav;
		}
	    }
	    rms_grav = sqrt(rms_grav/ngravdata);
	    if (verbose_level>=0) cerr << "done.\n";

	    if (printTransient || (printFinal && isFinal)){
		char transfn[MaxStr];
		sprintf(transfn, "%s.rgrav.%d", out_root, iter);
		ofstream os(transfn);
		for (int i=1; i<=ngravdata; i++){
		    os << grav_x(i) << " " << res_grav(i) << " " << res_grav(i)/weight_grav << '\n';
		}
	    }
	}
	
	int iset=0;
	for (double tmp_wsv=wsv_min; tmp_wsv<=wsv_max; tmp_wsv+=dwsv){
	    for (double tmp_wsd=wsd_min; tmp_wsd<=wsd_max; tmp_wsd+=dwsd){
		iset++;

		weight_s_v = logscale_vel ? pow(10.0,tmp_wsv) : tmp_wsv;
		weight_s_d = logscale_dep ? pow(10.0,tmp_wsd) : tmp_wsd;
		
		double wdv,wdd;
		int nlsqr=0, lsqr_iter=0;
		time_t start_t = time(NULL);
		if (damping_is_fixed){
		    fixed_damping(lsqr_iter,nlsqr,wdv,wdd);
		}else{
		    auto_damping(lsqr_iter,nlsqr,wdv,wdd);
		}
		double lsqr_time = difftime(time(NULL),start_t);

		if (do_filter){
		    if (verbose_level>=0) cerr << "filtering velocity perturbation...";
		    for (int i=1; i<=dmodel_total.size(); i++) filtered_model(i) = dmodel_total(i);
		    for (int i=1; i<=nx; i++){
			for (int k=kstart(i); k<=nz; k++){
			    int inode = smesh.nodeIndex(i,k);
			    Point2d p = smesh.nodePos(inode);
			    double Lh, Lv;
			    corr_vel_p->at(p,Lh,Lv);

			    double Lh2 = Lh*Lh;
			    double Lv2 = Lv*Lv;
			    double sum=0.0, beta_sum=0.0;
			    for (int ii=1; ii<=nx; ii++){
				int jnode = smesh.nodeIndex(ii,k);
				double dx = smesh.nodePos(jnode).x()-p.x();
				if (abs(dx)<=Lh){
				    double xexp = exp(-dx*dx/Lh2);
				    for (int kk=kstart(ii); kk<=nz; kk++){
					int knode = smesh.nodeIndex(ii,kk);
					double dz = smesh.nodePos(knode).y()-p.y();
					if (abs(dz)<=Lv){
					    double beta = xexp*exp(-dz*dz/Lv2);
					    sum += beta*dmodel_total(knode);
					    beta_sum += beta;
					}
				    }
				}
			    }
			    filtered_model(inode) = sum/beta_sum;
			}
		    }
		    if (verbose_level>=0) cerr << "done.\n";
		    
		    // conservative filtering (Deal and Nolet, GJI, 1996)
		    if (verbose_level>=0) cerr << "calculating a nullspace shuttle...";
		    dx_vec=0.0;
		    for (int i=1; i<=nnodev; i++) dx_vec(i) = filtered_model(i)-dmodel_total(i);
		    Array1d<double> Adm(ndata_valid);
		    SparseRectangular sparseA(A,tmp_data,tmp_node,nnode_total);
		    sparseA.Ax(dx_vec,Adm);

		    int lsqr_itermax=itermax_LSQR;
		    double chi;
		    dxr_vec=0.0;
		    iterativeSolver_LSQR(sparseA,Adm,dxr_vec,LSQR_ATOL,lsqr_itermax,chi);

		    // note: correction_vec for depth nodes should be zero, so I don't use them.
		    double orig_norm=0, new_norm=0, iprod=0;
		    for (int i=1; i<=nnodev; i++){
			orig_norm += filtered_model(i)*filtered_model(i);
			dmodel_total(i) = filtered_model(i)-dxr_vec(i);
			new_norm += dmodel_total(i)*dmodel_total(i);
			dxn_vec(i) = dx_vec(i)-dxr_vec(i);
			iprod += dxr_vec(i)*dxn_vec(i);
		    }
		    if (verbose_level>=0) cerr << "done.\n";
		    if (printLog){
			*log_os_p << "# a posteriori filter check: "
				  << sqrt(orig_norm/nnodev)
				  << " " << sqrt(new_norm/nnodev) 
				  << " " << iprod << endl;
		    }
		}
		
		// take stats
		double pred_chi = calc_chi();
		double dv_norm = calc_ave_dmv();
		double dd_norm = calc_ave_dmd();
		double Lmvh=-1, Lmvv=-1, Lmd=-1;
		calc_Lm(Lmvh,Lmvv,Lmd);

		if (jumping) dmodel_total_sum += dmodel_total;

		// scale to slowness perturbation
		tmp_modelv = modelv;
		tmp_modeld = modeld;
		for (int i=1; i<=nnodev; i++){
		    tmp_modelv(i) += mvscale(i)*dmodel_total(i);
		}
		smesh.set(tmp_modelv);
		if (nrefl>0){
		    for (int i=1; i<=nnoded; i++){
			int tmp_i = tmp_nodedc(i);
			tmp_modeld(i) += mdscale(i)*dmodel_total(tmp_i)*refl_weight;
		    }
		    reflp->set(tmp_modeld);
		}

		if (printTransient || (printFinal && isFinal)){
		    char transfn[MaxStr];
		    sprintf(transfn, "%s.smesh.%d.%d", out_root, iter, iset);
		    ofstream os(transfn);
		    smesh.outMesh(os);

		    if (nrefl>0){
			char transfn2[MaxStr];
			sprintf(transfn2, "%s.refl.%d.%d", out_root, iter, iset);
			ofstream os2(transfn2);
			os2 << *reflp;
		    }
		}
		if (printLog){
		    *log_os_p << iter << " " << iset << " "
			      << ndata-ndata_valid << " "
			      << rms_tres_total << " " << init_chi_total << " "
			      << ndata_in[0] << " " << rms_tres[0] << " " << init_chi[0] << " "
			      << ndata_in[1] << " " << rms_tres[1] << " " << init_chi[1] << " "
			      << graph_time << " " << bend_time << " " 
			      << weight_s_v << " " << weight_s_d << " "
			      << wdv << " " << wdd << " "
			      << nlsqr << " " << lsqr_iter << " " << lsqr_time << " "
			      << pred_chi << " " << dv_norm << " " << dd_norm << " "
			      << Lmvh << " " << Lmvv << " " << Lmd;
		    if (gravity){
			*log_os_p << " " << rms_grav;
		    }
		    *log_os_p << endl;
		}
	    }
	}

	if (isFinal) break;
	iter++;
    }

    // output DWS for the last iteration
    if (out_mask){
	dws.resize(nnodev);
	dws=0.0;
	typedef map<int,double>::iterator mapIterator;
	for (int i=1; i<=A.size(); i++){
	    for (mapIterator p=A(i).begin(); p!=A(i).end(); p++){
		int inode=p->first;
		if (inode<=dws.size()) dws(inode) += p->second;
	    }
	}
	smesh.printMaskGrid(*vmesh_os_p, dws);
    }
    if (out_grav_dws){
	dws.resize(nnodev);
	dws=0.0;
	typedef map<int,double>::iterator mapIterator;
	for (int i=1; i<=B.size(); i++){
	    for (mapIterator p=B(i).begin(); p!=B(i).end(); p++){
		int inode=p->first;
		if (inode<=dws.size()) dws(inode) += p->second/(weight_grav*mvscale(inode));
	    }
	}
	smesh.printMaskGrid(*grav_dws_osp, dws);
    }
}

void
TomographicInversion2d::addRefl(Interface2d* intfp)
{
    reflp = intfp;
    nrefl = 1;
    nnoded = reflp->numNodes();
    Rd.resize(nnoded);
    Td.resize(nnoded);
    modeld.resize(nnoded);
    mdscale.resize(nnoded);
    dmodel_total.resize(nnodev+nnoded);
    nnode_total = nnodev+nnoded;
    tmp_node.resize(nnodev+nnoded);
    tmp_nodedc.resize(nnoded);
    tmp_nodedr.resize(nnoded);
}

void TomographicInversion2d::doFullRefl()
{
    do_full_refl = true;
    graph.do_refl_downward();
}

void TomographicInversion2d::setReflWeight(double x)
{
    if (x>0){
	refl_weight = x;
    }else{
	cerr << "TomographicInversion2d::setReflWeight - non-positive value ignored.\n";
    }
}

void TomographicInversion2d::SmoothVelocity(const char* fn,
					    double start, double end, double d,
					    bool scale)
{
    smooth_velocity=true;
    wsv_min=start; wsv_max=end; dwsv=d;
    logscale_vel=scale;
    corr_vel_p = new CorrelationLength2d(fn);
}

void TomographicInversion2d::applyFilter(const char *fn)
{
    do_filter=true;
    if (fn[0] != '\0'){
	uboundp = new Interface2d(fn);
    }else{
	uboundp = new Interface2d(smesh); // use bathymetry as upperbound
    }
}

void TomographicInversion2d::SmoothDepth(const char* fn,
					 double start, double end, double d,
					 bool scale)
{
    SmoothDepth(start,end,d,scale);
    corr_dep_p = new CorrelationLength1d(fn);
}

void TomographicInversion2d::SmoothDepth(double start, double end, double d,
					 bool scale)
{
    smooth_depth=true;
    wsd_min=start; wsd_max=end; dwsd=d;
    logscale_dep=scale;
}

void TomographicInversion2d::DampVelocity(double a)
{
    damp_velocity = true;
    target_dv = a/100.0; // % -> fraction
}

void TomographicInversion2d::DampDepth(double a)
{
    damp_depth = true;
    target_dd = a/100.0; // % -> fraction 
}

void TomographicInversion2d::FixDamping(double v, double d)
{
    damping_is_fixed = true;
    damp_velocity = true;
    damp_depth = true;
    weight_d_v = v;
    weight_d_d = d;
}

void TomographicInversion2d::Squeezing(const char* fn)
{
    damping_wt_p = new DampingWeight2d(fn);
    do_squeezing = true;
}

void TomographicInversion2d::doJumping()
{
    jumping = true;
}

void TomographicInversion2d::reset_kernel()
{
    typedef map<int,double>::iterator mapIterator;
    
    for (int i=1; i<=A.size(); i++){
	for (mapIterator p=A(i).begin(); p!=A(i).end(); p++) A(i).erase(p);
    }
    for (int i=1; i<=Rv_h.size(); i++){
	for (mapIterator p=Rv_h(i).begin(); p!=Rv_h(i).end(); p++) Rv_h(i).erase(p);
    }
    for (int i=1; i<=Rv_v.size(); i++){
	for (mapIterator p=Rv_v(i).begin(); p!=Rv_v(i).end(); p++) Rv_v(i).erase(p);
    }
    for (int i=1; i<=Tv.size(); i++){
	for (mapIterator p=Tv(i).begin(); p!=Tv(i).end(); p++) Tv(i).erase(p);
    }
    if (nrefl>0){
	for (int i=1; i<=Rd.size(); i++){
	    for (mapIterator p=Rd(i).begin(); p!=Rd(i).end(); p++) Rd(i).erase(p);
	}
	for (int i=1; i<=Td.size(); i++){
	    for (mapIterator p=Td(i).begin(); p!=Td(i).end(); p++) Td(i).erase(p);
	}
    }
    if (gravity){
	for (int i=1; i<=B.size(); i++){
	    for (mapIterator p=B(i).begin(); p!=B(i).end(); p++) B(i).erase(p);
	}
    }
}

void TomographicInversion2d::add_kernel(int idata, const Array1d<Point2d>& cur_path)
{
    map<int,double>& A_i = A(idata);

    int np = cur_path.size();
    int nintp = betasp.numIntp();
    makeBSpoints(cur_path,pp);

    double path_len=0.0;
    Index2d guess_index = smesh.nodeIndex(smesh.nearest(*pp(1)));
    for (int i=1; i<=np+1; i++){
	int j1=i;
	int j2=i+1;
	int j3=i+2;
	int j4=i+3;
	betasp.interpolate(*pp(j1),*pp(j2),*pp(j3),*pp(j4),Q);
	
	for (int j=2; j<=nintp; j++){
	    Point2d midp = 0.5*(Q(j-1)+Q(j));
	    double dist= Q(j).distance(Q(j-1));
	    path_len+=dist;

	    int jUL, jLL, jLR, jUR;
	    double r, s, rr, ss;
	    int icell=smesh.locateInCell(midp,guess_index,
					 jUL,jLL,jLR,jUR,r,s,rr,ss);
	    if (icell>0){
		A_i[jUL] += rr*ss*dist;
		A_i[jLL] += rr*s*dist;
		A_i[jLR] += r*s*dist;
		A_i[jUR] += r*ss*dist;

		nodev_hit(jUL)++;
		nodev_hit(jLL)++;
		nodev_hit(jLR)++;
		nodev_hit(jUR)++;
	    }
	}
    }
    path_length(idata)=path_len;
}

void TomographicInversion2d::add_kernel_refl(int idata, const Array1d<Point2d>& cur_path,
					     int ir0, int ir1)
{
    map<int,double>& A_i = A(idata);

    int jL, jR;
    reflp->locateInSegment(path(ir0).x(),jL,jR);
    if (jL==jR){
	cerr << "TomographicInversion2d::add_kernel_refl - bottoming point out of bounds\n";
	return; // ignore this ray info
    }
    double x1 = reflp->x(jL);
    double x2 = reflp->x(jR);
    double z1 = reflp->z(x1);