seisinversionproceeds.cpp

石油勘探专业算法

		int den_num = 0;
		float* t0s = NULL;
		well_vel_data->GetTimeRange(st0,et0,dt0);
		vel_values = well_vel_data->GetValues(vel_num);
		cout<<"vel_values[0] = "<<vel_values[0]<<endl;
		if(vel_values[0] < -9000.0) vel_values[0] = 0.0;
		if(den_idx > 0) den_values = well_den_data->GetValues(den_num);
		if(den_idx == 0){
			for(int kkk = 0;kkk<vel_num;kkk++){
				temp = vel_values[kkk];
				if(linear_flag == 1){
					temp = temp * a_val * pow(temp,b_val);
				} else {
					temp = temp * default_value;
				}
				vel_values[kkk] = temp;
			}
		} else {
			den_values = well_den_data->GetValues(den_num);
			well_den_data->GetTimeRange(st01,et01,dt01);
			int kkk = 0;
			if(st01 > st0) {
				while(st01 > st0){
					temp = vel_values[kkk];
					if(linear_flag == 1){
						temp = temp * a_val * pow(temp,b_val);
					} else {
						temp = temp * default_value;
					}
					vel_values[kkk] = temp;
					st0 = st0 + dt0;
					kkk++;
				}
			}
			if(dt0 != dt01){
				cout<<"error"<<endl;
			}
			int kkk1 = 0;
			for(kkk1 = 0;kkk1<den_num;){
				temp = vel_values[kkk];
				if(linear_flag == 1){
					temp = temp * a_val * pow(temp,b_val);
				} else {
					temp = temp * default_value;
				}
				vel_values[kkk] = temp;
				st0 = st0 + dt0;
				kkk++;
				kkk1++;
			}
			while(et01 < et0){
				temp = vel_values[kkk];
				if(linear_flag == 1){
					temp = temp * a_val * pow(temp,b_val);
				} else {
					temp = temp * default_value;
				}
				vel_values[kkk] = temp;
				st0 = st0 + dt0;
				kkk++;
			}
		}

		well_vel_data->GetTimeRange(st0,et0,dt0);
		t0s = new float[vel_num];
		for(int j = 0;j<vel_num;j++) {
			t0s[j] = (float)st0 + (float)(dt0 * j);
		}

		cout<<"freq_flag = "<<freq_flag<<endl;
		if(freq_flag) {
			values1 = SeisDataFilter(vel_values,vel_num,dt0,_freq1,_freq2,0.0,0.0);
			curve->SetPoints(t0s,values1,vel_num);
		} else {
			curve->SetPoints(t0s,vel_values,vel_num);
		}
		delete t0s;
		delete vel_values;
		delete den_values;
		delete values1;
		data->AddOneCurve(curve);
		delete curve;
		OnDrawsReflash();
	}	
}

void SeisInversion::GetHorizonSurfaceData()
{
	int i = 0,j = 0;
	wURHorizonRestrictPars* hor_sestrict_par = NULL;
	int num = 0;
	hor_sestrict_par = _seis_tpt->GetHorizonPars(num);
	for(i = 0;i < num - 1;i++) {
		if(_hor_surface_data_file_name[i]) delete _hor_surface_data_file_name[i];
		_hor_surface_data_file_name[i] = new char[50];
		sprintf(_hor_surface_data_file_name[i],"temp_hor_surface_data%d.tmp",i);
	}
	float* data = NULL;
	int num1 = 0,num2 = 0;
	for(i = 0;i < num - 1;i++) {
		FILE* fp = NULL;
		fp = fopen(_hor_surface_data_file_name[i],"wb+");
		wDBSeisHorizonInSurvey* h = NULL;
		h = new wDBSeisHorizonInSurvey();
		h->OpenHorInSurvey(hor_sestrict_par[i]._bottom_hor_id,
			hor_sestrict_par[i]._survey_id);
		char *hs_id = NULL;
		hs_id = h->GetHorSurveyID();
		delete h;
		wDBSeisHorizonSurfaceInSurvey* hs = NULL;
		hs = new wDBSeisHorizonSurfaceInSurvey();
		int status = hs->OpenHorSurface(hs_id,SliceOfTimeDepth,"TimeDepthSurf");
		if(status != OK_STATUS) { 
			delete hs;
			delete hs_id;
			continue;
		}
		data = hs->ObtainSurfaceData(num1,num2);
		if( i == 0) {
			hs->GetLineRange(_slid,_elid,_lid_inv);
			hs->GetCDPRange(_scdp,_ecdp,_cdp_inv);
			_survey_id = hor_sestrict_par[i]._survey_id;
		}
		fwrite(&num1,sizeof(int),1,fp);
		fwrite(&num2,sizeof(int),1,fp);
		fwrite(data,sizeof(float),num1 * num2,fp);
		delete hs;
		delete hs_id;
		fclose(fp);
	}
	//Get Curves:
	wDBScatterDataInSurvey *scatter_data = NULL;
	scatter_data = new wDBScatterDataInSurvey();
	wURScatterDataInSection* section_data = NULL;
	section_data = _seis_tpt->ObtainScatterSection();
	char sc_id[16];
	section_data->GetScatterID(sc_id);
	scatter_data->SetScatterID(sc_id);
	
	scatter_data->OpenScatterData();
	wDBScatterCurveInSurvey** curves = NULL;
	curves = scatter_data->ListScatterCurves(_curves_in_survey_num);
	cout<<"_curves_in_survey_num = "<<_curves_in_survey_num<<endl;
	_curves_in_survey = new wDBInversionCurveInSurvey*[_curves_in_survey_num];
	int values_num = 0;
	_seis_tpt->Obtain3DTemplate()->ObtainSeis3D()->GetTimeRange(_st0,_et0,_dt0);
	_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;
} 


float SeisInversion::GetDeltaT0AndDataFromDisk(int idx,int num,wURHorizonRestrictPars* par,
		int& line_id_max_pos,int& cdp_max_pos,float& t01_max,float& t02_max)
{
	float delta_t0 = 0.0;
	if(idx == 0) {
		_hor_surface_data1 = NULL;

		FILE* fp = NULL;
		fp = fopen(_hor_surface_data_file_name[0],"rb+");
		fread(&_hor_surface_data2_num1,sizeof(int),1,fp);
		fread(&_hor_surface_data2_num2,sizeof(int),1,fp);
		_hor_surface_data2 = new float[_hor_surface_data2_num1 *
			_hor_surface_data2_num2];

		fread(_hor_surface_data2,sizeof(float),
			_hor_surface_data2_num1 * _hor_surface_data2_num2,fp);

		fclose(fp);

		_hor_surface_data1_num1 = _hor_surface_data2_num1;
		_hor_surface_data1_num2 = _hor_surface_data2_num2;
		float t1 = par->_t0_depth;
		t01_max = t1;
		int i = 0,j = 0;
		float* temp = NULL;
		for(i = 0;i<_hor_surface_data2_num1;i++) {
			temp = _hor_surface_data2 + i * _hor_surface_data2_num2;
			for(j = 0;j<_hor_surface_data2_num2;j++) {
				if( (*temp > 0)
					&&((*temp - t1) > delta_t0) ) {
					line_id_max_pos = _slid + i * _lid_inv;
					t02_max = *temp;
					delta_t0 = *temp - t1;
					cdp_max_pos = _scdp + j * _cdp_inv;
				}
				temp++;
			}
		}
		delta_t0 = delta_t0 / ((float)par->_layer_num);
		return delta_t0;
	}
	if(idx == num - 1)
	{
		FILE* fp = NULL;
		fp = fopen(_hor_surface_data_file_name[idx - 1],"rb+");

		fread(&_hor_surface_data1_num1,sizeof(int),1,fp);
		fread(&_hor_surface_data1_num2,sizeof(int),1,fp);
		_hor_surface_data1 = new float[_hor_surface_data1_num1 *
			_hor_surface_data1_num2];
		fread(_hor_surface_data1,sizeof(float),
			_hor_surface_data1_num1 * _hor_surface_data1_num2,fp);
		fclose(fp);

		_hor_surface_data2_num1 = _hor_surface_data1_num1;
		_hor_surface_data2_num2 = _hor_surface_data1_num2;
		_hor_surface_data2 = NULL;

		float t1 = par->_t0_depth;
		t02_max = t1;
		int i = 0,j = 0;
		float* temp = NULL;
		for(i = 0;i<_hor_surface_data1_num1;i++) {
			temp = _hor_surface_data1 + i * _hor_surface_data2_num2;
			for(j = 0;j<_hor_surface_data1_num2;j++) {
				if( (*temp > 0)
					&&((t1 - *temp) > delta_t0) ){
					line_id_max_pos = _slid + i * _lid_inv;
					cdp_max_pos = _scdp + j * _cdp_inv;
					t01_max = *temp;
					delta_t0 = t1 - *temp;
				}
				temp++;
			}
		}
		delta_t0 = delta_t0 / ((float)par->_layer_num);

		return delta_t0;
	}
	FILE* fp = NULL;
	fp = fopen(_hor_surface_data_file_name[idx - 1],"rb+");

	fread(&_hor_surface_data1_num1,sizeof(int),1,fp);
	fread(&_hor_surface_data1_num2,sizeof(int),1,fp);
	_hor_surface_data1 = new float[_hor_surface_data1_num1 *
		_hor_surface_data1_num2];
	fread(_hor_surface_data1,sizeof(float),
		_hor_surface_data1_num1 * _hor_surface_data1_num2,fp);
	fclose(fp);

	fp = fopen(_hor_surface_data_file_name[idx],"rb+");
	fread(&_hor_surface_data2_num1,sizeof(int),1,fp);
	fread(&_hor_surface_data2_num2,sizeof(int),1,fp);
	_hor_surface_data2 = new float[_hor_surface_data2_num1 *
	_hor_surface_data2_num2];
	fread(_hor_surface_data2,sizeof(float),
		_hor_surface_data2_num1 * _hor_surface_data2_num2,fp);
	fclose(fp);

	int i = 0,j = 0;
	float* temp1 = NULL;
	float* temp2 = NULL;
	for(i = 0;i<_hor_surface_data1_num1;i++) {
		temp1 = _hor_surface_data1 + i * _hor_surface_data2_num2;
		temp2 = _hor_surface_data2 + i * _hor_surface_data2_num2;
		for(j = 0;j<_hor_surface_data1_num2;j++) {
			if( (*temp1 > 0)&&(*temp2 > 0)
				&&( ( *temp2 - * temp1)> delta_t0) ){
				line_id_max_pos = _slid + i * _lid_inv;
				cdp_max_pos = _scdp + j * _cdp_inv;
				t01_max = *temp1;
				t02_max = *temp2;
				delta_t0 = t02_max - t01_max;
			}
			temp1++;
			temp2++;
		}
	}
	delta_t0 = delta_t0 / ((float)par->_layer_num + 1);
	return delta_t0;
}

void SeisInversion::GetLineIDAndCDPForInterpolate(int& line_id,int& cdp,float& t0,wURHorizonRestrictPars* par,
	float delta_t0,int line_id_max_pos,int cdp_max_pos,float t01_max,float t02_max,int j,int k,int m)
{
	line_id = _slid_proceed + _lid_inv * k;
	cdp = _scdp_proceed + _cdp_inv * m;

	float hor_max_pos1,hor_max_pos2;
	float hor_cur_pos1,hor_cur_pos2;

	int kk = 0;
	int mm = 0;
	kk = (line_id - _slid)/_lid_inv;
	mm = (cdp - _scdp)/_cdp_inv;

	if(_hor_surface_data1) {
		hor_cur_pos1 = _hor_surface_data1[kk * _hor_surface_data1_num2 + mm];
	} else {
		hor_cur_pos1 = par->_t0_depth;
	}
	
	if(_hor_surface_data2) {
		hor_cur_pos2 = _hor_surface_data2[kk * _hor_surface_data1_num2 + mm];
	} 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::GetCurveT0(float& hor_cur_pos1,float& hor_cur_pos2,int line_id,int cdp,wURHorizonRestrictPars* par)
{
	int k = 0,m = 0;

	k = (line_id - _slid) / _lid_inv;
	m = (cdp - _scdp) / _cdp_inv;

	if(_hor_surface_data1) {
		hor_cur_pos1 = _hor_surface_data1[k * _hor_surface_data1_num2 + m];
	} else {
		hor_cur_pos1 = par->_t0_depth;
	}
	
	if(_hor_surface_data2) {
		hor_cur_pos2 = _hor_surface_data2[k * _hor_surface_data1_num2 + m];
	} else {
		hor_cur_pos2 = par->_t0_depth;
	}

	return;
}
wDBInversionCurveInSurvey** SeisInversion::GetInterpolateCurvesForLayer(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;
		GetCurveT0(hor_cur_pos1,hor_cur_pos2,line_id1,cdp1,par);
		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::GetInterpolateValue(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,int*& interpolate_lids)
{
	float t0 = 0.0;
	int i = 0,j = 0;
	float hor_cur_pos1,hor_cur_pos2;