井眼轨迹数据的测量值是离散的,根据某些测斜公式,我们可以计算出离散的三维的井眼轨迹坐标,但是真实的井眼轨迹是一条平滑的曲线,这就需要我们对测斜数据进行插值,使井眼轨迹变得平滑,我暂时决定使用三次样条进行插值。(但是插值出来的点,并不是真实的测量值,并不能真实的反映经验轨迹的实际情况,仅供分析使用)

  三次样条函数:(函数是在网上找到的,测试可用)

  ThreeSpline.h

#pragma once
class ThreeSpline
{
public:
ThreeSpline(void);
~ThreeSpline(void);
// 利用求出的二阶导数求给定点值(结合Spline1,Spline2)
void Splint(double *xa,double *ya,double *m,int n,double &x,double &y);
//对一系列点求二阶偏导数,点横坐标单调递增(I型边界)(结合Spline)一阶偏导数
void Spline1(double *xa,double *ya,int n,double *&m,double bound1,double bound2);
//对一系列点求二阶偏导数,点横坐标单调递增(II型边界)(结合Spline)二阶偏导数
void Spline2(double *xa,double *ya,int n,double *&m,double bound1=,double bound2=); };

ThreeSpline.cpp

#include "ThreeSpline.h"

ThreeSpline::ThreeSpline(void)
{
} ThreeSpline::~ThreeSpline(void)
{
}
//================================================================
// 函数功能: 利用求出的二阶导数求给定点值(结合Spline1,Spline2)
// 输入参数: *xa 为横坐标值,ya为纵坐标值,n为点个数,m为二阶偏导数
// x为给定点,y接收插出来的值
// 返回值: 无返回值
//================================================================
void ThreeSpline::Splint(double *xa,double *ya,double *m,int n,double &x,double &y)
{
int klo,khi,k;
klo=; khi=n-;
double hh,bb,aa; while(khi-klo>) // 二分法查找x所在区间段
{
k=(khi+klo)>>;
if(xa[k]>x) khi=k;
else klo=k;
}
hh=xa[khi]-xa[klo]; aa=(xa[khi]-x)/hh;
bb=(x-xa[klo])/hh; y=aa*ya[klo]+bb*ya[khi]+((aa*aa*aa-aa)*m[klo]+(bb*bb*bb-bb)*m[khi])*hh*hh/6.0;
}
//===========================================================================
// 函数功能: 对一系列点求二阶偏导数,点横坐标单调递增(I型边界)(结合Spline)
// 输入参数: *xa 为横坐标值,ya为纵坐标值,n为点个数,m为二阶偏导数(输出值)
// bound1、bound2为边界点一阶偏导数
// 返回值: 无返回值
//
//===========================================================================
void ThreeSpline::Spline1(double *xa,double *ya,int n,double *&m,double bound1,double bound2)
{
// 追赶法解方程求二阶偏导数
double f1=bound1,f2=bound2; double *a=new double[n]; // a:稀疏矩阵最下边一串数
double *b=new double[n]; // b:稀疏矩阵最中间一串数
double *c=new double[n]; // c:稀疏矩阵最上边一串数
double *d=new double[n]; double *f=new double[n]; double *bt=new double[n];
double *gm=new double[n]; double *h=new double[n];
m=new double[n]; for(int i=;i<n;i++) b[i]=; // 中间一串数为2
for(int i=;i<n-;i++) h[i]=xa[i+]-xa[i]; // 各段步长
for(int i=;i<n-;i++) a[i]=h[i-]/(h[i-]+h[i]);
a[n-]=; c[]=;
for(int i=;i<n-;i++) c[i]=h[i]/(h[i-]+h[i]); for(int i=;i<n-;i++)
f[i]=(ya[i+]-ya[i])/(xa[i+]-xa[i]); d[]=*(f[]-f1)/h[];
d[n-]=*(f2-f[n-])/h[n-]; for(int i=;i<n-;i++) d[i]=*(f[i]-f[i-])/(h[i-]+h[i]); bt[]=c[]/b[]; // 追赶法求解方程
for(int i=;i<n-;i++) bt[i]=c[i]/(b[i]-a[i]*bt[i-]); gm[]=d[]/b[];
for(int i=;i<=n-;i++) gm[i]=(d[i]-a[i]*gm[i-])/(b[i]-a[i]*bt[i-]); m[n-]=gm[n-];
for(int i=n-;i>=;i--) m[i]=gm[i]-bt[i]*m[i+]; delete a;
delete b;
delete c;
delete d;
delete gm;
delete bt;
delete f;
delete h;
}
//===========================================================================
// 函数功能: 对一系列点求二阶偏导数,点横坐标单调递增(II型边界)(结合Spline)
// 输入参数: *xa 为横坐标值,ya为纵坐标值,n为点个数,m为二阶偏导数(输出值)
// bound1、bound2为边界点二阶偏导数,当bound1和bound2不给值时则使用
// 默认值0,即自然边界
// 返回值: 无返回值
//
// 作者: 蒋锦朋 1034378054@qq.com
// 单位: 中国地质大学(武汉) 地球物理与空间信息学院
// 日期: 2014/12/03
//===========================================================================
void ThreeSpline::Spline2(double *xa,double *ya,int n,double *&m,double bound1,double bound2)
{
// 追赶法解方程求二阶偏导数
double f11=bound1,f22=bound2; double *a=new double[n]; // a:稀疏矩阵最下边一串数
double *b=new double[n]; // b:稀疏矩阵最中间一串数
double *c=new double[n]; // c:稀疏矩阵最上边一串数
double *d=new double[n]; double *f=new double[n]; double *bt=new double[n];
double *gm=new double[n]; double *h=new double[n];
m=new double[n]; for(int i=;i<n;i++) b[i]=;
for(int i=;i<n-;i++) h[i]=xa[i+]-xa[i];
for(int i=;i<n-;i++) a[i]=h[i-]/(h[i-]+h[i]);
a[n-]=; c[]=;
for(int i=;i<n-;i++) c[i]=h[i]/(h[i-]+h[i]); for(int i=;i<n-;i++)
f[i]=(ya[i+]-ya[i])/(xa[i+]-xa[i]); //d[0]=6*(f[0]-f1)/h[0];
//d[n-1]=6*(f2-f[n-2])/h[n-2]; for(int i=;i<n-;i++) d[i]=*(f[i]-f[i-])/(h[i-]+h[i]); d[]=d[]-a[]*f11;
d[n-]=d[n-]-c[n-]*f22;
//bt[0]=c[0]/b[0];
//for(int i=1;i<n-1;i++) bt[i]=c[i]/(b[i]-a[i]*bt[i-1]); //gm[0]=d[0]/b[0];
//for(int i=1;i<=n-1;i++) gm[i]=(d[i]-a[i]*gm[i-1])/(b[i]-a[i]*bt[i-1]); //m[n-1]=gm[n-1];
//for(int i=n-2;i>=0;i--) m[i]=gm[i]-bt[i]*m[i+1];
// 追赶法求解方程
bt[]=c[]/b[];
for(int i=;i<n-;i++) bt[i]=c[i]/(b[i]-a[i]*bt[i-]); gm[]=d[]/b[];
for(int i=;i<=n-;i++) gm[i]=(d[i]-a[i]*gm[i-])/(b[i]-a[i]*bt[i-]); m[n-]=gm[n-];
for(int i=n-;i>=;i--) m[i]=gm[i]-bt[i]*m[i+]; m[]=f11;
m[n-]=f22; delete a;
delete b;
delete c;
delete d;
delete gm;
delete bt;
delete f;
delete h;
}

调用方法

void CMathTestView::OnSpline()
{
// TODO: 在此添加命令处理程序代码
double x[]={0.52,,17.95,28.65,50.65,104.6,156.6,260.7,364.4,,,};
double y[]={5.28794,13.84,20.2,24.9,31.1,36.5,36.6,,20.9,7.8,1.5,0.2}; double xd[]={,,,,,,,};
double f1=1.86548,f2=-0.046115;
double f11=-0.279319,f22=0.0111560;
CSpline spline;
double *m;
spline.Spline2(x,y,,m,f11,f22); CString showstr;
showstr.Empty();
for(int i=;i<;i++)
{
double yd;
spline.Splint(x,y,m,,xd[i],yd);
CString s;
s.Format(_T("%lf %lf\n"),xd[i],yd);
showstr+=s;
}
delete m;
MessageBox(showstr);
}

仅仅记录测方法的使用过程

 在使用过程中,要对横坐标进行从小到大的排序,还有去除重复点,否则结果出错

05-11 09:41
查看更多