controlwnd.cpp

是有关石油物理测井曲线的分层方法

	case IDNO:
		for(kk=0;kk<wh.nTotalNumPoint;kk++)
			wh.pLogData[kk][nCurId]=wh.pLogData[kk][wh.tnNumCur];
		this->m_pFr->SendMessage(WM_CLOSE);	
		break;
	case IDYES:
		OnFileSave();
		this->m_pFr->SendMessage(WM_CLOSE);	
		break;
	default://IDCANCEL
		return;
	}

}
void CControlWnd::OnSelchangeCmbCurnm() 
{
	// TODO: Add your control notification handler code here
	flag=false;	
}

//极值方差法分层
void CControlWnd::JZFCFC(int cur, int nStPoint, int nendPoint, int w)
{
	int i,j,k;
	int nPoint=nendPoint-nStPoint;
	//对测井曲线进行[0-1]归一化处理
	float max,min;//,pData[1100];
//	float *Q=new float[wh.nTotalNumPoint+1];
	float *Q=new float[wh.nTotalNumPoint+1];
	float *pData=new float[wh.nTotalNumPoint+1];
	max=min=wh.pLogData[nStPoint][cur];
	for(i=nStPoint;i<nendPoint;i++)
	{
		if(max<wh.pLogData[i][cur])
			max=wh.pLogData[i][cur];
		if(min>wh.pLogData[i][cur])
		{
			min=wh.pLogData[i][cur];
		}
	}
	float tt=max-min;
	for(i=nStPoint;i<nendPoint;i++)
		pData[i]=(wh.pLogData[i][cur]-min)/tt;

	//求取分层因数Q(n)
	float Tz=0,T1=0,T2=0;
	//float Q[1100];

//	int nBed=0,*pBed=new int[nPoint+1];
//	int nBed=0,pBed[60],pBed1[60];
	nBed=0;
	for(i=nStPoint;i<nendPoint;i++)
		Tz+=pData[i];
	for(i=nStPoint;i<nendPoint;i++)
	{
		T1+=pData[i];
		T2=Tz-T1;
		Q[i]=T1*T1/(float)(i-nStPoint+1)+T2*T2/(float)(nendPoint-i-1);
//		if(i>18900)
//			int aas=0;
	}
	//根据分层因数进行粗分层(分成nBed层)
	pBed[nBed++]=nStPoint;
	for(i=nStPoint+1;i<nendPoint;i++)
	{
		if((Q[i]>Q[i-1] && Q[i]>Q[i+1]))//||(Q[i]<Q[i-1] && Q[i]<Q[i+1]))
		{
			pBed[nBed++]=i;
		}
	}
	pBed[nBed]=nendPoint-1;
	//根据方差分层指标进行细分层
	int m=4,m1,m2;
	float *R=new float[wh.nTotalNumPoint+1];	
	float RR[9];
	float A,B;
	R[0]=1.0;

//Rmax细分层
	for(i=1;i<nBed;i++)
	{
		//求层间方差
		float Tt=0,Tp,Tp1,Tp2;//相邻两层测井总值、平均值,各层均值
		for(k=pBed[i-1];k<pBed[i+1];k++)
			Tt+=pData[k];
		Tp=Tt/(pBed[i+1]-pBed[i-1]);		
		B=0.0;	//相邻两层总方差
		for(k=pBed[i-1];k<pBed[i+1];k++)
			B+=(pData[k]-Tp)*(pData[k]-Tp);
		m=4;
		m1=m2=4;
		if(m1>pBed[i]-pBed[i-1]-1)
			m1=pBed[i]-pBed[i-1]-1;
		if(m2>pBed[i+1]-pBed[i]-1)
			m2=pBed[i+1]-pBed[i]-1;

		int kk=0;
		for(j=pBed[i]-m1;j<pBed[i]+m2;j++)
		{	
			T1=0;
			for(k=pBed[i-1];k<j;k++)
				T1+=pData[k];
			Tp1=T1/(j-pBed[i-1]);
			Tp2=(Tt-T1)/(pBed[i+1]-j);
			//求层内方差
			A=0;
			for(k=pBed[i-1];k<j;k++)
				A+=(pData[k]-Tp1)*(pData[k]-Tp1);
			for(k=j;k<pBed[i+1];k++)
				A+=(pData[k]-Tp2)*(pData[k]-Tp2);
			A=A/(pBed[i+1]-pBed[i-1]);
			RR[kk++]=1-A/B;
		}
		int kkk=0;
		max=RR[0];
		for(j=1;j<m2+m1;j++)
		{
			if(max<RR[j])
			{
				max=RR[j];
				kkk=j;
			}
		}
		R[i]=max;
		pBed[i]=pBed[i]-m1+kkk;
	}
	//根据方差确定实际分层界线值(A-层内，B-层间)
	float Tp=Tz/nPoint,Xp[1000];
	//求层内平均值Xp[i]
	for(i=0;i<nBed;i++)
	{
		Xp[i]=0;
//		if(i>1020)
//			int aaaa=0;
		for(j=pBed[i];j<pBed[i+1];j++)
			Xp[i]+=pData[j];
		Xp[i]/=(pBed[i+1]-pBed[i]+1);
	}
	//求层内方差A
	A=0;
	for(i=0;i<nBed;i++)
	{
		for(j=pBed[i];j<pBed[j+1];j++)
		{
			A+=ECF(pData[j]-Xp[i]);
		}
	}
	A/=(nPoint-nBed);
	//求层间方差B
	B=0;
	for(i=0;i<nBed;i++)
	{
		B+=(pBed[j+1]-pBed[i]-1)*ECF(Xp[i]-Tp);
	}
	B/=(nBed-1);
	float Ra=1-A/B;

	//对测井曲线进行分层，求取各层点的值，组成新的分层曲线
	float TT;
	for(i=0;i<=nBed;i++)
	{
		TT=0;
		for(j=pBed[i];j<pBed[i+1];j++)
			TT+=wh.pLogData[j][cur]; 
		TT=TT/(pBed[i+1]-pBed[i]);
		for(j=pBed[i];j<pBed[i+1];j++)
			wh.pLogData[j][cur]=TT;
	}

//输出验证
/*	CString str,str1;
	
	CFile ff("f:\\out.dat", CFile::modeCreate | CFile::modeWrite );
	for(i=0;i<nendPoint;i++)
	{
		str.Format("%10.3f %10.6f\n",wh.pLogData[i][0],Q[i]);
		ff.Write(str,str.GetLength());
	}
	ff.Close(); 

	str=wh.FileName;
	str=wh.FileName ;
	i=str.GetLength();
	str1="_BED_"+ wh.Curve[cur].strCurName; 
	str.Insert(i-4,str1); 

	CFile f(str, CFile::modeCreate | CFile::modeWrite );
	str.Format("%s\n","分层方法: 方差极值法");
	f.Write(str,str.GetLength());
	str.Format("%s","总层数 ：");
	f.Write(str,str.GetLength());
	str.Format("%4d\n\n",nBed);
	f.Write(str,str.GetLength());	

	str.Format("%s"," 序号  "); 
	f.Write(str,str.GetLength()); 
	str.Format("%s","顶界面深度    "); 
	f.Write(str,str.GetLength());
	str.Format("%s","底界面深度    "); 
	f.Write(str,str.GetLength());
	str.Format("%s","层厚度    ");
	f.Write(str,str.GetLength());	
	str.Format("%s","Q 极值    "); 
	f.Write(str,str.GetLength());
	str.Format("%s","R 极值    "); 
	f.Write(str,str.GetLength());
	str.Format("%s",wh.Curve[cur].strCurName); 
	f.Write(str,str.GetLength());
	f.Write("\n",1);

	for(i=0;i<nBed;i++)
	{
		str.Format("%4d ",i+1);
		f.Write(str,4);
		str.Format("%10.3f ",wh.pLogData[pBed[i]][0]);
		f.Write(str,str.GetLength());
		str.Format("%10.3f ",wh.pLogData[pBed[i+1]-1][0]);
		f.Write(str,str.GetLength());
		str.Format("%10.3f ",wh.pLogData[pBed[i+1]-1][0]-wh.pLogData[pBed[i]][0]);
		f.Write(str,str.GetLength());
		str.Format("%10.5f ",Q[i]);
		f.Write(str,str.GetLength());
		str.Format("%10.5f ",R[i]);
		f.Write(str,str.GetLength());
		str.Format("%10.3f\n",wh.pLogData[pBed[i]][cur]);
		f.Write(str,str.GetLength());
	}
	f.Close();
*/
	outJirGuo();
	delete []Q;
	delete []pData;
	delete []R;
}

//利用对称差斜率法进行分层处理
//n-窗长，窗内点数，sill-门限值，对称差大于该值的深度点为分层点
void CControlWnd::DCCXLFC(int cur, int nStPoint, int nendPoint, int n, float sill)
{
	int i,j,k;
	CString str,str1;
	int nPoint=nendPoint-nStPoint+1;
	//对测井曲线进行[0-1]归一化处理
	float max,min;//,pData[1100];
	float *T=new float[wh.nTotalNumPoint];
	//int *Tide=new int[1010]; 
	//int Tide[1010];
	//double T[1010];
	float *pData=new float[wh.nTotalNumPoint];
	max=min=wh.pLogData[nStPoint][cur];
	for(i=nStPoint;i<=nendPoint;i++)
	{
		if(max<wh.pLogData[i][cur])
			max=wh.pLogData[i][cur];
		if(min>wh.pLogData[i][cur])
			min=wh.pLogData[i][cur];
	}
	float tt=max-min;
	for(i=nStPoint;i<=nendPoint;i++)
		pData[i]=(wh.pLogData[i][cur]-min)/tt;
	//逐点求取对称差斜率T=tt*A/B和趋势码Tide
	T[0]=sill;
	tt=n*(n+1)*wh.RLEV; 
	double A=0,B=0;
	float pi=180/3.14159265;
	for(i=1;i<=n;i++)
		B+=(n-i)*(n-i);
	B/=n;
	for(k=nStPoint+n/2;k<nendPoint-n/2;k++)
	{
		A=0;
		for(i=1;i<=n/2;i++)
			A+=(pData[k+i]-pData[k-i])*(n-i)*(n-i);
		A/=n;
		T[k]=fabs((50*A/B)/tt);
		T[k]=atan(T[k]);
		T[k]*=pi;
	}
/*	CFile ff("f:\\oo.dat", CFile::modeCreate | CFile::modeWrite );
	str.Format("%s"," 序号  "); 
	ff.Write(str,str.GetLength()); 
	str.Format("%s"," 深度  "); 
	ff.Write(str,str.GetLength()); 
	str.Format("%s","对称差斜率值T  "); 
	ff.Write(str,str.GetLength());
	for(i=nStPoint;i<wh.nTotalNumPoint;i++)
	{
		str.Format("%4d ",i+1);
		ff.Write(str,4);
		str.Format("%10.3f ",wh.pLogData[i][0] );
		ff.Write(str,str.GetLength());
		str.Format("%10.3f\n",T[i-nStPoint]);
		ff.Write(str,str.GetLength());
	}
	ff.Close();
*/

//根据分层因数进行粗分层(分成nBed层)
//	int nBed=0,*pBed=new int[nPoint];
//	int nBed=0,pBed[1000],pBed1[1000];
	nBed=0;
	pBed[nBed++]=nStPoint;
	for(i=nStPoint+n/2+1;i<nendPoint-n/2;i++)
	{
		if(T[i]>sill && (T[i]>T[i-1] && T[i]>=T[i+1]))
		{
			pBed[nBed++]=i;
		}
	}
	pBed[nBed]=nendPoint;
	//对测井曲线进行分层，求取各层点的值，组成新的分层曲线
	float TT[81];
	for(i=0;i<=nBed;i++)
	{
		for(k=1;k<=wh.nNumCur;k++)
		{
			TT[k]=0;
			for(j=pBed[i];j<pBed[i+1];j++)
				TT[k]+=wh.pLogData[j][k]; 
			TT[k]=TT[k]/(pBed[i+1]-pBed[i]);
			for(j=pBed[i];j<pBed[i+1];j++)
				wh.pLogData[j][k+wh.nNumCur]=TT[k];
		}
	}
/*
//输出验证
	str=wh.FileName;
	i=str.GetLength();
	str1="_BED_"+ wh.Curve[cur].strCurName; 
	str.Insert(i-4,str1); 
	CFile f(str, CFile::modeCreate | CFile::modeWrite );
	str.Format("%s","窗长: ");
	f.Write(str,str.GetLength());
	str.Format("%4d\n",n);
	f.Write(str,str.GetLength());
	str.Format("%s","截止值 ："); 
	f.Write(str,str.GetLength());
	str.Format("%4.2f\n",sill);
	f.Write(str,str.GetLength());
	str.Format("%s","总层数 ：");
	f.Write(str,str.GetLength());
	str.Format("%4d\n\n",nBed);
	f.Write(str,str.GetLength());

	str.Format("%s"," 序号  "); 
	f.Write(str,str.GetLength()); 
	str.Format("%s","顶界面深度    "); 
	f.Write(str,str.GetLength());
	str.Format("%s","底界面深度    "); 
	f.Write(str,str.GetLength());
	str.Format("%s","层厚度    "); 
	f.Write(str,str.GetLength());
	str.Format("%s","对称差斜率值T  "); 
	f.Write(str,str.GetLength());
	str.Format("%s",wh.Curve[cur].strCurName); 
	f.Write(str,str.GetLength());
	f.Write("\n",1);
	for(i=0;i<nBed;i++)
	{
		str.Format("%4d ",i+1);
		f.Write(str,4);
		str.Format("%10.3f ",wh.pLogData[pBed[i]][0]);
		f.Write(str,str.GetLength());
		str.Format("%10.3f ",wh.pLogData[pBed[i+1]-1][0]);
		f.Write(str,str.GetLength());
		str.Format("%10.5f ",wh.pLogData[pBed[i+1]-1][0]-wh.pLogData[pBed[i]][0]);
		f.Write(str,str.GetLength());
		str.Format("%10.3f ",T[pBed[i]]);
		f.Write(str,str.GetLength());
		str.Format("%10.3f\n",wh.pLogData[pBed[i]][cur]);
		f.Write(str,str.GetLength());
	}
	f.Close();
*/	outJirGuo();
//	delete []pData;
	delete []T;
}

void CControlWnd::outJirGuo()
{
	int i,j,k;
	float sum,c;
	CString str,tmpStr,str1;
	str=wh.FileName;
	i=str.GetLength();
	str1="_SSBP"; 
	str.Insert(i-4,str1); 
	CFile f(str, CFile::modeCreate | CFile::modeWrite);
	//文本格式文件
	f.Write(wh.WellName,wh.WellName.GetLength());
	f.Write("\n",1);							//井名
	tmpStr.Format("%3d\n",wh.nNumCur+2);
	f.Write(tmpStr,tmpStr.GetLength());							//曲线条数
	tmpStr.Format("%d,%d\n",1,nBed);			//起止深度信息
	f.Write(tmpStr,tmpStr.GetLength());	
	tmpStr.Format("%d\n",1);							//采样间隔
	f.Write(tmpStr,tmpStr.GetLength());	
	f.Write("END\n\n",5);
	//曲线名称
	str.Format("%s"," 序号  "); 
	f.Write(str,str.GetLength()); 
	str.Format("%s","dep1    "); 
	f.Write(str,str.GetLength());
	str.Format("%s","dep2    "); 
	f.Write(str,str.GetLength());
	for(i=1;i<=wh.nNumCur;i++)
	{
		tmpStr.Format("%14s ",wh.Curve[i].strCurName) ;
		f.Write(tmpStr,tmpStr.GetLength() );
	}	
	f.Write("\n",1);
	for(i=0;i<nBed-1;i++)
	{
		str.Format("%4d ",i+1);
		f.Write(str,4);
		str.Format("%10.3f ",wh.pLogData[pBed[i]][0]);
		f.Write(str,str.GetLength());
		str.Format("%10.3f ",wh.pLogData[pBed[i+1]-1][0]);
		f.Write(str,str.GetLength());
		for(j=1;j<=wh.nNumCur;j++)
		{
			sum=0.0;
			for(k=pBed[i];k<=pBed[i+1];k++)
				sum+=wh.pLogData[k][j];
			sum/=(pBed[i+1]-pBed[i]+1);
			tmpStr.Format("%14.5f ",sum) ;
			f.Write(tmpStr,tmpStr.GetLength());
		}
		f.Write("\n",1);
	}
	f.Close() ;
}

//层内差异法进行分层，
//cur-曲线号,nStPoint-开始采样点数,nendPoint-终止采样点数,B-误差倍数参数(分层系数)
void CControlWnd::CNCYFFC(int cur, int nStPoint, int nendPoint, float B)
{
	int i,j,k;
	CString str,str1;
	int nPoint=nendPoint-nStPoint+1;
	//对测井曲线进行[0-1]归一化处理
	float max,min;//,pData[1100];
	float *pData=new float[wh.nTotalNumPoint];
	max=min=wh.pLogData[nStPoint][cur];
	for(i=nStPoint;i<=nendPoint;i++)
	{
		if(max<wh.pLogData[i][cur])
			max=wh.pLogData[i][cur];
		if(min>wh.pLogData[i][cur])
			min=wh.pLogData[i][cur];
	}
	float tt=max-min;
	for(i=nStPoint;i<=nendPoint;i++)
		pData[i]=(wh.pLogData[i][cur]-min)/tt;

//假设分成BED层(分成nBed层)
//	int nBed=0,*pBed=new int[nPoint];
//	int ,pBed[1000],pBed1[1000];
	nBed=0;
	float xp,dd,S;//层内均值，均方根差，方差
	int st,ed,n;//当前层的起止点(st<=p<=ed),总点数n=ed-st+1
	float sum=0.0;
	k=nStPoint;
	n=0;
	pBed[nBed++]=nStPoint;
	do{
		if(n==0)
		{
			st=k;
			sum=pData[k++];
		}
		sum+=pData[k++];
		n++;

		xp=sum/(n+1);
		S=0;
		for(i=st;i<k;i++)
			S+=ECF(pData[i]-xp);
		S/=(k-st-1);
		dd=sqrt(S);
		if(fabs(pData[k]-xp)>B*dd)
		{
			pBed[nBed++]=k;
			n=0;
			sum=0.0;
		}
	}
	while(k<nendPoint);
	pBed[nBed]=nendPoint;
	//对测井曲线进行分层，求取各层点的值，组成新的分层曲线
	float TT;
	for(i=0;i<=nBed;i++)
	{
		TT=0;
		for(j=pBed[i];j<pBed[i+1];j++)
			TT+=wh.pLogData[j][cur]; 
		TT=TT/(pBed[i+1]-pBed[i]);
		for(j=pBed[i];j<pBed[i+1];j++)
			wh.pLogData[j][cur]=TT;
	}
	outJirGuo();
}