inverse.cc

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

    double z2 = reflp->z(x2);
    double x = path(ir0).x();

    double dx = x2-x1;
    double dz = z2-z1;
    double cos_alpha = dx/sqrt(dx*dx+dz*dz);

    double p1 = smesh.at(path(ir0));

    Point2d a = path(ir0-1)-path(ir0); a /= a.norm();
    Point2d b = path(ir1+1)-path(ir1); b /= b.norm();
    Point2d c = a+b; c /= c.norm();
    double cos_theta = a.inner_product(c);
    
    double com_factor = 2.0*cos_theta*cos_alpha*p1/dx;
    A_i[nnodev+jL] = com_factor*(x2-x);
    A_i[nnodev+jR] = com_factor*(x-x1);

    add_kernel(idata,path);
}    

void TomographicInversion2d::calc_averaging_matrix()
{
    for (int i=1; i<=nx; i++){
	for (int k=1; 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 beta_sum = 0.0;
	    for (int ii=1; ii<=nx; ii++){
		int jnode = smesh.nodeIndex(ii,k);
		if (jnode!=inode){
		    double dx = smesh.nodePos(jnode).x()-p.x();
		    if (abs(dx)<=Lh){
			double dxL2 = dx*dx/Lh2;
			double beta = exp(-dxL2);
			Rv_h(inode)[jnode] = beta*mvscale(jnode);
			beta_sum += beta;
		    }
		}
	    }
	    Rv_h(inode)[inode] = -beta_sum*mvscale(inode);

	    beta_sum = 0.0;
	    for (int kk=1; kk<=nz; kk++){
		int jnode = smesh.nodeIndex(i,kk);
		if (jnode!=inode){
		    double dz = smesh.nodePos(jnode).y()-p.y();
		    if (abs(dz)<=Lv){
			double dzL2 = dz*dz/Lv2;
			double beta = exp(-dzL2);
			Rv_v(inode)[jnode] = beta*mvscale(jnode);
			beta_sum += beta;
		    }
		}
	    }
	    Rv_v(inode)[inode] = -beta_sum*mvscale(inode);
	}
    }
}

void TomographicInversion2d::calc_refl_averaging_matrix()
{
    for (int i=1; i<=nnoded; i++){
	double x = reflp->x(i);
	double Lh;
	if (corr_dep_p != 0){
	    Lh = corr_dep_p->at(x);
	}else if (corr_vel_p != 0){
	    double Lv;
	    corr_vel_p->at(Point2d(x,reflp->z(x)),Lh,Lv);
	}else{
	    error("TomographicInversion2d::calc_refl_averaging_matrix - no correlation length available");
	}
	double Lh2 = Lh*Lh;
	double beta_sum = 0.0;
	for (int ii=1; ii<=nnoded; ii++){
	    if (ii!=i){
		double dx = reflp->x(ii)-x;
		if (abs(dx)<=Lh){
		    double beta = exp(-dx*dx/Lh2);
		    Rd(i)[ii] = beta*mdscale(ii)*refl_weight;
		    beta_sum += beta;
		}
	    }
	}
	Rd(i)[i] = -beta_sum*mdscale(i)*refl_weight;
    }
}

void TomographicInversion2d::calc_damping_matrix()
{
    Array2d<double> local_T(4,4);
    Array1d<int> j(4);
    
    for (int i=1; i<=smesh.numCells(); i++){
	smesh.cellNodes(i, j(1), j(2), j(3), j(4));
	smesh.cellNormKernel(i, local_T);

	if (do_squeezing){
	     Point2d p = 0.25*(smesh.nodePos(j(1))+smesh.nodePos(j(2))
			       +smesh.nodePos(j(3))+smesh.nodePos(j(4)));
	     double w;
	     damping_wt_p->at(p,w);
	     local_T *= w;
	}
	add_global(local_T, j, Tv);
    }
}

void TomographicInversion2d::calc_refl_damping_matrix()
{
    double dx = reflp->x(2)-reflp->x(1); // assumes uniform interval
    double fac = sqrt(dx);

    for (int i=1; i<=nnoded; i++){
	Td(i)[i] = fac;
    }
}

void TomographicInversion2d::add_global(const Array2d<double>& a,
					const Array1d<int>& j, sparseMat& global)
{
    for (int m=1; m<=4; m++){
	for (int n=1; n<=4; n++){
	    global(j(m))[j(n)] += a(m,n);
	}
    }
}

double TomographicInversion2d::calc_ave_dmv()
{
    double dm_norm=0.0;
    for (int i=1; i<=nnodev; i++){
	dm_norm += dmodel_total(i)*dmodel_total(i);
    }
    return sqrt(dm_norm/nnodev);
}

double TomographicInversion2d::calc_ave_dmd()
{
    if (nrefl==0) return 0.0;
	
    double dm_norm=0.0;
    for (int i=1+nnodev; i<=nnoded+nnodev; i++){
	dm_norm += dmodel_total(i)*dmodel_total(i);
    }
    return sqrt(dm_norm/nnoded)*refl_weight;
}

void TomographicInversion2d::calc_Lm(double& lmh, double& lmv, double& lmd)
{
    if (smooth_velocity){
	Array1d<double> Rvm(nnodev), dmv_vec(nnodev);
	SparseRectangular sparseRv_h(Rv_h,tmp_nodev,tmp_nodev);
	for (int i=1; i<=nnodev; i++) dmv_vec(i) = dmodel_total(i);

	sparseRv_h.Ax(dmv_vec,Rvm);
	double val=0.0;
	for (int i=1; i<=nnodev; i++) val += Rvm(i)*Rvm(i);
	lmh = sqrt(val/nnodev);

	SparseRectangular sparseRv_v(Rv_v,tmp_nodev,tmp_nodev);
	sparseRv_v.Ax(dmv_vec,Rvm);
	val = 0.0;
	for (int i=1; i<=nnodev; i++) val += Rvm(i)*Rvm(i);
	lmv = sqrt(val/nnodev);
    }

    if (nrefl>0 && smooth_depth){
	Array1d<double> Rdm(nnoded), dmd_vec(nnoded);
	SparseRectangular sparseRd(Rd,tmp_nodedr,tmp_nodedr);
	for (int i=1+nnodev, j=1; i<=nnoded+nnodev; i++) dmd_vec(j++) = dmodel_total(i);

	sparseRd.Ax(dmd_vec,Rdm);
	double val = 0.0;
	for (int i=1; i<=nnoded; i++) val += Rdm(i)*Rdm(i);
	lmd = sqrt(val/nnoded);
    }
}
 
double TomographicInversion2d::calc_chi()
{
    Array1d<double> Adm(ndata_valid);
    SparseRectangular sparseA(A,tmp_data,tmp_node,nnode_total);
    sparseA.Ax(dmodel_total,Adm);
    double val=0.0;
    for (int i=1; i<=ndata; i++){
	int j=tmp_data(i);
	if (j>0){
	    double res=Adm(j)-data_vec(i);
	    val += res*res;
	}
    }
    return val/ndata_valid;
}

int TomographicInversion2d::_solve(bool sv, double wsv, bool sd, double wsd,
				   bool dv, double wdv, bool dd, double wdd)
{
    // construct total kernel
    Array1d<const sparseMat*> As;
    Array1d<SparseMatAux> matspec;
    As.push_back(&A);
    matspec.push_back(SparseMatAux(1.0,&tmp_data,&tmp_node));
    int ndata_plus=ndata_valid;
    if (gravity){
	As.push_back(&B);
	matspec.push_back(SparseMatAux(1.0,&tmp_gravdata,&tmp_node));
	ndata_plus += ngravdata;
    }
    if (sv){
	if (wsv<0) error("TomographicInversion2d::_solve - negative wsv");
	As.push_back(&Rv_h);
	matspec.push_back(SparseMatAux(wsv,&tmp_nodev,&tmp_nodev));
	ndata_plus += nnodev;
	As.push_back(&Rv_v);
	matspec.push_back(SparseMatAux(wsv,&tmp_nodev,&tmp_nodev));
	ndata_plus += nnodev;
    }
    if (nrefl>0 && sd){
	if (wsd<0) error("TomographicInversion2d::_solve - negative wsd");
	As.push_back(&Rd);
	matspec.push_back(SparseMatAux(wsd,&tmp_nodedr,&tmp_nodedc));
	ndata_plus += nnoded;
    }
    if (dv){
	if (wdv<0) error("TomographicInversion2d::_solve - negative wdv");
	As.push_back(&Tv);
	matspec.push_back(SparseMatAux(wdv,&tmp_nodev,&tmp_nodev));
	ndata_plus += nnodev;
    }
    if (nrefl>0 && dd){
	if (wdd<0) error("TomographicInversion2d::_solve - negative wdd");
	As.push_back(&Td);
	matspec.push_back(SparseMatAux(wdd,&tmp_nodedr,&tmp_nodedc));
	ndata_plus += nnoded;
    }
    SparseRectangular B(As,matspec,nnode_total);

    // construct total data vector
    total_data_vec.resize(ndata_plus);
    int idata=1;
    for (int i=1; i<=ndata; i++){
	if (tmp_data(i)>0){
	    total_data_vec(idata++) = data_vec(i);
	}
    }
    if (gravity){
	for (int i=1; i<=ngravdata; i++){
	    total_data_vec(idata++) = res_grav(i);
	}
    }
    if (sv){
	if (jumping){
	    Array1d<double> Rvm(nnodev), dmv_vec(nnodev);
	    SparseRectangular sparseRv_h(Rv_h,tmp_nodev,tmp_nodev);
	    for (int i=1; i<=nnodev; i++) dmv_vec(i) = dmodel_total_sum(i);
	    sparseRv_h.Ax(dmv_vec,Rvm);
	    for (int i=1; i<=nnodev; i++) total_data_vec(idata++) = -wsv*Rvm(i); // Lh

	    SparseRectangular sparseRv_v(Rv_v,tmp_nodev,tmp_nodev);
	    sparseRv_v.Ax(dmv_vec,Rvm);
	    for (int i=1; i<=nnodev; i++) total_data_vec(idata++) = -wsv*Rvm(i); // Lv
	}else{
	    for (int i=1; i<=nnodev; i++) total_data_vec(idata++) = 0.0; // Lh
	    for (int i=1; i<=nnodev; i++) total_data_vec(idata++) = 0.0; // Lv
	}
    }
    if (nrefl>0 && sd){
	if (jumping){
	    Array1d<double> Rdm(nnoded), dmd_vec(nnoded);
	    SparseRectangular sparseRd(Rd,tmp_nodedr,tmp_nodedr);
	    for (int i=1+nnodev, j=1; i<=nnoded+nnodev; i++) dmd_vec(j++) = dmodel_total_sum(i);
	    sparseRd.Ax(dmd_vec,Rdm);
	    for (int i=1; i<=nnoded; i++) total_data_vec(idata++) = -wsd*Rdm(i); // Ld
	}else{
	    for (int i=1; i<=nnoded; i++) total_data_vec(idata++) = 0.0;
	}
    }
    if (dv){
	if (jumping){
	    Array1d<double> Tvm(nnodev), dmv_vec(nnodev);
	    SparseRectangular sparseTv(Tv,tmp_nodev,tmp_nodev);
	    for (int i=1; i<=nnodev; i++) dmv_vec(i) = dmodel_total_sum(i);
	    sparseTv.Ax(dmv_vec,Tvm);
	    for (int i=1; i<=nnodev; i++) total_data_vec(idata++) = -wdv*Tvm(i); 
	}else{
	    for (int i=1; i<=nnodev; i++) total_data_vec(idata++) = 0.0;
	}
    }
    if (nrefl>0 && dd){
	if (jumping){
	    Array1d<double> Tdm(nnoded), dmd_vec(nnoded);
	    SparseRectangular sparseTd(Td,tmp_nodedr,tmp_nodedr);
	    for (int i=1+nnodev, j=1; i<=nnoded+nnodev; i++) dmd_vec(j++) = dmodel_total_sum(i);
	    sparseTd.Ax(dmd_vec,Tdm);
	    for (int i=1; i<=nnoded; i++) total_data_vec(idata++) = -wdd*Tdm(i); 
	}else{
	    for (int i=1; i<=nnoded; i++) total_data_vec(idata++) = 0.0;
	}
    }

    int lsqr_itermax=itermax_LSQR;
    double chi;

    int istop=iterativeSolver_LSQR(B,total_data_vec,dmodel_total,LSQR_ATOL,lsqr_itermax,chi);

    return lsqr_itermax;
}

void TomographicInversion2d::fixed_damping(int& iter, int& n, double& wdv, double& wdd)
{
    wdv = weight_d_v;
    wdd = weight_d_d;
    
    iter += _solve(smooth_velocity,weight_s_v,smooth_depth,weight_s_d,
		   true,weight_d_v,true,weight_d_d);
    n++;

    if (robust){
	if (verbose_level>=0){
	    cerr << "TomographicInversion2d:: check outliers... ";
	}
	int ndata_valid_orig=ndata_valid;
	removeOutliers();
	if (verbose_level>=0){
	    cerr << ndata_valid_orig-ndata_valid << " found\n";
	}
	if (ndata_valid < ndata_valid_orig){
	    if (verbose_level>=0){
		cerr << "TomographicInversion2d:: re-inverting without outliers...\n";
	    }
	    iter += _solve(smooth_velocity,weight_s_v,smooth_depth,weight_s_d,
			   true,weight_d_v,true,weight_d_d);
	    n++;
	}
    }
}

void TomographicInversion2d::auto_damping(int& iter, int& n, double& wdv, double& wdd)
{
    // check if damping is necessary
    iter += _solve(smooth_velocity,weight_s_v,smooth_depth,weight_s_d,
		   false,0.0,false,0.0);
    n++;

    if (robust){
	if (verbose_level>=0){
	    cerr << "TomographicInversion2d:: check outliers... ";
	}
	int ndata_valid_orig=ndata_valid;
	removeOutliers();
	if (verbose_level>=0){
	    cerr << ndata_valid_orig-ndata_valid << " found\n";
	}
	if (ndata_valid < ndata_valid_orig){
	    if (verbose_level>=0){
		cerr << "TomographicInversion2d:: re-inverting without outliers...\n";
	    }
	    iter += _solve(smooth_velocity,weight_s_v,smooth_depth,weight_s_d,
			   false,0.0,false,0.0);
	    n++;
	}
    }

    double ave_dmv0 = calc_ave_dmv();
    double ave_dmd0 = calc_ave_dmd();
    if (verbose_level>=0){
	cerr << "\t\tave_dm = " << ave_dmv0*100
	     << "%, " << ave_dmd0*100 << "% with no damping\n";
    }
    wdv = wdd = 1.0;
    if (ave_dmv0<=target_dv || !damp_velocity) wdv = 0.0;
    if (ave_dmd0<=target_dd || !damp_depth) wdd = 0.0; 
    if (wdv==0.0 && wdd==0.0) return;

    if (wdd>0.0) wdd = auto_damping_depth(wdv,iter,n);
    if (wdv>0.0) wdv = auto_damping_vel(wdd,iter,n);
}

double TomographicInversion2d::auto_damping_depth(double wdv, int& iter, int& n)
{
    // secant and bisection search
    if (verbose_level>=0) cerr << "\t\tsearching weight_d_depth...\n";
    const double wd_max = 1e7; // absolute bound
    double wd1 = 1.0; // initial guess
    double wd2 = 1e2;
    const double ddm = target_dd*0.3; // target accuracy = 30 %
    const int secant_itermax=10;
    double rts, xl, swap, dx, fl, f;

    iter += _solve(smooth_velocity,weight_s_v,smooth_depth,weight_s_d,
		   damp_velocity,wdv,true,wd1); n++;
    fl = calc_ave_dmd()-target_dd;
    if (verbose_level>=0){
	cerr << "\t\tave_dm = " << (fl+target_dd)*100 << "% at wd1("
	     << wd1 << ")\n";
    }
    iter += _solve(smooth_velocity,weight_s_v,smooth_depth,weight_s_d,
		   damp_velocity,wdv,true,wd2); n++;
    f = calc_ave_dmd()-target_dd;
    if (verbose_level>=0){
	cerr << "\t\tave_dm = " << (f+target_dd)*100 << "% at wd2("
	     << wd2 << ")\n";
    }

    if (abs(fl) < abs(f)){
	rts = wd1;
	xl = wd2;
	swap=fl; fl=f; f=swap;
    }else{
	xl = wd1; rts = wd2;
    }
    for (int j=1; j<=secant_itermax; j++){
	dx = (xl-rts)*f/(f-fl);
	if (rts+dx<=0){ // switch to bisection
	    xl=rts; fl=f;
	    rts *= 0.5;
	}else if(rts+dx>wd_max){
	    xl=rts; fl=f;
	    rts += (wd_max-rts)*0.5;
	}else{
	    xl=rts; fl=f;
	    rts+=dx;
	}
	iter += _solve(smooth_velocity,weight_s_v,smooth_depth,weight_s_d,
		       damp_velocity,wdv,true,rts); n++;
	f = calc_ave_dmd()-target_dd;
	if (verbose_level>=0){
	    cerr << "\t\tave_dm = " << (f+target_dd)*100
		 << "% at w_d of " << rts << '\n';
	}
	if (rts > wd_max) break;
	if (abs(f) < ddm) break;
    }
    return rts;
}

double TomographicInversion2d::auto_damping_vel(double wdd, int& iter, int& n)
{
    // secant and bisection search
    if (verbose_level>=0) cerr << "\t\tsearching weight_d_vel...\n";
    const double wd_max = 1e5; // absolute bound
    double wd1 = 1.0; // initial guess
    double wd2 = 1e2;
    const double ddm = target_dv*0.3; // target accuracy = 30 %
    const int secant_itermax=10;
    double rts, xl, swap, dx, fl, f;

    iter += _solve(smooth_velocity,weight_s_v,smooth_depth,weight_s_d,
		   true,wd1,damp_depth,wdd); n++;
    fl = calc_ave_dmv()-target_dv;
    if (verbose_level>=0){
	cerr << "\t\tave_dm = " << (fl+target_dv)*100 << "% at wd1("
	     << wd1 << ")\n";
    }
    iter += _solve(smooth_velocity,weight_s_v,smooth_depth,weight_s_d,
		   true,wd2,damp_depth,wdd); n++;
    f = calc_ave_dmv()-target_dv;
    if (verbose_level>=0){
	cerr << "\t\tave_dm = " << (f+target_dv)*100 << "% at wd2("
	     << wd2 << ")\n";
    }

    if (abs(fl) < abs(f)){
	rts = wd1;
	xl = wd2;
	swap=fl; fl=f; f=swap;
    }else{
	xl = wd1; rts = wd2;
    }
    for (int j=1; j<=secant_itermax; j++){
	dx = (xl-rts)*f/(f-fl);
	if (rts+dx<=0){ // switch to bisection
	    xl=rts; fl=f;
	    rts *= 0.5;
	}else{
	    xl=rts; fl=f;
	    rts+=dx;
	}
	iter += _solve(smooth_velocity,weight_s_v,smooth_depth,weight_s_d,
		       true,rts,damp_depth,wdd); n++;
	f = calc_ave_dmv()-target_dv;
	if (verbose_level>=0){
	    cerr << "\t\tave_dm = " << (f+target_dv)*100
		 << "% at w_d of " << rts << '\n';
	}
	if (rts > wd_max) break;
	if (abs(f) < ddm) break;
    }
    return rts;
}