空间曲线的线性参数插值

在断层曲面拟合的过程中，发现当解释的空间数据点过于稀疏的化，其断层面拟合的效果较差，我们采用空间曲线线性插值加密的算法，增加插值控制点的数量，改善插值的效果。

1.1 问题描述即算法描述

已知空间三维离散折线 \(l=(p_1,p_2,...,p_i,...,p_n)\) 上有 \(n\)个点 \(p_i(x_i,y_i,z_i)\)，假设点与点之间是线性渐变的，我们将在折线段中两点中添加\(m\)个点,即 \({...,p_{i},q_{1}^{i},q_2^{i},q_{3}^i,...,q_{j}^{i},...,q_{m}^{i},p_{i+1},...}\) , 其中\(q_{j}^i=(x_{j}^i,y_{j}^i,z_j^i)\) 为待插值点，其中引入参数\(t_j\), 其待插值点的坐标如下：

\[\begin{equation} \begin{aligned} x_j^{i}=(1-t_j)x_i+t_jx_{i+1} \\ y_j^{i}=(1-t_j)y_i+t_jx_{i+1} \\ z_j^{i}=(1-t_j)z_i+t_jz_{i+1} \\ \end{aligned} \tag{1.1.1} \space (i=1,2,3,...,n-1) and(j=1,2,3,..,m) \end{equation} \]

其中插值参数\(t_j\)：

\[\begin{aligned} \begin{split} t_j &=j\Delta{t}\\ \Delta{t}&=1/(m+1)\\ j&=1,2,3,...,m \end{split} \end{aligned} \tag{1.2.1} \]

1.2 算法代码

1.断层插值算法：

%函数的功能：绘制断层线的解释网格
%函数的描述：利用线性插值断层线解释控制点加密
%函数的使用：y=func(input1,input2)
%输入：
%      X-待插值的断层解释控制点 X坐标向量 （n）
%      Y-待插值的断层解释控制点 Y坐标向量 （n） 
%      T-待插值的断层解释控制点 T(或Z)坐标向量 （n）
% CDP-待插值的断层解释控制点 CDPIndex 向量 (n)
% LINE-待插值的断层解释控制点 LINEIndex向量 （n）
% faultIndex-待插值的解释控制点 断层编号向量 （n）
%       n-每个折线加密点数
%输出：
%      X_inter-线性插值后的断层解释控制点X坐标向量
%      Y_inter-线性插值后的断层解释控制点Y坐标向量
%      T_inter -线性插值后的断层解释控制点T坐标向量
%      CDP_inter-线性插值后的断层的插值控制点CDPIndex向量
%     LINE_inter-线性插值后的断层的插值控制点LINEIndex向量
%     faultIndex_inter-线性插值后的断层的插值控制点断层编号向量 
%注意事项：利用函数的适用范围。
%文档日期：
%标签：
%创建日期：
%最后更新日期：
function [X_inter,Y_inter,T_inter,CDP_inter,LINE_inter,faultIndex_inter]=interfaultInline(X,Y,T,CDP,LINE,faultIndex,n)

faultIndexArray=unique(faultIndex);
nfaultline=length(faultIndexArray);

X_inter=[];
Y_inter=[];
T_inter=[];
CDP_inter=[];
LINE_inter=[];
faultIndex_inter=[];

t=linspace(0,1,n+2);
for ifault=1:nfaultline
      faltIndexTemp=faultIndexArray(ifault);
      indexSubTemp=find(faultIndex==faltIndexTemp);
       X_temp=X(indexSubTemp);
       Y_temp=Y(indexSubTemp);
       T_temp=T(indexSubTemp);
       CDP_temp=CDP(indexSubTemp);
       LINE_temp=LINE(indexSubTemp);

       nSample=length(indexSubTemp);
       nInterSample=(nSample-1)*n+nSample

       Xinter_temp=zeros(nInterSample,1);
       Yinter_temp=zeros(nInterSample,1);
       Tinter_temp=zeros(nInterSample,1);
       CDPinter_temp=zeros(nInterSample,1);
       LINEinter_temp=zeros(nInterSample,1);
       faultIndexinter_temp=zeros(nInterSample,1);

       for isample=1:nSample-1
          for iter=1:n+2
              Xinter_temp((isample-1)*(n+1)+iter)=(1-t(iter))*X_temp(isample)+t(iter)*X_temp(isample+1);
              Yinter_temp((isample-1)*(n+1)+iter)=(1-t(iter))*Y_temp(isample)+t(iter)*Y_temp(isample+1);
              Tinter_temp((isample-1)*(n+1)+iter)=(1-t(iter))*T_temp(isample)+t(iter)*T_temp(isample+1);
              CDPinter_temp((isample-1)*(n+1)+iter)=(1-t(iter))*CDP_temp(isample)+t(iter)*CDP_temp(isample+1);
              LINEinter_temp((isample-1)*(n+1)+iter)=(1-t(iter))*LINE_temp(isample)+t(iter)*LINE_temp(isample+1);
              faultIndexinter_temp((isample-1)*(n+1)+iter)=ifault;
          end
       end
     
       X_inter=[X_inter;Xinter_temp];
       Y_inter=[Y_inter;Yinter_temp];
       T_inter=[T_inter;Tinter_temp];
       CDP_inter=[CDP_inter;CDPinter_temp];
       LINE_inter=[LINE_inter;LINEinter_temp];
       faultIndex_inter=[faultIndex_inter;faultIndexinter_temp];
end



end

2.断层解释线绘制代码

function drawfaultline(X,Y,T,Index,method)
%函数的功能：绘制断层线的解释网格
%函数的描述：
%函数的使用：y=func(input1,input2)
%输入：
%     input1:
%     input2:
%输出：
%     Y:
%例子：y=func(1,'type1');
%注意事项：利用函数的适用范围。
%文档日期：
%标签：
%创建日期：
%最后更新日期：
maxIndex=max(Index);
minIndex=min(Index);
k=1;
for index=minIndex:maxIndex
     indexfalut=find(Index==index);
     faultlineX=X(indexfalut);
     faultlineY=Y(indexfalut);
     faultlineT=T(indexfalut);
     k=k+1;
     plot3(faultlineX,faultlineY,faultlineT,method);
     hold on
end
set(gca,'Zdir','reverse');
daspect([1,1,1]);
hold off
end

1.3 算法效果及案例

断层线插值案例

clc
clear all
close all

%% 读取断层面散点数据
pointdata=load("./data/F10_GeoInline.txt");
pointdata1=load("./data/F10_GeoXline.txt");

maxIndexLine=max(pointdata(:,6));

%% 断层线数据加密
Xpoint=[pointdata(:,1);pointdata1(:,1)];
Ypoint=[pointdata(:,2);pointdata1(:,2)];
Tpoint=[pointdata(:,3);pointdata1(:,3)];
LINEpoint=[pointdata(:,4);pointdata1(:,4)];
CDPpoint=[pointdata(:,5);pointdata1(:,5)];
INDEXpoint=[pointdata(:,6);pointdata1(:,6)+1+maxIndexLine];
n=20;
[X_inter,Y_inter,T_inter,CDP_inter,LINE_inter,faultIndex_inter]=interfaultInline(Xpoint,Ypoint,Tpoint,CDPpoint,LINEpoint,INDEXpoint,n);

Xpoint1=X_inter;
Ypoint1=Y_inter;
Tpoint1=T_inter;
LINEpoint1=LINE_inter;
CDPpoint1=CDP_inter;
INDEXpoint1=faultIndex_inter;

%%
figure(1)
subplot(1,2,1)
drawfaultline(Xpoint,Ypoint,Tpoint,INDEXpoint,'-or');
xlabel("X");
ylabel("Y");
zlabel("Z");
title("插值前断层解释控制点")
grid on
subplot(1,2,2)
drawfaultline(Xpoint,Ypoint,Tpoint,INDEXpoint,'or');
hold on
drawfaultline(Xpoint1,Ypoint1,Tpoint1,INDEXpoint1,'.b');
xlabel("X");
ylabel("Y");
zlabel("Z");
title("插值后断层解释控制点")
grid on

断层插值后的效果

posted @ 2024-12-29 15:19 GeoFXR 阅读(128) 评论(0) 收藏举报

刷新页面返回顶部

GeoFXR

空间曲线的线性参数插值