seisinversionproceeds1.cpp

石油勘探专业算法：地震资料波阻抗反演后期算法

		memcpy(bufx,bufo,_t0s_num * sizeof(float));
		result = NULL;
		result = SeisDataFilter(bufx,kpoint,_dt0,_freq1,_freq2,0.0,0.0);
		for(int mmk = 0;mmk<_t0s_num;mmk++) bufo[mmk] = bufx[mmk] - result[mmk];
		delete result;

		for(i2 = 0;i2<kpoint;i2++){
			lfvdata[i2] = lfvdata[i2] + _Vlog_vmark * bufo[i2];
		}
		if(_Vlog_flag == 1){
			memcpy(rimp2d_data + temp * _t0s_num,bufo,_t0s_num * sizeof(float));
			memcpy(rimp2d_data1 + temp * _t0s_num,bufx,_t0s_num * sizeof(float));
		}
		memcpy(imp2d_data + temp * _t0s_num,lfvdata,_t0s_num * sizeof(float));
	}
	if(_Vlog_flag == 1) {
		rimp2d_section_data->WriteSection(rimp2d_data);
		rimp2d_section_data1->WriteSection(rimp2d_data1);
		rimp2d_section_data2->WriteSection(rimp2d_data2);
	}
	imp2d_section_data->WriteSection(imp2d_data);
	if(_Vlog_flag == 1){
		delete rimp2d_section_data2;
		delete rimp2d_section_data1;
		delete rimp2d_section_data;
	}
	delete rimp2d_data;
	delete rimp2d_data1;
	delete rimp2d_data2;
	delete imp2d_data;
	delete imp2d_section_data;
	delete lfv2d_section_data;
	delete seis2d_section_data;

	delete kdimdata;
	delete bufx;
	delete bufo;
	delete lfvdata;
	delete dlg;
    return(1);
}


int SeisInversion::V_logBatchProcess2()
{
	char* dname = NULL;
	dname = _seis_tpt->Obtain3DTemplate()->GetDataName();
	if( (!dname)||(strcmp(dname,"Amp") != 0) ) {
		cout<<"Read Data ERROR!"<<endl;
		return 0;
	}
	delete dname;

	int i2,i3,ist,iet,il=61;
	int kcdpi , knumpoint , kpoint,kline,status;
	float ymax;
	float *kdimdata ,*proceeddimdata , *bufx , *bufo , *lfvdata , *zerodata;
	int init_status;
	kcdpi = (_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1;
	knumpoint = (_et0 - _st0) / _dt0 + 1;
	_t0s_num = knumpoint;
	kpoint = knumpoint;
	kcdpi = knumpoint;
	lfvdata = new float[kcdpi];

        if(!_seis_attr_pro) {
	        wDBSeis3DVolumeRange* seis_3d_range = NULL;
	        seis_3d_range = _seis_tpt->Obtain3DTemplate()->ObtainSeis3D();
        	_seis_attr_pro = new wPRSeisAttrsBase();
	        _seis_attr_pro->SetOrigDataName("Amp");
	        _seis_attr_pro->SetSurveyID(_survey_id);
	        _seis_attr_pro->SetDomain(_data_domain);
	        _seis_attr_pro->SetResultDataName("Lfv");
	        _seis_attr_pro->SetProcessType(Seis3DVolumeProcessing);
	        _seis_attr_pro->Set3DRange(seis_3d_range);
	        _seis_attr_pro->SetCreateFlag(SeisLineVolume);
	}
        wURSeis3D* seis3d_volume_data = NULL;
	_seis_attr_pro->SetOrigDataName("Amp");
        seis3d_volume_data = _seis_attr_pro->GetCurrentVolume("Amp");

	_seis_attr_pro->SetOrigDataName("Lfv");
	wURSeis3D* lfv3d_volume_data = NULL;
	lfv3d_volume_data = _seis_attr_pro->GetCurrentVolume("Lfv");

	wDBSeis3DVolumeData* imp3d_volume_data = NULL;
	imp3d_volume_data = _seis_attr_pro->Create3DVolume("Imp");

	wDBSeis3DVolumeData* ting3d_volume_data = NULL;
	if(_Vlog_flag2 == 1)
		ting3d_volume_data = _seis_attr_pro->Create3DVolume("Ting");

	float* seis3d_data = NULL;
	int seis3d_data_num1,seis3d_data_num2;
	float* lfv3d_data = NULL;
	int lfv_data_num1,lfv_data_num2;
	float* imp3d_data = new float[((_ecdp_proceed - _scdp_proceed)/_lid_inv + 1) * _t0s_num];
	float* ting3d_data = new float[((_ecdp_proceed - _scdp_proceed)/_lid_inv + 1) * _t0s_num];

	int kdimdata_num = 1;
	while(1) {
		kdimdata_num = 2 * kdimdata_num;
		if(_t0s_num <= kdimdata_num) break;
	}
	kdimdata_num = 2 * kdimdata_num;
	kdimdata = new float[kdimdata_num];
	proceeddimdata = new float[kdimdata_num];
	bufx = new float[kdimdata_num];
	bufo = new float[kdimdata_num];
	zerodata = new float[kdimdata_num];
	for(int p = 0;p<kdimdata_num;p++){
		kdimdata[p] = 0.0;
		proceeddimdata[p] = 0.0;
		bufx[p] = 0.0;
		bufo[p] = 0.0;
		zerodata[p] = 0.0;
	}


	SeisInversionProgressDialog1* dlg = NULL;
	dlg = new SeisInversionProgressDialog1(this);
	dlg->Manage();
	float pos1,pos2;
	for(kline = _slid_proceed;kline <= _elid_proceed;kline += _lid_inv){
		pos1 = 100.0 * (float)(kline - _slid_proceed)/(float)(_elid_proceed - _slid_proceed);
		dlg->SetCurrentValues((int)pos1,-1);

		seis3d_volume_data->SetCurrentLine(kline);
		seis3d_data = seis3d_volume_data->ObtainValues(seis3d_data_num1,seis3d_data_num2);
		lfv3d_volume_data->SetCurrentLine(kline);
		lfv3d_data = lfv3d_volume_data->ObtainValues(lfv_data_num1,lfv_data_num2);

	        for( kcdpi = _scdp_proceed;kcdpi <= _ecdp_proceed;kcdpi += _cdp_inv) {
			pos2 = 100.0 * (float)(kcdpi - _scdp_proceed)/(float)(_ecdp_proceed - _scdp_proceed);
			dlg->SetCurrentValues((int)pos1,(int)pos2);

			memcpy(kdimdata,zerodata,kdimdata_num * sizeof(float));
			memcpy(bufo,zerodata,kdimdata_num * sizeof(float));
			memcpy(bufx,zerodata,kdimdata_num * sizeof(float));
			int temp = (kcdpi - _scdp_proceed)/_cdp_inv;
			memcpy(kdimdata,seis3d_data + temp * _t0s_num,_t0s_num * sizeof(float));
			memcpy(lfvdata,lfv3d_data + temp * _t0s_num,_t0s_num * sizeof(float));
			float sum = 0.0;
			int tempi = 0;
			for(tempi = 0;tempi<_t0s_num;tempi++) {
				if(kdimdata[tempi] > 0)
				sum = kdimdata[tempi];
			}
			if(sum <= 0.0) {
				memcpy(imp3d_data + temp * _t0s_num,lfvdata,_t0s_num * sizeof(float));
				continue;
			}
			SeismicTraceIntegration(_t0s_num,(float)_dt0 / 1000.0,_min_freq_value,
					kdimdata_num,kdimdata,proceeddimdata,bufx,bufo);
			for(int tempj = 0;tempj<_t0s_num;tempj++){
				float factor = exp(_Vlog_vmark2 * proceeddimdata[tempj]);
				if(factor < 0.01) factor = 0.01;
				if(factor > 1.5) factor = 1.5;
				lfvdata[tempj] = lfvdata[tempj] * factor;
			}
			memcpy(ting3d_data + temp * _t0s_num,proceeddimdata,_t0s_num * sizeof(float));
			memcpy(imp3d_data + temp * _t0s_num,lfvdata,_t0s_num * sizeof(float));
		}
		if(_Vlog_flag2 == 1)
			ting3d_volume_data->WriteInLine(kline,ting3d_data);
		imp3d_volume_data->WriteInLine(kline,imp3d_data);
	}
	delete ting3d_data;
	delete imp3d_data;
	
	if(_Vlog_flag2 == 1)
		delete ting3d_volume_data;
		
	delete imp3d_volume_data;
	delete lfv3d_volume_data;
	delete seis3d_volume_data;

	delete kdimdata;
	delete proceeddimdata;
	delete bufx;
	delete bufo;
	delete lfvdata;
	delete zerodata;
	delete dlg;
	return(1);
}


int SeisInversion::V_logBatch2DProcess2()
{
        float max_value;
        float min_value;

	char* dname = NULL;
	dname = _seis_tpt->Obtain2DTemplate()->GetDataName();
	if( (!dname)||(strcmp(dname,"Amp") != 0) ) {
		cout<<"Read Data ERROR!"<<endl;
		return 0;
	}
	delete dname;

	int i2,i3,ist,iet,il=61;
	int kcdpi , knumpoint , kpoint,kline,status;
	float ymax;
	float *kdimdata , *proceeddimdata , *bufx , *bufo , *lfvdata , *zerodata;
	int init_status;
	kcdpi = (_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1;
	knumpoint = (_et0 - _st0) / _dt0 + 1;
	_t0s_num = knumpoint;
	kpoint = knumpoint;
	kcdpi = knumpoint;
	lfvdata = new float[kcdpi];

        if(!_seis_attr_pro) {
		wDBSeis2DSectionRange* seis_2d_range = NULL;
		seis_2d_range = _seis_tpt->Obtain2DTemplate()->ObtainSeis2D();
		_seis_attr_pro = new wPRSeisAttrsBase();
		_seis_attr_pro->SetOrigDataName("Amp");
		_seis_attr_pro->SetSurveyID(_survey_id);
		_seis_attr_pro->SetDomain(_data_domain);
		_line_id = _seis_tpt->Obtain2DTemplate()->GetLineID();
		_seis_attr_pro->SetLineID(_line_id);
		_seis_attr_pro->SetResultDataName("Lfv");
		_seis_attr_pro->SetProcessType(Seis2DSectionProcessing);
		_seis_attr_pro->Set2DRange(seis_2d_range);
		_seis_attr_pro->SetCreateFlag(SeisLineVolume);
	}
        wURSeis2D* seis2d_section_data = NULL;
	_seis_attr_pro->SetOrigDataName("Amp");
        seis2d_section_data = _seis_attr_pro->GetCurrentSection("Amp");

	wURSeis2D* lfv2d_section_data = NULL;
	_seis_attr_pro->SetOrigDataName("Lfv");
	lfv2d_section_data = _seis_attr_pro->GetCurrentSection("Lfv");

	wDBSeis2DSection* imp2d_section_data = NULL;
	imp2d_section_data = _seis_attr_pro->Create2DSection("Imp");

	wDBSeis2DSection* ting2d_section_data = NULL;
	if(_Vlog_flag2 == 1) {
		ting2d_section_data = _seis_attr_pro->Create2DSection("Ting");
	}

	float* seis2d_data = NULL;
	int seis2d_data_num1,seis2d_data_num2;
	float* lfv2d_data = NULL;
	int lfv2d_data_num1,lfv2d_data_num2;
	float* imp2d_data = new float[((_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1) * _t0s_num];
	float* ting2d_data = new float[((_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1) * _t0s_num];

	seis2d_section_data->SetCurrentLine(kline);
	seis2d_data = seis2d_section_data->ObtainValues(seis2d_data_num1,seis2d_data_num2);
	lfv2d_section_data->SetCurrentLine(kline);
	lfv2d_data = lfv2d_section_data->ObtainValues(lfv2d_data_num1,lfv2d_data_num2);

	int kdimdata_num = 1;
	while(1) {
		kdimdata_num = 2 * kdimdata_num;
		if(_t0s_num <= kdimdata_num) break;
	}
	kdimdata_num = 2 * kdimdata_num;
	kdimdata = new float[kdimdata_num];
	for(int p = 0;p<kdimdata_num;p++){
		kdimdata[p] = 0.0;
	}
	proceeddimdata = new float[kdimdata_num];
	bufx = new float[kdimdata_num];
	bufo = new float[kdimdata_num];
	zerodata = new float[kdimdata_num];
	for(int p = 0;p<kdimdata_num;p++){
		kdimdata[p] = 0.0;
		proceeddimdata[p] = 0.0;
		bufx[p] = 0.0;
		bufo[p] = 0.0;
		zerodata[p] = 0.0;
	}
	
	SeisInversionProgressDialog2* dlg = NULL;
	dlg = new SeisInversionProgressDialog2(this);
	dlg->Manage();
	float pos1;

        for( kcdpi = _scdp_proceed;kcdpi <= _ecdp_proceed;kcdpi += _cdp_inv) {
		pos1 = 100.0 * (float)(kcdpi - _scdp_proceed)
			/(float)(_ecdp_proceed - _scdp_proceed);
		dlg->SetCurrentValues((int)pos1);

		memcpy(kdimdata,zerodata,kdimdata_num * sizeof(float));
		memcpy(bufo,zerodata,kdimdata_num * sizeof(float));
		memcpy(bufx,zerodata,kdimdata_num * sizeof(float));
		int temp = (kcdpi - _scdp_proceed)/_cdp_inv;
		memcpy(kdimdata,seis2d_data + temp * _t0s_num,_t0s_num * sizeof(float));
		memcpy(lfvdata,lfv2d_data + temp * _t0s_num,_t0s_num * sizeof(float));
		float sum = 0.0;
		int tempi = 0;
		for(tempi = 0;tempi<_t0s_num;tempi++) {
			if(kdimdata[tempi] > 0)
			sum = kdimdata[tempi];
		}
		if(sum <= 0.0) {
			memcpy(imp2d_data + temp * _t0s_num,lfvdata,_t0s_num * sizeof(float));
			continue;
		}
		SeismicTraceIntegration(_t0s_num,(float)_dt0 / 1000.0,_min_freq_value,
				kdimdata_num,kdimdata,proceeddimdata,bufx,bufo);
		for(int tempj = 0;tempj<_t0s_num;tempj++){
			float factor = exp(_Vlog_vmark2 * proceeddimdata[tempj]);
			if(factor < 0.01) factor = 0.01;
			if(factor > 1.5) factor = 1.5;
			lfvdata[tempj] = lfvdata[tempj] * factor;
		}
		memcpy(ting2d_data + temp * _t0s_num,proceeddimdata,_t0s_num * sizeof(float));
		memcpy(imp2d_data + temp * _t0s_num,lfvdata,_t0s_num * sizeof(float));
	}
	if(_Vlog_flag2 == 1) {
		ting2d_section_data->WriteSection(ting2d_data);
	}
	imp2d_section_data->WriteSection(imp2d_data);
	if(_Vlog_flag2 == 1){
		delete ting2d_section_data;
	}
	delete ting2d_data;
	delete imp2d_data;
	delete ting2d_section_data;
	delete imp2d_section_data;
	delete lfv2d_section_data;
	delete seis2d_section_data;

	delete kdimdata;
	delete proceeddimdata;
	delete bufx;
	delete bufo;
	delete lfvdata;
	delete zerodata;
	delete dlg;
    return(1);
}


/************************************************************************
 * 
 * 	Seismic trace integration by use of Fourier Transformation
 *	        
 * 
 ************************************************************************/
int SeisInversion::SeismicTraceIntegration(int nsample,float deltat,float freqmin,
		/*******************************************/
		int nfft,
		float *straceo,/* original seismic traces */
		float *stracep,/* processed seismic traces */
		float* xxr,float* xxi)
		/*******************************************/
{
	int 	i/*,nfft*/;
	float 	f,ff;
	memcpy(xxr,straceo,nfft * sizeof(float));

	fft1d(nfft,1,xxr,xxi);
	float temp = 0.0;
	ff=(float)nfft*deltat*freqmin;
	if (ff<1.0) ff=1.0;
	for (i=0;i<nfft/2;i++)	{
		f=(float)i;
		if (f<ff) f=f/ff/ff;
		else	f=1.0/f;
		temp = xxr[i];
		xxr[i]=xxi[i]*f;
		xxi[i]=-temp*f;
		temp = xxr[nfft-1-i];
		xxr[nfft-1-i]=-xxi[nfft-1-i]*f;
		xxi[nfft-1-i]=temp*f;
	}

	fft1d(nfft,-1,xxr,xxi);
	memcpy(stracep,xxr,nsample * sizeof(float));
	return 0;
}

/*-----------------------------------------------------------
	        One-dimensional Fast Fourier Transform
-------------------------------------------------------------
	n----------Number of samples
	inv--------Flag:  = 1  for forward FFT
	                  =-1  for inverse FFT
	Usage:
		int n,inv;
		float xreal[],ximag[];		  
		fft1d(n,inv,xreal,ximag);
------------------------------------------------------------*/
int SeisInversion::fft1d(int n,int inv,float* xreal,float* ximag)
{	int m,nv2,nm1,i,j,k,l,le,le1,ip;	
	float treal,timag,p1,ureal,uimag,p1le1,wreal,wimag;
	float xr,xi,xrea,xima,yreal,yimag,flt;
	m=(int)(log((double)n)/log(2.0)+0.1);
	nv2=n/2;	nm1=n-1;	j=1;
	for(i=1;i<=nm1;i++)	{
		if (i<j) { 
			treal=xreal[j-1];
			timag=ximag[j-1];
			xreal[j-1]=xreal[i-1];
			ximag[j-1]=ximag[i-1];
			xreal[i-1]=treal;
 			ximag[i-1]=timag;
		}
		k=nv2;
		while (k<j) {j=j-k;	k=k/2;} 
		j=j+k;
	}
		p1=4.0*atan(1.0);
	for (l=1;l<=m;l++)	{
		le=1;
		for (i=0;i<l;i++) le*=2;
		le1=le/2;
		ureal=1.0;
		uimag=0.0;
		p1le1=p1/(float)le1;
		wreal=+cos(p1le1);
		wimag=-sin(p1le1);
		if (inv==-1) wimag=-wimag;
		for (j=1;j<=le1;j++)	{
				for (i=j;i<=n;i+=le)	{
					ip=i+le1;
					xr=xreal[ip-1];
					xi=ximag[ip-1];
					xrea=xreal[i-1];
					xima=ximag[i-1];
					treal=xr*ureal-xi*uimag;
					timag=xi*ureal+xr*uimag;
					xreal[ip-1]=xrea-treal;
					ximag[ip-1]=xima-timag;
					xreal[i-1]=+xreal[i-1]+treal;
					ximag[i-1]=+ximag[i-1]+timag;
				}	
			yreal=ureal*wreal-uimag*wimag;
			yimag=uimag*wreal+ureal*wimag;
			ureal=yreal;
			uimag=yimag;
		}
	}
	if (inv==1) return (0);
	for (i=0;i<n;i++)	{
		flt=(float)n;
		xreal[i]=xreal[i]/flt;
		ximag[i]=ximag[i]/flt;
	}
	return 0;
}