断层面拟合(Matlab)
一、算法目的
断层面的可视化对于地质构造分析、油气勘探开发具有重要作用。常用的断层面重构方法是利用断层数据进行曲面拟合，然而的断层数据一般用离散不规则分布的数据散点云来反映，需要基于这些散点云数据进行三维地质建模。建模过程需要利用插值算法与计算机图形学以获得光滑、连续的断层面。
算法将实现散点云\(\Longrightarrow\)断层面的过程。
二、算法流程
三、程序结构
四、程序分析
0 读取断层数据
待处理的数据格式如下表所示，为文本格式，程序首先从origin_point.txt中获取断层散点
| 大地坐标X | 大地坐标Y | 时间T | 线号 | 道号 | 所属断层线索引 |
|---|---|---|---|---|---|
| 598449.6 | 5132452 | 1473.383 | 786 | 1682 | 0 |
| 598470 | 5132452 | 1535.572 | 786 | 1685 | 0 |
| 598489.1 | 5132452 | 1560.448 | 786 | 1686 | 0 |
| 598619.9 | 5132532 | 1455.97 | 790 | 1699 | 1 |
| 598655.4 | 5132532 | 1516.916 | 790 | 1703 | 1 |
| 598698.9 | 5132532 | 1569.154 | 790 | 1707 | 1 |
| 598719.4 | 5132612 | 1454.726 | 794 | 1709 | 2 |
| 598749.4 | 5132612 | 1526.866 | 794 | 1712 | 2 |
faultorg=load("origin_point.txt");
X=faultorg(:,1);
Y=faultorg(:,2);
Z=-faultorg(:,3);
LineIndex=faultorg(:,4);
CdpIndex=faultorg(:,5);
Index=faultorg(:,6);
1 归一化、中心化
归一化是为了平衡各维度数据的量纲，减小方差大的特征的影响，并加快算法计算速度；
本次采用的归一化方法为线性归一化(Min-Max scaling)，将原始数据转到\([0,1]\)的范围内，公式为：
中心化的操作则是对各列数据减去列平均值，目的是将坐标原点放在数据的中心。
[Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin]=normfaultdata(X,Y,Z);
normfaultdata.m
function [Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin]=normfaultdata(X,Y,Z)
%{
*函数功能：* ***完成断层数据中心化、归一化***
*输入：*
- X：断层面的X坐标
- Y：断层面的Y坐标
- Z：断层面的Z坐标
*输出：*
- Xnorm/Ynorm/Znorm：空间坐标的归一化值
- Xmean/Ymean/Zmean：空间坐标的平均值
- Xmax/Ymax/Zmax：空间坐标的最大值
- Xmin/Ymin/Zmin：空间坐标的最小值
%}
Xmax=max(X);
Xmin=min(X);
Ymax=max(Y);
Ymin=min(Y);
Zmax=max(Z);
Zmin=min(Z);
Xnorm1=(X-Xmin)/(Xmax-Xmin); %归一化
Ynorm1=(Y-Ymin)/(Ymax-Ymin);
Znorm1=(Z-Zmin)/(Zmax-Zmin);
Xmean=mean(Xnorm1); % 数据中心化
Ymean=mean(Ynorm1);
Zmean=mean(Znorm1);
Xnorm=Xnorm1-Xmean;
Ynorm=Ynorm1-Ymean;
Znorm=Znorm1-Zmean;
end
2 拟合投影平面
投影的目的是将三维空间降低一个维度，从而简化空间中的插值问题，投影的实现将采用主成分分析(PCA)算法。
首先简要介绍PCA算法。
考虑这样一种情况
相当于把二维的信息降到一维，所以PCA的目的就是找这样一个坐标系，使得数据降维后信息损失最小，即信息保留最多。
假设我们手上的数据已服从正态分布，在经过中心化、归一化之后，那它可以由白数据经过拉伸、旋转变换得到。
可以看出求解\(\boldsymbol R\)矩阵是获取新坐标系的关键。
概括PCA的求解过程：
(1)协方差矩阵
协方差表示的是两个向量在变化过程中是同方向变化还是反方向变化，同向或反向的程度如何，协方差为正时，\(x\)，\(y\)为正相关。而协方差矩阵表示为：
白数据中\(x\)，\(y\)不相关，其协方差矩阵为\(C_0=\left[ {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right]\)，于是我们手上的数据的协方差：
又有拉伸变换\(S^T=S\)，旋转变换\(R^T=R^{-1}\)，令\(L=SS^T=\left[ {\begin{array}{*{20}{c}} a^2&0\\ 0&b^2 \end{array}} \right]\)，于是原式化为
(2)协方差矩阵的特征值
\(L=\left[ {\begin{array}{*{20}{c}} \lambda_1&0\\ 0&\lambda_2 \end{array}} \right]=SS^T=\left[ {\begin{array}{*{20}{c}} a^2&0\\ 0&b^2 \end{array}} \right]\)，\(a^2\)和\(b^2\)是两个轴方向的方差（拉伸倍数的平方），同时又是协方差矩阵的特征值。
(3)协方差矩阵的特征向量
\(R\)矩阵由\({\vec v}_1\)、\({\vec v}_2\)组成，\({\vec v}_1\)为新坐标系的\(X\)轴方向，\({\vec v}_2\)为新坐标系的\(Y\)轴方向，或者说这两个向量就是PCA的主成分。
(4)PCA与空间变换
从上文的介绍可以看出，PCA并未对数据做任何改变，它仅切换了坐标轴。其主成分（特征向量）则恰好表示新的基向量，例如从一开始的\(\left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&1&0\\ 0&0&1 \end{array}} \right]\)变为\(\left[ {\begin{array}{*{20}{c}} { - 0.63}&{0.43}&{0.64}\\ {0.35}&{0.90}&{ - 0.25}\\ { - 0.68}&{0.07}&{ - 0.72} \end{array}} \right]\)，则表示X轴由通常理解的方向变为向量\([-0.63,0.35,-0.68]\)指向的方向。但这一过程仅变动了坐标轴指向，并未发生线性变换，即并未带动数据一起转动。
当仅讨论三维数据时，PCA的第一、第二主成分即为方差最大（\(X\)轴方向）、次大（\(Y\)轴方向）的方向，而最小特征值对应的特征向量则是\(Z\)轴，这个新的\(Z\)轴是降维平面\(X-Y\)的法向量。设其为\(v_1=[a,b,c]\)，中心点(Xmean,Ymean,Zmean)过平面，易得截距\(d\)，于是PCA投影所成的平面可表示为：
本节中将用到PCA拟合投影面和获取投影坐标点两个函数：
[a,b,c,d]=planeSurfaceFitting(Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean);
planeSurfaceFitting.m
function [a,b,c,d,V2]=planeSurfaceFitting(Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean)
%{
*函数功能：* ***断层面数据拟合断层平面***
*输入：*
- Xnorm/Ynorm/Znorm：空间坐标的归一化值
- Xmean/Ymean/Zmean：空间坐标的平均值
*输出：*
- a/b/c/d：断层面多项式系数
- V2：得到的特征向量
%}
C=[Xnorm,Ynorm,Znorm]; % 输入的矩阵为归一化后的空间坐标
Cov=C'*C; % CT*C得到方阵
[eigenV,eigenU]=eig(Cov); % 求出特征值
[~,IndexValue]=min(diag(eigenU)); %最小的特征值对应的向量即为法向量
[~,Index]=sort(diag(eigenU),'descend');
V2=eigenV(:,Index); %V2中存放了所有的特征向量
V1=eigenV(:,IndexValue);
a=V1(1);
b=V1(2);
c=V1(3);
d=-dot([a b c],[Xmean,Ymean,Zmean]); %均值点过投影平面，算出截距
end
接下来根据原坐标点和投影平面，获取投影点坐标。
[x_projection,y_projection,z_projection]=projectplane(Xnorm,Ynorm,Znorm,[a,b,c,d]);
projectplane.m
function [x_projection,y_projection,z_projection]=projectplane(X,Y,Z,cfitt)
%{
*函数功能 * ***获取投影点坐标***
*输入：*
- X/Y/Z：原始空间坐标
- cfitt：断层面多项式系数列表
*输出：*
- x_projection/y_projection/z_projection：拟合面上投影点的空间坐标
%}
x_projection=zeros(1,length(X));
y_projection=zeros(1,length(Y));
z_projection=zeros(1,length(Z));
A0=cfitt(1);
B0=cfitt(2);
C0=cfitt(3);
D0=cfitt(4);
for i=1:length(X)
temp=-(A0*X(i)+B0*Y(i)+C0*Z(i)+D0)/(C0^2+A0^2+B0^2); %点到平面的距离公式
x_projection(i)=X(i)+A0*temp; %投影之后的新的x,y,z
y_projection(i)=Y(i)+B0*temp;
z_projection(i)=Z(i)+C0*temp;
end
end
随后汇总和调整新坐标系下的数据
data=[x_projection;y_projection;z_projection]; %新坐标存入data中
ori = (data(:,1)+data(:,end))/2; %取始末点的坐标中点为新的坐标原点
xorder=data(:,1)-data(:,end);%选取原始始末点相减为投影后x轴方向
xor=xorder/norm(xorder);%2维坐标下的x轴方向向量
3 旋转投影平面
矩阵与线性变换
回到线性变换与矩阵的本质，一个二维线性变换可以仅由四个数字完全确定，变换后的\(i\)的两个坐标和变换后的\(j\)的两个坐标，将这两个坐标封装在2×2的矩阵中。现在有一个描述线性变换的2×2矩阵\(\left[ {\begin{array}{*{20}{c}} a&b\\ c&d \end{array}} \right]\)，以及任意初始向量\(\vec v=(x,y)\)，我们将很容易得到变换之后的新向量，并通过矩阵乘法描述
这时的\(\left[ {\begin{array}{*{20}{c}} a&b\\ c&d \end{array}} \right]\)就可以看作一种函数，其逆矩阵\(\left[ {\begin{array}{*{20}{c}} a&b\\ c&d \end{array}} \right]^{-1}\)则相当于“重置”操作，这一思想将在后续程序中加以运用。
基向量的改变引起空间变换，而矩阵记录了变换之后基向量的落点。在上一步PCA中得到的基向量为\(\left[ {\begin{array}{*{20}{c}} { - 0.63}&{0.43}&{0.64}\\ {0.35}&{0.90}&{ - 0.25}\\ { - 0.68}&{0.07}&{ - 0.72} \end{array}} \right]\)，而接下来要执行的旋转投影平面，则需要将这一基向量“校正”至“初始状态”，即使用拉伸(\(S\))、旋转(\(R\))变换将基向量\(\vec i\)重新指向\([1,0,0]\)，基向量\(\vec j\)重新指向\([0,1,0]\)，基向量\(\vec k\)重新指向\([0,0,1]\)。反投影则只需对变换矩阵\(RS\)求逆，反变换即可。
将三维空间中的投影平面通过旋转、平移至水平
[data_xyzor,R1]=Convert3Dto2D([a,b,c,d],xor,ori,norm_data);
function [data_xyzor,R1]=Convert3Dto2D(cfitt,xor,ori,data)
%{
*函数功能：*将三维空间中的X-Y平面旋转、平移至水平面
*输入：*
- cfitt：断层面多项式系数列表
- xor：投影面 x轴方向向量
- ori：投影面原点坐标
- data：待旋转、投影的数据/散点
*输出：*
- data_xyzor：旋转、平移后的散点
- R1：变换矩阵（旋转*平移）
%}
aall=cfitt(1);ball=cfitt(2); call=cfitt(3); %从输入列表中提取出a,b,c
zor=[aall,ball,call]/norm([aall,ball,call]); %汇总到法向量z，并标准化z方向
yor=cross(xor,zor);yor=yor/norm(yor); %现在有了x方向和z方向的方向向量，用叉乘算出y方向向量
move=[1,0,0,0;0,1,0,0;0,0,1,0;-ori',1]; %平移矩阵，以新的坐标原点去描述平移的尺度
%旋转矩阵，以xor为例，其为新的空间中的一个坐标轴，上面已经标准化为单位长度，即xor只剩下旋转的分量，表示从最初的x轴--->xor所需要的变换向量
rotate=[xor(1),yor(1),zor(1),0;xor(2),yor(2),zor(2),0;xor(3),yor(3),zor(3),0;0,0,0,1];
R1=move*rotate; %R1为整个变换矩阵
data(4,:)=1; % 将数据补充至四维，增加一个平移维度
data_xyzor=(data')*R1;% 乘变换矩阵，所有散点跟随坐标系旋转到X、Y轴水平的状态
end
4 限制插值区域
这一步主要是为了裁剪掉后续多余的插值点，使用Matlab中的凸包convhull工具，根据现有的基础散点，作出其凸包络线，对插值范围进行约束。
Xr=data_xyzor(:,1); %获得数据点新的x,y,z坐标
Yr=data_xyzor(:,2);
Zr=data_xyzor(:,3);
[k,a1]=convhull(Xr,Yr);
xv=Xr(k);
yv=Yr(k);
此时已经获得二维平面的散点，\(T\)值将是待插属性。
经此步骤，得到断层散点的边界线，如图；
5 二维平面插值
现在手上的数据变成了(Xr，Yr)以及对应坐标上的属性值Tr，目标是在圈定的范围内插出新的、过渡的属性值。这一环节已简化为二维的插值任务。
本次插值采用的是双调和样条插值(biharmonic spline interpolation)算法（matlab中的V4算法），是不规则间隔的二维数据点的插值，简要介绍该算法：
插值平面是以每个数据点为中心的格林函数的线性组合。 通过求解线性方程组可以找到格林函数的幅度。
插值平面\(S(\boldsymbol{x})\)的表达式为：
其中\(\boldsymbol x\)是待插值的样本点坐标，\(g(\boldsymbol{x},\boldsymbol{x}_{j})\)为格林函数，\(\boldsymbol{x}_{i}\)、\(\boldsymbol{x}_{j}\)分别是第个\(i、j\)个已知数据，\(n\)是平面的总已知点数，\(\boldsymbol{w}_{j}\)是\(\boldsymbol{x}与\boldsymbol{x}_{j}\)的相关性的权重矩阵。
算法步骤如下：
①假设插值平面所需的输出节点是\(p_0\)，最接近\(p_0\)的插值点有\(k\)个，\(p_0\)到第\(j\)个已知点的距离按下式计算：
其中\(\boldsymbol{x}_{0}\)是\(p_0\)的位置向量；所有距离都使用合并排序算法按升序排序，在排序之后，获取前\(k\)个数据点的索引和坐标；
②计算这\(k\)个数据点的格林函数矩阵\(\boldsymbol G\)：
根据格林函数矩阵和属性矩阵，权重矩阵可以计算：
③计算点\(p_0\)的\(1×k\)的格林函数矩阵$ \mathbf{G}_{p_0}$：
易得，点\(p_0\)的属性值\(z_p\)：
④重复步骤"①~③"，计算其他点\(p_i\)插值后的属性值\(z_{p_i}\).
程序将使用Matlab中的griddata函数中的v4算法，在平面上对\(T\)值进行插值。
Xrmax=max(Xr); %最大最小值坐标
Xrmin=min(Xr);
Yrmax=max(Yr);
Yrmin=min(Yr);
h=0.01; %计算步长
[Xin,Yin]=meshgrid(Xrmin:h:Xrmax,Yrmin:h:Yrmax); %搭建插值网格，假设X方向上的采样点为100个，Y方向上50个，那么Xin和Yin都将是100×50的矩阵，Xin与Yin组合，方能表示矩阵中任一位置的坐标
Zin=griddata(Xr,Yr,Zr,Xin,Yin,'v4'); %进行插值，Zin的大小将和Xin、Yin相同，Zin也可延申为在这个平面网格上的其他属性值
插值完成后，剔除范围外的插值点：
% 先将Xin,Yin,Zin展成一维，因为要配合裁剪的xv与yv也是一维的
[ny,nx]=size(Zin);
Xin_row1=zeros(1,ny*nx);
Yin_row1=zeros(1,ny*nx);
Zin_row1=zeros(1,ny*nx);
for iy=1:ny
for ix=1:nx
index=(iy-1)*nx+ix;
Xin_row1(index)=Xin(iy,ix);
Yin_row1(index)=Yin(iy,ix);
Zin_row1(index)=Zin(iy,ix);
end
end
% 同样是Matlab工具箱中的函数
index_cut=inpolygon(Xin_row1,Yin_row1,xv,yv);
% 得到插值、裁剪后的一维数据
Xin_row=Xin_row1(index_cut);
Yin_row=Yin_row1(index_cut);
Zin_row=Zin_row1(index_cut);
插值后的效果为：
6 散点云三角剖分
插值后的数据间隔已经非常小，但依然是空间中的三维散点，需要将这些散点连接成面，才算达到构建断层面的效果。此时将用到点集的三角剖分技术(Triangulation)，三角剖分研究的是如何把散点集合剖分成不均匀的三角形网格。其中Delaunay三角剖分算法较为成熟，广泛应用与各领域图像和点云的处理上。
使用matlab中的库函数delaunay实现
%进行Delauny三角剖分算法
%输入为点坐标P，以向量/矩阵形式指定
tri=delaunay(Xin_row,Yin_row); % 输出tri的每行都包含P的行号（索引），它们构成三角剖分中的单个三角形
三角剖分后的平面如图，此时只需按照Z值进行渲染，便能够得到最终的插值断层面。
7 反投影
pr_in=[Xin_row',Yin_row',Zin_row',ones(length(Zin_row),1)]; %同样增加一维
Faultnorm_in=pr_in*inv(R1); %乘变换矩阵的逆，返回原坐标系
8 反归一化
Xin_norm=Faultnorm_in(:,1);
Yin_norm=Faultnorm_in(:,2);
Zin_norm=Faultnorm_in(:,3);
[Xorg,Yorg,Zorg]=inversenormfaultdata(Xin_norm,Yin_norm,Zin_norm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin);
function [Xorg,Yorg,Zorg]=inversenormfaultdata(Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin)
%{
*函数功能：*断层数据反归一化
*输入：*
- Xnorm/Ynorm/Znorm：断层面上的归一化空间坐标
- Xmean/Ymean/Zmean：断层面上空间坐标的平均值
- Xmax/Ymax/Zmax：断层面上空间坐标的最大值
- Xmin/Ymin/Zmin：断层面上空间坐标的最小值
*输出：*
- Xorg/Yorg/Zorg：反归一化后的空间坐标，即原始坐标系下的坐标值
%}
Xorg1=Xnorm+Xmean;
Yorg1=Ynorm+Ymean;
Zorg1=Znorm+Zmean;
Xorg=Xorg1*(Xmax-Xmin)+Xmin;
Yorg=Yorg1*(Ymax-Ymin)+Ymin;
Zorg=Zorg1*(Zmax-Zmin)+Zmin;
end
得到插值、拟合后的最终断层面
main函数完整代码
clc
clear all
close all
%% 导入断层散点
faultorg=load("origin_point.txt");
X=faultorg(:,1);
Y=faultorg(:,2);
Z=-faultorg(:,3);
LineIndex=faultorg(:,4);
CdpIndex=faultorg(:,5);
Index=faultorg(:,6);
%% 数据归一化、中心化
[Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin]=normfaultdata(X,Y,Z);
%% Step three: 断层面PCA算法拟合及投影
[a,b,c,d]=planeSurfaceFitting(Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean);
%坐标旋转投影
[x_projection,y_projection,z_projection]=projectplane(Xnorm,Ynorm,Znorm,[a,b,c,d]);
data=[x_projection;y_projection;z_projection];
ori=(data(:,1)+data(:,end))/2;%投影后的坐标原点，即转换为2维坐标后的原点
xorder=data(:,1)-data(:,end);%选取原始始末点相减为投影后x轴方向
xor=xorder/norm(xorder);%2维坐标下的x轴方向向量
norm_data = [Xnorm';Ynorm';Znorm']
[data_xyzor,R1]=Convert3Dto2D([a,b,c,d],xor,ori,norm_data);
Xr=data_xyzor(:,1);
Yr=data_xyzor(:,2);
Zr=data_xyzor(:,3);
%% 形成凸多边形
[k,a1]=convhull(Xr,Yr);
xv=Xr(k);
yv=Yr(k);
%% 插值
Xrmax=max(Xr);
Xrmin=min(Xr);
Yrmax=max(Yr);
Yrmin=min(Yr);
h=0.01; %计算步长
[Xin,Yin]=meshgrid(Xrmin:h:Xrmax,Yrmin:h:Yrmax); %搭建插值网格
Zin=griddata(Xr,Yr,Zr,Xin,Yin,'v4'); %进行插值
[ny,nx]=size(Zin);
Xin_row1=zeros(1,ny*nx);
Yin_row1=zeros(1,ny*nx);
Zin_row1=zeros(1,ny*nx);
for iy=1:ny
for ix=1:nx
index=(iy-1)*nx+ix;
Xin_row1(index)=Xin(iy,ix);
Yin_row1(index)=Yin(iy,ix);
Zin_row1(index)=Zin(iy,ix);
end
end
%剔除部分点
index_cut=inpolygon(Xin_row1,Yin_row1,xv,yv);
Xin_row=Xin_row1(index_cut);
Yin_row=Yin_row1(index_cut);
Zin_row=Zin_row1(index_cut);
%进行Delauny三角剖分算法
tri=delaunay(Xin_row,Yin_row);
%反向旋转投影变化
pr_in=[Xin_row',Yin_row',Zin_row',ones(length(Zin_row),1)];
Faultnorm_in=pr_in*inv(R1);
%反归一化操作恢复原始坐标
Xin_norm=Faultnorm_in(:,1);
Yin_norm=Faultnorm_in(:,2);
Zin_norm=Faultnorm_in(:,3);
[Xorg,Yorg,Zorg]=inversenormfaultdata(Xin_norm,Yin_norm,Zin_norm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin);