forward.cpp

该程序用于地球物理电测深的正演模拟

//=======================================================//
//函数名称：Forward(int LayerNum,double *LayerResistivty,//
//			        double *LayerThickness,double *Rs,   //
//                  int min,int max,int sampling)        //
//函数目的：用滤波系数法进行一维正演模拟  	             //
//参数说明：       LayerNum：层数						 //
//			LayerResistivty：存储每层的电阻率  		     //
//			 LayerThickness: 存储每层的厚度			     //
//			             Rs：存储正演的理论视电阻率		 //
//                      min：正演理论值的最小极距的序号  //
//                      max: 正演理论值的最小极距的序号  //
//=======================================================//
#include"math.h"
void Forward(int LayerNum,double *LayerResistivty,double *LayerThickness,
			 double *Rs,int min,int max,int sampling)
{	
	int totalnum=max-min+1;
	///////////////////////////////////////////////////////////////////////////
	//滤波系数
	//采样率为1/3lg10的滤波系数
	double FactorThree[9]={0.0148,-0.0814,0.4018,-1.5716,1.9720,
				           0.1854,0.1064,-0.0499,0.0225};

	//采样率为1/4lg10的滤波系数
	double FactorFour[15]={0.0055,-0.0086,0.0136,-0.0225,0.0387,-0.0700,
						   0.1402,-0.3470,1.0570,-2.4820,1.8371,0.6448,
						   0.1686,0.0168,0.0098};
	//采样率为1/5lg10的滤波系数
	double FactorFive[13]={-0.000109308,0.001133902,0.003441750,0.019546163,
						   0.128387722,0.599970171,1.730157072,-2.028810605,
						   0.652450612,-0.125506562,0.023042366,-0.004341549,
						   0.000633436};

	//采样率为1/6lg10的滤波系数
	double FactorSix[20]={-0.000318,0.002072,-0.004978,0.01125,
						  -0.02521,0.05812,-0.1436,0.393,-1.1324,
						  2.7044,-3.4507,0.4248,1.1817,0.6194,0.2374,
						  0.08688,0.0235,0.01284,-0.001198,0.00304};

	//采样率为1/8lg10的滤波系数
	double FactorEight[15]={0.003813057,-0.012133073,0.018649510,-0.020403095,
							0.026125439,-0.008410902,0.158681610,0.511621613,
							1.897114021,-2.172056478,0.722719685,-0.153045037,
							0.034972594,-0.009234553,0.0015868564};

	//采样率为1/10lg10的滤波系数
	double FactorTen[141]={6174e-8,-12484e-8,12726e-8,-12975e-8,13231e-8,
							-13494e-8,13765e-8,14043e-8,14330e-8,-14625e-8,
							14930e-8,-15244e-8,15567e-8,-15901e-8,16246e-8,
							-16602e-8,16971e-8,-17352e-8,17746e-8,-18154e-8,
							18577e-8,-19015e-8,19469e-8,-19941e-8,20429e-8,
							-20936e-8,21463e-8,-22009e-8,22577e-8,-23166e-8,
							23779e-8,-24416e-8,25079e-8,-25768e-8,26487e-8,
						    -27235e-8,28016e-8,-28830e-8,29680e-8,-30568e-8,
							31496e-8,-32467e-8,33484e-8,-34549e-8,35666e-8,
							-36838e-8,38069e-8,-39363e-8,40724e-8,-42156e-8,
							43666e-8,-45259e-8,46940e-8,-48717e-8,50596e-8,
							-52587e-8,54697e-8,-56936e-8,59314e-8,-61845e-8,
							64540e-8,-67414e-8,70484e-8,-73767e-8,77284e-8,
							-81057e-8,85111e-8,-89475e-8,94183e-8,-99267e-8,
						    104775e-8,-110741e-8,117248e-8,-124303e-8,132085e-8,
							-140416e-8,149959e-8,-159826e-8,171917e-8,-182946e-8,
							199955e-8,-209469e-8,239052e-8,-234543e-8,304916e-8,
							-234124e-8,453990e-8,-106745e-8,899282e-8,550573e-8,
							2442523e-8,3250077e-8,7926675e-8,13023345e-8,
							25610307e-8,41150741e-8,64231809e-8,72803988e-8,
							36118538e-8,-100406442e-8,-242172543e-8,20052460e-8,
							444506381e-8,-489348908e-8,294899398e-8,-137791072e-8,
							61285163e-8,-29362551e-8,15817356e-8,-9504597e-8,
							6226174e-8,-4353505e-8,3198475e-8,-2441493e-8,1920840e-8,
						    -1548505e-8,1273595e-8,-1065148e-8,903512e-8,-775750e-8,
							673079e-8,-589375e-8,520264e-8,-462558e-8,413891e-8,
							-372478e-8,336951e-8,-306251e-8,279543e-8,-256168e-8,
							235594e-8,-217394e-8,201216e-8,-186773e-8,173826e-8,
						    -162176e-8,151657e-8,-142126e-8,133463e-8,-125568e-8,
							60905e-8};	
	/////////////////////////////////////////////////////////////////
	//根据滤波系数的个数动态定义数组
	int FactorNum,TNum,End;
	switch(sampling)
	{
	case 6:                      //采样率为1/6lg10
		FactorNum=20;            //滤波系数的个数
		TNum=totalnum+FactorNum; //核函数值的个数 	
		End=totalnum+5;
		break;
	case 10:                     //采样率为1/10lg10
		FactorNum=141;           //滤波系数的个数
		TNum=totalnum+FactorNum; //核函数值的个数 	
		End=totalnum+101;
		break;
	case 4:                      //采样率为1/4lg10
		FactorNum=15;            //滤波系数的个数
		TNum=totalnum+FactorNum; //核函数值的个数 	
		End=totalnum+2;
		break;
	case 3:                      //采样率为1/3lg10
		FactorNum=9;             //滤波系数的个数
		TNum=totalnum+FactorNum; //核函数值的个数 	
		End=totalnum+3;
		break;
	}
	int i,j;
	double *T=new double[TNum];	
	/////////////////////////////////////////////////////////////////
	///开始电阻率正演模拟
	switch(sampling)
	{
	case 6:		
		/////////////////////////////////////////////////////////////
	    //计算核函数
		for(i=-14;i<=End;i++)
		{
			int t=i+14;
			double Ttemp=LayerResistivty[LayerNum-1];
			double Ftemp=1.1396*pow(10.0,-i/6.0);
			for(j=LayerNum-2;j>=0;j--)
			{
				double temp=(Ttemp-LayerResistivty[j])
							*exp(-2*Ftemp*LayerThickness[j])
							/(Ttemp+LayerResistivty[j]);
				Ttemp=(1+temp)*LayerResistivty[j]/(1-temp);	
			}
			T[t]=Ttemp;	
		}
		///////////////////////////////////////////////////////////////
		//计算理论视电阻率
		for(i=min;i<=max;i++)
		{
			double sum=0.0;
			for(int j=0;j<FactorNum;j++)
			{
				sum=sum+T[j+i-min]*FactorSix[j];				
			}
			Rs[i-min]=sum;				
		}		
		break;
	case 10:		
		/////////////////////////////////////////////////////////////
	    //计算核函数
		for(i=-39;i<=End;i++)
		{
			int t=i+39;
			double Ttemp=LayerResistivty[LayerNum-1];
			double Ftemp=pow(10.0,-i/10.0)/exp(-1.7239458);
			for(j=LayerNum-2;j>=0;j--)
			{
				double temp=(Ttemp-LayerResistivty[j])
							*exp(-2*Ftemp*LayerThickness[j])
							/(Ttemp+LayerResistivty[j]);
				Ttemp=(1+temp)*LayerResistivty[j]/(1-temp);	
			}
			T[t]=Ttemp;
		}		
		///////////////////////////////////////////////////////////////
		//计算理论视电阻率
		for(i=min;i<=max;i++)
		{
			double sum=LayerResistivty[LayerNum-1]*FactorTen[0]
						+LayerResistivty[0]*FactorTen[140];
			for(int j=-99;j<=39;j++)
			{
				sum=sum+T[i-j+39-min]*FactorTen[100+j];				
			}		
			Rs[i-min]=sum;	
		}		
		break;
	case 4:		
		/////////////////////////////////////////////////////////////
	    //计算核函数
		for(i=-12;i<=End;i++)
		{
			int t=i+12;
			double Ttemp=LayerResistivty[LayerNum-1];
			double Ftemp=pow(10.0,-i/4.0)/exp(-0.1950);
			for(j=LayerNum-2;j>=0;j--)
			{
				double temp=(Ttemp-LayerResistivty[j])
							*exp(-2*Ftemp*LayerThickness[j])
							/(Ttemp+LayerResistivty[j]);
				Ttemp=(1+temp)*LayerResistivty[j]/(1-temp);	
			}
			T[t]=Ttemp;			
		}
		///////////////////////////////////////////////////////////////
		//计算理论视电阻率
		for(i=min;i<=max;i++)
		{
			double sum=0.0;
			for(int j=0;j<FactorNum;j++)
			{
				sum=sum+T[j+i-min]*FactorFour[j];				
			}		
			Rs[i-min]=sum;	
		}		
		break;
	case 3:		
		/////////////////////////////////////////////////////////////
	    //计算核函数
		for(i=-5;i<=End;i++)
		{
			int t=i+5;
			double Ttemp=LayerResistivty[LayerNum-1];
			double Ftemp=pow(10.0,-i/3.0)/exp(-0.0488);
				for(j=LayerNum-2;j>=0;j--)
			{
				double temp=(Ttemp-LayerResistivty[j])
							*exp(-2*Ftemp*LayerThickness[j])
							/(Ttemp+LayerResistivty[j]);
				Ttemp=(1+temp)*LayerResistivty[j]/(1-temp);	
			}
			T[t]=Ttemp;			
		}
		///////////////////////////////////////////////////////////////
		//计算理论视电阻率
		for(i=min;i<=max;i++)
		{
			double sum=0.0;
			for(int j=0;j<FactorNum;j++)
			{
				sum=sum+T[j+i-min]*FactorThree[j];				
			}		
			Rs[i-min]=sum;		
		}		
		break;
	}
	delete T;
}
//=======================================================//
//函数名称：ReadModelFile() 							 //
//函数目的：读取正演模型文件，调用其它函数进行正演模拟   //
//=======================================================//  	
#include "iostream.h"
#include "iomanip.h"
#include "fstream.h"
int ReadModelFile(int sample)
{
	int i;
	const int SoundNum=100;              //正演理论值的个数
	double *Rs=new double[SoundNum];
	double *AB=new double[SoundNum];

	////////////////////////////////////////////////////////////
	//读取高程数据文件
	char file1[20];
	cout<<endl;
	cout<<"         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
	cout<<"         !!                                !!"<<endl;
	cout<<"         !!      请输入正演数据文件名:     !!"<<endl;
	cout<<"         !!                                !!"<<endl;
	cout<<"         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
	cin>>file1;

	ifstream infile;
	infile.open(file1,ios::in|ios::nocreate);
	if(infile.fail())
	{
		cout<<endl;
		cout<<"         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
		cout<<"         !!                                !!"<<endl;
		cout<<"         !!     此文件格式不对或不存在     !!"<<endl;
		cout<<"         !!                                !!"<<endl;
		cout<<"         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
		return 0;
	}
	int LayerNum;
	float minAB,maxAB;
	infile>>minAB;
	infile>>maxAB;
	infile>>LayerNum;

	double *LayerR=new double[LayerNum];
    double *LayerT=new double[LayerNum-1];
	for(i=0;i<LayerNum;i++)
	{
		infile>>LayerR[i];
		if(i==LayerNum-1)break;	
		infile>>LayerT[i];	
	}	
	infile.close();
	int Min,Max,Number;		
	switch(sample)
	{
    ////////////////////////////////////////////////////////////
	//采样率为1/3
	case 3:
		Min=int(3*log10(minAB));
		Max=int(3*log10(maxAB));		
		Number=Max-Min+1;
		cout<<Min<<ends<<Max<<ends<<Number<<endl;
		////////////////////////////////////////////////////////////////////
		//算出电极距
		for(i=Min;i<=Max;i++)
		{
			AB[i-Min]=(pow(10.0,i/3.0));
			cout<<i-Min<<ends<<AB[i-Min]<<endl;
		}
		break;
	//////////////////////////////////////////////////////////
	//采样率为1/4
	case 4:
		Min=int(4*log10(minAB));
		Max=int(4*log10(maxAB));		
		Number=Max-Min+1;
		////////////////////////////////////////////////////////////////////
		//算出电极距
		for(i=Min;i<=Max;i++)
		{
			AB[i-Min]=(pow(10.0,i/4.0));		
		}
		break;
	//////////////////////////////////////////////////////////
	//采样率为1/6
	case 6:
		Min=int(6*log10(minAB));
		Max=int(6*log10(maxAB));
		Number=Max-Min+1;
		
		////////////////////////////////////////////////////////////////////
		//算出电极距
		for(i=Min;i<=Max;i++)
		{
			AB[i-Min]=(pow(10.0,i/6.0));			
		}
		break;
	//////////////////////////////////////////////////////////
	//采样率为1/10
	case 10:
		Min=int(10*log10(minAB));
		Max=int(10*log10(maxAB));		
		Number=Max-Min+1;
		////////////////////////////////////////////////////////////////////
		//算出电极距
		for(i=Min;i<=Max;i++)
		{
			AB[i-Min]=(pow(10.0,i/10.0));					
		}
		break;	
	}	
	///////////////////////////////////////////////////////////////////
	//调用正演模拟函数	
	Forward(LayerNum,LayerR,LayerT,Rs,Min,Max,sample);	
	////////////////////////////////////////////////////////////
	//输出修改后的实测数据到文件
	char file2[20];
	cout<<endl;
	cout<<"         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
	cout<<"         !!                                !!"<<endl;
	cout<<"         !!      请输入输出数据文件名:     !!"<<endl;
	cout<<"         !!                                !!"<<endl;
	cout<<"         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
	cin>>file2;

	ofstream out;
	out.open(file2,ios::out||ios::right);
	out<<"         电 测 深 正 演 结 果          "<<endl;
	out<<"======================================="<<endl;
	out<<setw(15)<<"电极距"	   
	   <<setw(15)<<"模拟Res"
	   <<endl<<endl;
	for(i=0;i<Number;i++)
	{
		out<<setiosflags(ios::fixed)
	       <<setprecision(2)
		   <<setw(15)<<AB[i]		   
		   <<setw(15)<<Rs[i]		   
		   <<endl;
	}
	return 0;
}	
//=======================================================//
//函数名称：main()                                       //
//函数目的：进行一维电测深反演                           //
//程序作者：刘海飞                                       //
//编程时间：2004.8                                       //
//=======================================================//

void main()
{
	int num,sample;
	cout<<endl<<endl;
	while(1)
	{
		cout<<"         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
		cout<<"         !!============正演模拟============!!"<<endl;
		cout<<"         !!           选择采样率           !!"<<endl;
		cout<<"         !!                                !!"<<endl;
		cout<<"         !!         1、采样率为1/3         !!"<<endl;
		cout<<"         !!         2、采样率为1/4         !!"<<endl;
		cout<<"         !!         3、采样率为1/6         !!"<<endl;
		cout<<"         !!         4、采样率为1/10        !!"<<endl;
		cout<<"         !!         5、返回!!!             !!"<<endl;
		cout<<"         !!                                !!"<<endl;
		cout<<"         !!         请输入:<<1—5>>        !!"<<endl;
		cout<<"         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
		cin>>num;
		switch(num)
		{
        ////////////////////////////////////////////////////////////
		//采样率为1/3
		case 1:
			sample=3;
			ReadModelFile(sample);
			break;
		//////////////////////////////////////////////////////////
		//采样率为1/4
		case 2:
			sample=4;
			ReadModelFile(sample);
			break;
		//////////////////////////////////////////////////////////
		//采样率为1/6
		case 3:
			sample=6;
			ReadModelFile(sample);
			break;
		//////////////////////////////////////////////////////////
		//采样率为1/10
		case 4:
			sample=10;
			ReadModelFile(sample);
			break;
		//////////////////////////////////////////////////////////
		//程序结束
		default:		
			return;
		}
	}
	return ;
}