seisinversionproceeds.cpp

石油勘探专业算法

	_t0s_num = (_et0 - _st0)/_dt0 + 1;
	_t0s = new float[_t0s_num];
	for(int m = 0;m<_t0s_num;m++) {
		_t0s[m] = (float)(_st0 + _dt0 * m);
	}
	for(i = 0;i<_curves_in_survey_num;i++) { 
		curves[i]->OpenScatterCurve();
		_curves_in_survey[i] = new wDBInversionCurveInSurvey(curves[i],_st0,_et0,_dt0);
		delete curves[i];
	}
	delete[] curves;
	delete scatter_data;
	return 0;
}

float SeisInversion::Get2DDeltaT0(int idx,int num,wURHorizonRestrictPars* par,
		int& cdp_max_pos,float& t01_max,float& t02_max)
{
	float delta_t0 = 0.0;
	int val_num = 0;
	float* temp1 = NULL;
	if(idx == 0) {
		t01_max = par->_t0_depth;
		temp1 = _2d_hor_curves[0].ObtainValues(val_num);
		for(int k = 0;k<_2d_cdp_num;k++) {
			if( (temp1[k] - t01_max) > delta_t0) {
				cdp_max_pos = _scdp + _cdp_inv * k;
				delta_t0 = temp1[k] - t01_max;
				t02_max = temp1[k];
			}
		}
		delta_t0 = delta_t0 / ((float)par->_layer_num + 1);
		return delta_t0;
	}
	if(idx == num - 1) {
		t02_max = par->_t0_depth;
		temp1 = _2d_hor_curves[num - 2].ObtainValues(val_num);
		for(int k = 0;k<_2d_cdp_num;k++) {
			if( (t02_max - temp1[k]) > delta_t0) {
				cdp_max_pos = _scdp + _cdp_inv * k;
				delta_t0 = t02_max - temp1[k];
				t01_max = temp1[k];
			}
		}
		delta_t0 = delta_t0 / ((float)par->_layer_num + 1);
		return delta_t0;
	}	

	int i = 0,j = 0;
	temp1 = _2d_hor_curves[idx - 1].ObtainValues(val_num);
	float* temp2 = _2d_hor_curves[idx].ObtainValues(val_num);
	
	for(i = 0;i<_2d_cdp_num;i++) {
		if( ( temp1[i] > 0)&&(temp2[i] > 0)
			&&( (temp2[i] - temp1[i])> delta_t0) ){
			cdp_max_pos = _scdp + i * _cdp_inv;
			t01_max = temp1[i];
			t02_max = temp2[i];
			delta_t0 = t02_max - t01_max;
		}
	}
	delta_t0 = delta_t0 / ((float)par->_layer_num + 1);
	return delta_t0;
}

void SeisInversion::Get2DCDPForInterpolate(int& cdp,float& t0,wURHorizonRestrictPars* par,
	float delta_t0,int cdp_max_pos,float t01_max,float t02_max,int i,int j,int k)
{
	cdp = _scdp + _cdp_inv * k;
	float *temp = NULL;
	int num = 0;

	float hor_max_pos1,hor_max_pos2;
	float hor_cur_pos1,hor_cur_pos2;

	if(par->_top_hor_id[0] != '\0') {
		temp = _2d_hor_curves[i - 1].ObtainValues(num);
		hor_cur_pos1 = temp[k];
	} else {
		hor_cur_pos1 = par->_t0_depth;
	}
	
	if(par->_bottom_hor_id[0] != '\0') {
		temp = _2d_hor_curves[i].ObtainValues(num);
		hor_cur_pos2 = temp[k];
	} else {
		hor_cur_pos2 = par->_t0_depth;
	}

	switch(par->_type) {
		case NonRestrictLayer: {
			t0 = t01_max + (j + 1) * delta_t0;
			break;
		}
		case ParallelTopHorizon: {
			t0 = hor_cur_pos1 + (j + 1) * delta_t0;
			break;
		}
		case ParallelBottomHorizon: {
			t0 = hor_cur_pos2 - (par->_layer_num - j) * delta_t0;
			break;
		}
		case IsohedralBetweenBorders: {
			t0 = hor_cur_pos1 + (j + 1) * delta_t0 / (t02_max - t01_max ) * (hor_cur_pos2 - hor_cur_pos1);
			break;
		}
	}
	if((hor_cur_pos1 > 0)&&(t0 < hor_cur_pos1)) t0 = -1.0;
	if((hor_cur_pos2 > 0)&&(t0 > hor_cur_pos2)) t0 = -1.0;
}

void SeisInversion::Get2DCurveT0(float& hor_cur_pos1,float& hor_cur_pos2,int cdp,wURHorizonRestrictPars* par,int i)
{
	int k = 0;
	k = (cdp - _scdp) / _cdp_inv;
	float *temp = NULL;
	int num = 0;

	if(par->_top_hor_id[0] != '\0') {
		temp = _2d_hor_curves[i - 1].ObtainValues(num);
		hor_cur_pos1 = temp[k];
	} else {
		hor_cur_pos1 = par->_t0_depth;
	}
	
	if(par->_bottom_hor_id[0] != '\0') {
		temp = _2d_hor_curves[i].ObtainValues(num);
		hor_cur_pos2 = temp[k];
	} else {
		hor_cur_pos2 = par->_t0_depth;
	}
	return;
}

wDBInversionCurveInSurvey** SeisInversion::Get2DInterpolateCurvesForLayer(int layer_idx,
		wURHorizonRestrictPars* par,int& num)
{
	int i = 0;
	wDBInversionCurveInSurvey** result = NULL;
	float hor_cur_pos1,hor_cur_pos2;
	int line_id1;
	int cdp1;

	result = new wDBInversionCurveInSurvey*[_curves_in_survey_num];
	for(i = 0;i<_curves_in_survey_num;i++) {
		result[i] = new wDBInversionCurveInSurvey(_curves_in_survey[i]);
		result[i]->GetLineAndCDP(line_id1,cdp1);
		float* values = NULL;
		float* t0s = NULL;
		int values_num = 0;
		values_num = result[i]->GetPoints(t0s,values);
		delete[] t0s;
		Get2DCurveT0(hor_cur_pos1,hor_cur_pos2,cdp1,par,layer_idx);
		for(int n = 0;n<values_num;n++) {
			if( (hor_cur_pos1 > 0)&&(_t0s[n] < hor_cur_pos1 - (float)_dt0) ) values[n] = -1.0;
			if( (hor_cur_pos2 > 0)&&(_t0s[n] > hor_cur_pos2 + (float)_dt0) ) values[n] = -1.0;
		}
		result[i]->SetPoints(_t0s,values,values_num);
		delete[] values;
	}
	num = _curves_in_survey_num;
	return result;
}

int SeisInversion::Get2DInterpolateValue(wDBInversionCurveInSurvey** curves,int num,int zone_idx,
	wURHorizonRestrictPars* par,int layer_idx,
	float delta_t0,int cdp_max_pos,
	float t01_max,float t02_max,
	float*& interpolate_values,int*& interpolate_cdps)
{
	float t0 = 0.0;
	int i = 0,j = 0;
	float hor_cur_pos1,hor_cur_pos2;
	int line_id1;
	int cdp1;
	float* values = new float[num];
	int * cdps = new int[num];
	float temp = 0.0;
	
	for(i = 0;i<num;i++) {
		curves[i]->GetLineAndCDP(line_id1,cdp1);
		Get2DCurveT0(hor_cur_pos1,hor_cur_pos2,cdp1,par,zone_idx);
		switch(par->_type) {
			case NonRestrictLayer: {
				t0 = t01_max + (layer_idx + 1) * delta_t0;
				break;
			}
			case ParallelTopHorizon: {
				t0 = hor_cur_pos1 + (layer_idx + 1) * delta_t0;
				break;
			}
			case ParallelBottomHorizon: {
				t0 = hor_cur_pos2 - (par->_layer_num - layer_idx) * delta_t0;
				break;
			}
			case IsohedralBetweenBorders: {
				t0 = hor_cur_pos1 + (layer_idx + 1) * delta_t0 / (t02_max - t01_max ) *
					(hor_cur_pos2 - hor_cur_pos1);
				break;
			}
		}
		temp = curves[i]->GetValue(t0);
		if(temp > 0){
			values[j] = temp;
			cdps[j] = cdp1;
			j++;
		}
	}
	if(j > 0){
		interpolate_values = new float[j];
		interpolate_cdps = new int[j];
		memcpy(interpolate_values,values,j * sizeof(float));
		memcpy(interpolate_cdps,cdps,j * sizeof(int));
	}
	delete[] values;
	delete[] cdps;
	return j;
}

float SeisInversion::Get2DWeighValue(int cdp,float* interpolate_values,int* interpolate_cdps,
		int interpolate_values_num)
{
	float value = 0.0;
	int dd = 0;
	float val = -1.0;
	float temp1 = 0.0,temp2 = 0.0;

	if(interpolate_values_num == 1) return interpolate_values[0];
	for(int i = 0;i < interpolate_values_num;i++) {
		dd = abs(interpolate_cdps[i] - cdp);
		if(dd == 0) {
			return interpolate_values[i];
		}
		temp1 = temp1 + interpolate_values[i] / (float)dd;
		temp2 = temp2 + 1.0 / (float)dd;
	}
	if( (temp1 > 0)&&(temp2 > 0) )
	{
		val = temp1 / temp2;
	}
	return val;
}

void SeisInversion::Interpolate2DOperation()
{
	time_t clock;

	time(&clock);
	char*time_str = NULL;
	time_str = ctime(&clock);
	cout<<"Interpolate Operation Begin Time: "<<time_str<<endl;

	int rt_code = Get2DHorizonSurfaceData();
	if(rt_code != 0){
		cout<<"Get2DHorizonSurfaceData() ERROR!"<<endl;
		cout<<"Open wURSeisHorizonCurve Data ERROR"<<endl;
		return;
	}

	int i = 0,j = 0,k = 0,m = 0;
	float delta_t0 = 0.0;
	wURHorizonRestrictPars* hor_sestrict_par = NULL;
	int num = 0;
	hor_sestrict_par = _seis_tpt->GetHorizonPars(num);

	SeisInversionProgressDialog* dlg = NULL;
	dlg = new SeisInversionProgressDialog(this,"SeisInversion2DProgressDialog");
	dlg->Manage();

	int cdp_max_pos = 0;
	float t01_max = 0.0;
	float t02_max = 0.0;
	float pos1,pos2,pos3,pos4;

	int buffer_idx = 0;
	for(i = 0;i<num;i++) {
		delta_t0 = Get2DDeltaT0(i,num,&hor_sestrict_par[i],cdp_max_pos,t01_max,t02_max);
		pos1 = 100.0 * (float)(i)/(float)num;
		dlg->SetCurrentValues((int)pos1,-1,-1,-1);
		wDBInversionCurveInSurvey** curves = NULL;
		int curves_num = 0;
		curves = Get2DInterpolateCurvesForLayer(i,
			&hor_sestrict_par[i],curves_num);
		for(j = 0;j<hor_sestrict_par[i]._layer_num + 1;j++) {
			pos2 = 100.0 * (float)(j)/(float)(hor_sestrict_par[i]._layer_num + 1);
			dlg->SetCurrentValues((int)pos1,(int)pos2,-1,-1);
			float *interpolate_values = NULL;
			int *interpolate_cdps = NULL;
			int interpolate_values_num = 0;
			interpolate_values_num = Get2DInterpolateValue(curves,curves_num,i,
				&hor_sestrict_par[i],j,
				delta_t0,cdp_max_pos,t01_max,t02_max,
				interpolate_values,interpolate_cdps);
			for(k = 0;k<_2d_cdp_num;k++) {
				pos3 = 100.0 * (float)(k + 1)/(float)_2d_cdp_num;
				dlg->SetCurrentValues((int)pos1,(int)pos2,(int)pos3,-1);
				int line_id = 0,cdp = 0;
				float t0 = 0.0;
				float value = 0.0;
				Get2DCDPForInterpolate(cdp,t0,&hor_sestrict_par[i],
					delta_t0,cdp_max_pos,t01_max,t02_max,i,j,k);
				if(t0 > 0) {
					value = Get2DWeighValue(cdp,interpolate_values,interpolate_cdps,
						interpolate_values_num);
				} else {
					value = -1.0;
				}
				if( (k > 0)&&(t0 > 0)&&(value < 0) ){
					_2d_buffers_t0[buffer_idx][k] = _2d_buffers_t0[buffer_idx][k - 1];
					_2d_buffers_v0[buffer_idx][k] = _2d_buffers_v0[buffer_idx][k - 1];
				} else {
					_2d_buffers_t0[buffer_idx][k] = t0;
					_2d_buffers_v0[buffer_idx][k] = value;
				}
			}//for(k = 0;k<_hor_surface_data2_num1;k++){
			buffer_idx++;
			delete[] interpolate_values;
			delete[] interpolate_cdps;
		}//for(j = 0;j<hor_sestrict_par[i]._layer_num;j++) {
		for(int mmm = 0;mmm < curves_num;mmm++) delete curves[mmm];
		delete[] curves;
	}//for(i = 0;i<num;i++) {

	time(&clock);
	time_str = ctime(&clock);
	cout<<"Interpolate Operation End Time: "<<time_str<<endl;

	int status = 0;
	if(_seis_attr_pro) delete _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);
	wDBSeis2DSection* seis2d_section_data = NULL;
	seis2d_section_data = _seis_attr_pro->Create2DSection("Lfv");
	float* inter_data_t0 = new float[_2d_buffers_num];
	float* inter_data_v0 = new float[_2d_buffers_num];
	_t0s_num = (_et0 - _st0)/_dt0 + 1;
	float* section_data = new float[_t0s_num * _2d_cdp_num];

	for(i = 0;i<_2d_cdp_num;i++) {
		pos4 = 100.0 * (float)(i + 1)/(float)_2d_cdp_num;
		dlg->SetCurrentValues(100,100,100,(int)pos4);
		int count = 0;
		for(j = 0;j<_2d_buffers_num;j++) {
			if( (_2d_buffers_t0[j][i] > 0)&&(_2d_buffers_v0[j][i] > 0) ){
				inter_data_t0[count] = _2d_buffers_t0[j][i];
				inter_data_v0[count] = _2d_buffers_v0[j][i];
				if( (count > 0)&&(inter_data_t0[count] > inter_data_t0[count - 1]) ) count++;
				if(count == 0) count++;
			}
		}
		if(count == 0) {
//			if(i > 0) memcpy(section_data + i * _t0s_num,
//				section_data + (i - 1) * _t0s_num,_t0s_num * sizeof(float));
			for(int ppk = 0;ppk<_t0s_num;ppk++){
				*(section_data + i * _t0s_num + ppk) = 0.0;
			}
			continue;
		}
		float* result = NULL;
		int result_num = 0;
		result = GetLinearInterpolation(inter_data_t0,inter_data_v0,count,
			(float)_st0,(float)_et0,(float)_dt0,result_num);
		memcpy(section_data + i * _t0s_num,result,result_num * sizeof(float));
		delete result;
	}
	status = seis2d_section_data->WriteSection(section_data);
	delete section_data;
	delete inter_data_t0;
	delete inter_data_v0;
	delete seis2d_section_data;
	delete dlg;
	time(&clock);
	time_str = ctime(&clock);
	cout<<"End Time: "<<time_str<<endl;
}