?? forward.cpp
字號(hào):
//=======================================================//
//函數(shù)名稱：Forward(int LayerNum,double *LayerResistivty,//
// double *LayerThickness,double *Rs, //
// int min,int max,int sampling) //
//函數(shù)目的：用濾波系數(shù)法進(jìn)行一維正演模擬 //
//參數(shù)說(shuō)明： LayerNum：層數(shù) //
// LayerResistivty：存儲(chǔ)每層的電阻率 //
// LayerThickness: 存儲(chǔ)每層的厚度 //
// Rs：存儲(chǔ)正演的理論視電阻率 //
// min：正演理論值的最小極距的序號(hào) //
// max: 正演理論值的最小極距的序號(hào) //
//=======================================================//
#include"math.h"
void Forward(int LayerNum,double *LayerResistivty,double *LayerThickness,
double *Rs,int min,int max,int sampling)
{
int totalnum=max-min+1;
///////////////////////////////////////////////////////////////////////////
//濾波系數(shù)
//采樣率為1/3lg10的濾波系數(shù)
double FactorThree[9]={0.0148,-0.0814,0.4018,-1.5716,1.9720,
0.1854,0.1064,-0.0499,0.0225};
//采樣率為1/4lg10的濾波系數(shù)
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的濾波系數(shù)
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的濾波系數(shù)
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的濾波系數(shù)
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的濾波系數(shù)
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};
/////////////////////////////////////////////////////////////////
//根據(jù)濾波系數(shù)的個(gè)數(shù)動(dòng)態(tài)定義數(shù)組
int FactorNum,TNum,End;
switch(sampling)
{
case 6: //采樣率為1/6lg10
FactorNum=20; //濾波系數(shù)的個(gè)數(shù)
TNum=totalnum+FactorNum; //核函數(shù)值的個(gè)數(shù)
End=totalnum+5;
break;
case 10: //采樣率為1/10lg10
FactorNum=141; //濾波系數(shù)的個(gè)數(shù)
TNum=totalnum+FactorNum; //核函數(shù)值的個(gè)數(shù)
End=totalnum+101;
break;
case 4: //采樣率為1/4lg10
FactorNum=15; //濾波系數(shù)的個(gè)數(shù)
TNum=totalnum+FactorNum; //核函數(shù)值的個(gè)數(shù)
End=totalnum+2;
break;
case 3: //采樣率為1/3lg10
FactorNum=9; //濾波系數(shù)的個(gè)數(shù)
TNum=totalnum+FactorNum; //核函數(shù)值的個(gè)數(shù)
End=totalnum+3;
break;
}
int i,j;
double *T=new double[TNum];
/////////////////////////////////////////////////////////////////
///開始電阻率正演模擬
switch(sampling)
{
case 6:
/////////////////////////////////////////////////////////////
//計(jì)算核函數(shù)
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;
}
///////////////////////////////////////////////////////////////
//計(jì)算理論視電阻率
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:
/////////////////////////////////////////////////////////////
//計(jì)算核函數(shù)
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;
}
///////////////////////////////////////////////////////////////
//計(jì)算理論視電阻率
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:
/////////////////////////////////////////////////////////////
//計(jì)算核函數(shù)
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;
}
///////////////////////////////////////////////////////////////
//計(jì)算理論視電阻率
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:
/////////////////////////////////////////////////////////////
//計(jì)算核函數(shù)
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;
}
///////////////////////////////////////////////////////////////
//計(jì)算理論視電阻率
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;
}
//=======================================================//
//函數(shù)名稱：ReadModelFile() //
//函數(shù)目的：讀取正演模型文件，調(diào)用其它函數(shù)進(jìn)行正演模擬 //
//=======================================================//
#include "iostream.h"
#include "iomanip.h"
#include "fstream.h"
int ReadModelFile(int sample)
{
int i;
const int SoundNum=100; //正演理論值的個(gè)數(shù)
double *Rs=new double[SoundNum];
double *AB=new double[SoundNum];
////////////////////////////////////////////////////////////
//讀取高程數(shù)據(jù)文件
char file1[20];
cout<<endl;
cout<<" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
cout<<" !! !!"<<endl;
cout<<" !! 請(qǐng)輸入正演數(shù)據(jù)文件名: !!"<<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<<" !! 此文件格式不對(duì)或不存在 !!"<<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;
}
///////////////////////////////////////////////////////////////////
//調(diào)用正演模擬函數(shù)
Forward(LayerNum,LayerR,LayerT,Rs,Min,Max,sample);
////////////////////////////////////////////////////////////
//輸出修改后的實(shí)測(cè)數(shù)據(jù)到文件
char file2[20];
cout<<endl;
cout<<" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
cout<<" !! !!"<<endl;
cout<<" !! 請(qǐng)輸入輸出數(shù)據(jù)文件名: !!"<<endl;
cout<<" !! !!"<<endl;
cout<<" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
cin>>file2;
ofstream out;
out.open(file2,ios::out||ios::right);
out<<" 電 測(cè) 深 正 演 結(jié) 果 "<<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;
}
//=======================================================//
//函數(shù)名稱：main() //
//函數(shù)目的：進(jìn)行一維電測(cè)深反演 //
//程序作者：劉海飛 //
//編程時(shí)間：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<<" !! 請(qǐng)輸入:<<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;
//////////////////////////////////////////////////////////
//程序結(jié)束
default:
return;
}
}
return ;
}
?? 快捷鍵說(shuō)明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -