seisinversionproceeds1.cpp

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

// SeisInversion.cpp: implementation of the TravInter class.
//
//////////////////////////////////////////////////////////////////////


#include "SeisInversion.h"
#include "newsGeneralHeader.hh"
#include "wDGColorPlateDialog.hh"
#include "wDBScatterCurveInSurvey.hh"
#include "wURSeis3DTemplate.hh"
#include "wURSeis3D.hh"
#include "wDBTempHorDataVolume.hh"
#include "wPRSeisAttrsBase.hh"

#include "stdio.h"
#include "string.h"
#include "stdlib.h"
#include "SeisInversion.h"
#include "SeisInversionDlg1.h"
#include "X11/Xos.h"
#include "time.h"
#include "jqconvlv.h"

void SeisInversion::OnProceedInversion()
{
	RemoveMouseEvent(-1);
	UnloadPopupMenu(-1);
	char inter_version[18];
	if(_seis_tpt->GetSurveyType() == Seis3Dsurvey){
		_seis_tpt->Obtain3DTemplate()->ObtainControlPar()->GetVersion(inter_version);
	}
	if(_seis_tpt->GetSurveyType() == Seis2Dsurvey){
		_seis_tpt->Obtain2DTemplate()->ObtainControlPar()->GetVersion(inter_version);
	}
	SeisInversionProceedParsDialog1* dlg = new SeisInversionProceedParsDialog1(this,inter_version);
	dlg->Manage();
	if(dlg->GetUserAnwser() == 1) {
		dlg->GetPars(_inversion_flag,_Vlog_flag,_Vlog_vmark,_wavelet_id,_Vlog_flag2,_min_freq_value,_Vlog_vmark2);
		if(_seis_tpt->GetSurveyType() == Seis3Dsurvey){
			_seis_tpt->Obtain3DTemplate()->ObtainSeis3D()->GetTimeRange(_st0,_et0,_dt0);
			_seis_tpt->Obtain3DTemplate()->ObtainSeis3D()->GetLineRange(_slid_proceed,_elid_proceed,_lid_inv);
			_seis_tpt->Obtain3DTemplate()->ObtainSeis3D()->GetCDPRange(_scdp_proceed,_ecdp_proceed,_cdp_inv);
			delete dlg;
			SetCurrentCursor(WatchShapeCursor);
			if(_inversion_flag == 1){
				if(strlen(_wavelet_id) == 0){
					delete dlg;
					return;
				}
				V_logBatchProcess1();
			}
			if(_inversion_flag == 2)
				V_logBatchProcess2();
			RestoreDefaultCursor();
		}
		if(_seis_tpt->GetSurveyType() == Seis2Dsurvey) {
			_seis_tpt->Obtain2DTemplate()->ObtainSeis2D()->GetT0DepthRange(_st0,_et0,_dt0);
			_seis_tpt->Obtain2DTemplate()->ObtainSeis2D()->GetCDPRange(_scdp_proceed,_ecdp_proceed,_cdp_inv);
			delete dlg;
			SetCurrentCursor(WatchShapeCursor);
			if(_inversion_flag == 1){
				V_logBatch2DProcess1();
				if(strlen(_wavelet_id) == 0){
					delete  dlg;
					return;
				}
			}
			if(_inversion_flag == 2)
				V_logBatch2DProcess2();
			RestoreDefaultCursor();
		}
	}
}

void V_logPreBatchProcess(float* kdimdata,int t0s_num,float max_value,float min_value)
{
	int i = 0;
	float max_value1,min_value1;
	max_value1 = (abs((int)min_value) + abs((int)max_value) ) / 2.0;
	min_value1 = -1.0 * max_value1;

	for(i = 0;i<t0s_num;i++){
		if(kdimdata[i] < min_value1){
			kdimdata[i] = -1.0;
			continue;
		}
		if(kdimdata[i] > max_value1){
			kdimdata[i] = 1.0;
			continue;
		}
		kdimdata[i] = 2.0 * (kdimdata[i] - min_value)/(max_value - min_value) - 1.0;
	}
}


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

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

	_seis_tpt->Obtain3DTemplate()->SetSeisData();
	_seis_tpt->Obtain3DTemplate()->GetMinAndMax(min_value,max_value);
	int i2,i3,ist,iet,il=61;
	int kcdpi , knumpoint , kpoint,kline,status;
	float ymax;
	float *kdimdata , *bufx , *bufo , *lfvdata;
	int init_status;
	kcdpi = (_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1;
	knumpoint = (_et0 - _st0) / _dt0 + 1;
	_t0s_num = knumpoint;
	kpoint = knumpoint;
	kcdpi = knumpoint;
	bufx = new float[kcdpi];
	bufo = new float[kcdpi];
	for(int kkk = 0;kkk < kcdpi;kkk++) {
		bufo[kkk] = 0.0;
		bufx[kkk] = 0.0;
	}
	lfvdata = new float[kcdpi];
	
	wURSeis3DTemplate* temp_seis_3d_lfv = new wURSeis3DTemplate();
	temp_seis_3d_lfv->SetUserID(_seis_tpt->Obtain3DTemplate()->GetUserID());
	temp_seis_3d_lfv->SetDataPars(_seis_tpt->Obtain3DTemplate()->GetSurveyID(),
			"Lfv",_seis_tpt->Obtain3DTemplate()->GetDataDomain());
	temp_seis_3d_lfv->SetSeisData();
	float min1,max1;
	temp_seis_3d_lfv->GetMinAndMax(min1,max1);
	delete temp_seis_3d_lfv;

        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* rimp3d_volume_data = NULL;
	wDBSeis3DVolumeData* rimp3d_volume_data1 = NULL;
	if(_Vlog_flag == 1) {
		rimp3d_volume_data = _seis_attr_pro->Create3DVolume("Rimp");
		rimp3d_volume_data1 = _seis_attr_pro->Create3DVolume("AbsImp");
	}

	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* rimp3d_data = new float[((_ecdp_proceed - _scdp_proceed)/_lid_inv + 1) * _t0s_num];
	float* rimp3d_data1 = new float[((_ecdp_proceed - _scdp_proceed)/_lid_inv + 1) * _t0s_num];

	wDBSeisWavelet* wavelet = new wDBSeisWavelet();
	char inter_version[18];
	_seis_tpt->Obtain3DTemplate()->ObtainControlPar()->GetVersion(inter_version);
	wavelet->SetInterVersion(inter_version);
	wavelet->OpenWavelet(_wavelet_id);
	float * wavelet_vals = NULL;
	int vals_len = 0;
	wavelet_vals = wavelet->GetValues(vals_len);

	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;
	}


	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);

		int temp = (kcdpi - _scdp_proceed)/_cdp_inv;
		memcpy(kdimdata + 50 / _dt0,seis3d_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));
			if(_Vlog_flag == 1) {
				memcpy(rimp3d_data + temp * _t0s_num,bufo,_t0s_num * sizeof(float));
				memcpy(rimp3d_data1 + temp * _t0s_num,bufx,_t0s_num * sizeof(float));
			}
			continue;
		}
		V_logPreBatchProcess(kdimdata  + 50 / _dt0,_t0s_num,max_value,min_value);

		float* convlv_data = NULL;
		convlv_data = convlv(kdimdata,kdimdata_num,wavelet_vals,vals_len,-1);
		memcpy(bufx,convlv_data,_t0s_num * sizeof(float));
		delete convlv_data;

                for(int pppm = 0;pppm < _t0s_num;pppm++) {
                        bufx[pppm] = bufx[pppm] * 1000.0;
                }
		float* result = NULL;
		result = SeisDataFilter(bufx,kpoint,_dt0,0.0,80.0,0.0,0.0);
		memcpy(bufx,result,_t0s_num * sizeof(float));
		delete result;
		memcpy(lfvdata,lfv3d_data + temp * _t0s_num,_t0s_num * sizeof(float));
		bufo[0] = min1;
                for(i2 = 1;i2 < kpoint;i2++){
                        bufo[i2] = (1000.0 + bufx[i2 - 1]) / (1000.0 - bufx[i2 - 1]) * bufo[i2 - 1];
                }
                memcpy(bufx,bufo,_t0s_num * sizeof(float));

		result = NULL;
		result = SeisDataFilter(bufx,_t0s_num,_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(rimp3d_data + temp * _t0s_num,bufo,_t0s_num * sizeof(float));
			memcpy(rimp3d_data1 + temp * _t0s_num,bufx,_t0s_num * sizeof(float));
		}
		memcpy(imp3d_data + temp * _t0s_num,lfvdata,_t0s_num * sizeof(float));
		}
	if(_Vlog_flag == 1) {
		rimp3d_volume_data->WriteInLine(kline,rimp3d_data);
		rimp3d_volume_data1->WriteInLine(kline,rimp3d_data1);
	}
	imp3d_volume_data->WriteInLine(kline,imp3d_data);
	}
	if(_Vlog_flag == 1) {
		delete rimp3d_volume_data;
		delete rimp3d_volume_data1;
	}
	delete rimp3d_data;
	delete rimp3d_data1;
	delete imp3d_data;
	delete imp3d_volume_data;
	delete lfv3d_volume_data;
	delete seis3d_volume_data;

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


int SeisInversion::V_logBatch2DProcess1()
{
        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;

        _seis_tpt->Obtain2DTemplate()->GetMinAndMax(min_value,max_value);
	int i2,i3,ist,iet,il=61;
	int kcdpi , knumpoint , kpoint,kline,status;
	float ymax;
	float *kdimdata , *bufx , *bufo , *lfvdata;
	int init_status;
	kcdpi = (_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1;
	knumpoint = (_et0 - _st0) / _dt0 + 1;
	_t0s_num = knumpoint;
	kpoint = knumpoint;
	kcdpi = knumpoint;
	kdimdata = NULL;
	bufx = new float[kcdpi];
	bufo = new float[kcdpi];
	for(int kkk = 0;kkk < kcdpi;kkk++) {
		bufo[kkk] = 0.0;
		bufx[kkk] = 0.0;
	}
	lfvdata = new float[kcdpi];
	
	wURSeis2DTemplate* temp_seis_2d_lfv = new wURSeis2DTemplate();
	temp_seis_2d_lfv->SetUserID(_seis_tpt->Obtain2DTemplate()->GetUserID());
	temp_seis_2d_lfv->SetDataPars(_seis_tpt->Obtain2DTemplate()->GetSurveyID(),
		_seis_tpt->Obtain2DTemplate()->GetLineID(),
		"Lfv",_seis_tpt->Obtain2DTemplate()->GetDataDomain());
	temp_seis_2d_lfv->SetSeisData();
	float min1,max1;
	temp_seis_2d_lfv->GetMinAndMax(min1,max1);
	delete temp_seis_2d_lfv;

        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* rimp2d_section_data = NULL;
	wDBSeis2DSection* rimp2d_section_data1 = NULL;
	wDBSeis2DSection* rimp2d_section_data2 = NULL;
	if(_Vlog_flag == 1) {
		rimp2d_section_data = _seis_attr_pro->Create2DSection("Rimp");
		rimp2d_section_data1 = _seis_attr_pro->Create2DSection("AbsImp");
		rimp2d_section_data2 = _seis_attr_pro->Create2DSection("ReflectFactor");
	}

	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* rimp2d_data = new float[((_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1) * _t0s_num];
	float* rimp2d_data1 = new float[((_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1) * _t0s_num];
	float* rimp2d_data2 = 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);

	wDBSeisWavelet* wavelet = new wDBSeisWavelet();
	char inter_version[18];
	_seis_tpt->Obtain2DTemplate()->ObtainControlPar()->GetVersion(inter_version);
	wavelet->SetInterVersion(inter_version);
	wavelet->OpenWavelet(_wavelet_id);
	float * wavelet_vals = NULL;
	int vals_len = 0;
	wavelet_vals = wavelet->GetValues(vals_len);

	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;
	}

	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);
		int temp = (kcdpi - _scdp_proceed)/_cdp_inv;
		memcpy(kdimdata + 50 / _dt0,seis2d_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];
		}
		memcpy(lfvdata,lfv2d_data + temp * _t0s_num,_t0s_num * sizeof(float));
		if(sum <= 0.0) {
			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(rimp2d_data2 + temp * _t0s_num,bufo,_t0s_num * sizeof(float));
			}
			memcpy(imp2d_data + temp * _t0s_num,lfvdata,_t0s_num * sizeof(float));
			continue;
		}
		V_logPreBatchProcess(kdimdata + 50 / _dt0,_t0s_num,max_value,min_value);
		float* convlv_data = NULL;
		convlv_data = convlv(kdimdata,kdimdata_num,wavelet_vals,vals_len,-1);
		memcpy(bufx,convlv_data,_t0s_num * sizeof(float));
		delete convlv_data;
		for(int pppm = 0;pppm < _t0s_num;pppm++) {
			bufx[pppm] = bufx[pppm] * 1000.0;
		}
		float* result = NULL;
		result = SeisDataFilter(bufx,kpoint,_dt0,0.0,80.0,0.0,0.0);
		memcpy(bufx,result,_t0s_num * sizeof(float));
		delete result;
		
		if(_Vlog_flag == 1){
			memcpy(rimp2d_data2 + temp * _t0s_num,bufx,_t0s_num * sizeof(float));
		}
		
		bufo[0] = min1;
		for(i2 = 1;i2 < kpoint;i2++){
			bufo[i2] = (1000.0 + bufx[i2 - 1]) / (1000.0 - bufx[i2 - 1]) * bufo[i2 - 1];
		}