Surfe’s Inputs and Outputs 输入和输出
Definitions 定义
Interface: a boundary surface separating two different regions. E.g. interface that separates two lithologies, a fault surface that separates two continuous domains/volumes.
地层界面:分隔两个不同区域的边界面。例如,分隔两个不同岩性的界面,断层面分隔两个连续域/体。
CPD function: Conditionally Positive Definite function. Conditional in the sense that the interpolation matrix is positive definite subject to the orthogonality constraints provided by the polynomial in the null space.
条件正定函数:条件正定函数。插值矩阵是正定的,受零空间中多项式提供的正交性约束的条件。
SPD function: Strictly Positive Definite function. A Gramm matrix (e.g. our interpolation matrices) generated by such a function/kernel will have all of items eigenvalues positive.
严格正定函数:严格正定函数。由这种函数/核生成的Gramm矩阵(例如我们的插值矩阵)的所有项特征值都是正的。
All input parameters and constraints as well as outputs to the Surfe interpolant are made through the Surfe_API class.
所有输入参数和约束以及Surfe 插值的输出都是通过Surfe_API类生成的。
Input Parameters
Modelling method
The modelling method parameter initializes the Surfe_API object.
使用建模方法参数来初始化Surfe_API
对象。
Set by
Surfe_API(const int &modelling_method)
Usage example:
Surfe_API surfe(1); // Single Surface method
Surfe_API *surfe = new Surfe_API(2); // Lajaunie method
Acceptable integer values:
-
1 : Single_surface
- Models a single interface surface. No increment points are used in this method and therefore is faster to obtain the interpolant and evaluate it at a list of points.
- 建立单一界面模型。此方法不使用增量点,因此获得插值点并在(重采样)点集中(重采样)对其求值的速度更快。
-
2 : Lajaunie_approach
- Models multiple conformal interface surfaces. Increment points are used in this method.
- 建立多个共形界面模型。这种方法使用增量点。
-
3 : Vector_field
- Models vectors fields from gradient and tangent constraints. Does not incorporate interface or increment constraints.
- 从梯度和切线约束中建模向量场模型。不包含界面或增量约束。
-
4: Stratigraphic_horizons
- Models multiple conformal interface surfaces with additional constraints on the order in which the layers/interfaces where deposited/created. E.g. Layer A is on top/younger of Layer B, Layer B is on top/younger of Layer C. This method is very useful in situations where the data does not sample the volume optimally. E.g. Interface data obtain from outcrops (horizon sampling bias) – Lajaunie method fails in these situation: The increment constraint only indicate that points on the same interface have the same scalar value – nothing about the order of these interfaces are imposed by such constraints.
- 对多个共形界面曲面建模,并对放置(排列)/创建层/界面层的顺序进行附加约束。例如A层位于B层的顶部(较年轻层),B层位于C层的顶部(较年轻层)。这种方法在数据没有最佳采样体的情况下非常有用。例如从露头获得的界面数据(层位采样偏差)- Lajaunie方法在这种情况下失败:增量约束仅表示同一界面上的点具有相同的标量值-这些界面的顺序不受此类约束的影响。
-
5: Continuous_property
- Models scalar data like assay data, points sampling a function like a fourier series. Does not incorporate gradient constraints.
- 标量数据模型,如分析数据,点采样函数,如傅立叶级数。不包含渐变约束。
RBF Kernel
The type of RBF/Kernel used in the interpolant.
Note: If CPD kernels are used special care must be taken when using inequality constraints. This is because a convex optimization problem must be solved in this case. As such, these methods require a SPD matrices. To make our matrices PD we use Lagrangian polynomial basis to ensure our functional space resides within a Reproducing Hilbert Kernel Space. Currently, the only CPD function that is supported is the cubic in these cases.
用于插值的径向基函数(核)的类型。
注意:如果使用CPD
(条件正定函数)内核,在使用不等式约束时必须特别小心。这是因为在这种情况下必须解决凸优化问题。因此,这些方法需要SPD
(严格正定函数)矩阵。为了使我们的矩阵PD
(正定),我们使用拉格朗日多项式基来确保我们的函数空间位于Reproducing Hilbert Kernel Space
内。目前,在这些情况下,唯一受支持的CPD
函数是cubic
。
Acceptable strings:
-
“r3” :
(CPD)(
cubic
三次立方体) -
“Gaussian” :
(SPD)(高斯函数)
-
“Multiquadrics”:
(CPD)
-
“Inverse Multiquadrics”:
(SPD)
-
“Thin Plate Spline”:
(CPD)
-
“r”:
(CPD)
-
“WendlandC2”:
(SPD)
-
“MaternC4”:
(SPD)
Important: For SPD kernels, a shape parameter is also needed. This is set by
对于SPD
内核,还需要一个形状参数。
surfe.SetRBFShapeParameter(const double &shape_parameter);
Polynomial order
Polynomials of the correct minimum order are required for CPD kernels to ensure an unique minimum norm interpolant. For cubic kernels the minimum order is 1. These polynomials are the drift functions in co-kriging.
CPD
核需要正确的最小阶多项式,以确保最小范数插值的唯一性。对于立方(cube)核,最小顺序为1。这些多项式是协同-克里格中的漂移函数。
Acceptable integer values:
- 0 order:
- 1 order:
- 2 order:
Regression Smoothing
Regression smoothing creates an approximate interpolant instead of exact fitting. Exact fitting with real (non-synthetic data) will likely produce topological errors especially with noisy or highly varying data. If a threshold for fitting data constraints is given effectively a least squares minimization is applied to residuals. In the co-kriging world this is the nugget effect.
回归平滑创建一个近似的插值而不是精确的拟合。与真实数据(非合成数据)的精确拟合可能会产生拓扑误差,特别是在噪声或高度变化的数据中。如果有效地给出了拟合数据约束的阈值,则对残差应用最小二乘法最小化。在协同-克里金中,这就是nugget effect
效应。
Set by
surfe.SetRegressionSmoothing(const bool &use_regression_smoothing, const double &amount);
Note: Depending on the kernel used and the structural complexity sampled by the data, the value specified does not necessarily correspond to physical interpretable meaning: e.g. distance
注意:根据使用的内核和数据采样的结构复杂性,指定的值不一定对应于物理可解释的含义:例如距离
Greedy Algorithm
Reduces the number of data constraints (a.k.a centers in approximation theory) used by the interpolant while also accurately approximating the interpolant. Effectively this is a massive speed boost in evaluation of the interpolant. This algorithm will start with the minimum number of constraints needed to obtain an interpolant, then the interpolant will be evaluated at every data constraint which was not included into the interpolant to measure residuals (how much discrepancy there is). A small number of residuals beyond the user specified residual threshold will then be added to the interpolant until all residuals are below this threshold. A threshold for interface data and orientation data are required.
在精确逼近插值点的同时,减少插值点使用的数据约束(a.k.a 中心近似理论)。实际上,这是一个巨大的速度提升在评估插值。此算法将从获得插值所需的最小约束数开始,然后在每个未包含在插值中的数据约束处对插值求值,以测量残差(存在多大差异)。然后,将超出用户指定的残差阈值的少量残差添加到插值中,直到所有残差都低于此阈值。界面数据和方向数据的阈值是必需的。
Set by
surfe.SetRegressionSmoothing(const bool &use_greedy, const double &interface_uncertainty, const double &angular_uncertainty);
Note: The disadvantage of this method is that a bias is placed on outliers. Also note, that support for this method is very limited. I believe only single surface method is support.
注:这种方法的缺点是对异常值有偏差。还要注意的是,这种方法的支持非常有限。我相信只有单一的表面方法是支持。
Interface residual at point not included in interpolant:
where f is the scalar field constraint for the interface point. For single surface methods, this is set to 0.
Orientation residual at point not included in interpolant:
Restricted Range
This option is somewhat like the greedy algorithm but does not suffer from the bias of outliers. All data constraints are included into the interpolant, and every data constraint has a bound on it where the smoothest interpolant is to be constrained within. To obtain such an interpolant a quadratic optimization problem is solved, which introduces computational overhead. But only finding the interpolant compute time is affected, the evaluation of the interpolant at user specified points is not affected.
约束范围
这个选项有点像贪婪算法,但不受异常值的影响。所有数据约束都包含在插值中,并且每个数据约束都有一个边界,其中最平滑的插值将被约束在其中。为了得到这样一个插值函数,需要求解一个二次优化问题,这会增加计算开销。但只影响插值计算时间的确定,而不影响插值在用户指定点的计算。
For interface points:
For orientation onts:
Set by
surfe.SetRestrictedRange(const bool &use_restricted_range, const double &interface_uncertainty, const double &angular_uncertainty);
Global Anisotropy
This option derives principle directions of anisotropy using all supplied planar/gradient constraints. These directions are used to modified how distances are computed. Works well if modelled structures are very global. However, works very poorly when plunge of structures largely vary.
此选项使用所有提供的平面/渐变约束导出各向异性的主方向。这些方向用于修改计算距离的方式。如果建模的结构是非常全局的,则可以很好地工作。然而,当结构的倾斜度变化很大时,效果很差。
Set by
surfe.SetGlobalAnisotropy(const bool &g_anisotropy);
Input Constraints
There are 4 types on data constraints that can be supplied into surfe:
有4种类型的数据约束可以提供给surfe
:
Interface constraints
- 3D points sampling an interface that separates two volumetric domains (e.g. lithologies, fault domains)
- (x, y, z, level)
- Level is a float value that organizes the stratigraphic order of the interfaces and lithologies (via inequalities). Important: Larger level values are younger (on top) than smaller level values (on bottom). Points sampling the same interface must have the same level value.
Level
(层位属性)是一个浮动值,用于组织界面和岩性的地层顺序(通过不等式)。较大的级别值(顶部)比较小的级别值(底部)年轻。采样相同界面的点必须具有相同的级别值。- Note: If only one interface is being model, any level value can be chosen. Normally 0 is specified in this case.
- 注意:如果只有一个界面模型,则可以选择任何级别值。在这种情况下,通常指定0。
- Added by
surfe.AddInterfaceConstraint(const double &x,const double &y,const double &z,const double &level);
- Added by
surfe.SetInterfaceConstraints(const MatrixXd &interface_constraints);
Accepted structure for MatrixXd:
Nx4 matrix, N = number of interface points
Columns: x, y, z, level
Note: In python, MatrixXd is a numpy array
Planar constraints
-
3D points that have normal (Younging) information attributed to them. The normal indicates the direction in which younger stratigraphy is found. If fault surfaces are being modelled, the polarity of the normal does not matter.
-
具有表面地层(年轻)信息的三维点。法线表示发现较年轻地层的方向。如果断层表面正在建模,则法向的极性无关紧要。
-
There are four ways planar constraints can be added:
- Supply normal vector
surfe.AddPlanarConstraintwNormal(const double &x,const double &y,const double &z,
const double &nx,const double &ny,const double &nz);
- Supply strike, dip, polarity
- Note polarity has two acceptable values: 0 == upright, 1 == overturned (points down)
surfe.AddPlanarConstraintwStrikeDipPolarity(const double &x,const double &y,const double &z,
const double &strike, const double &dip, const int &polarity);
- Supply azimuth, dip, polarity
- Note polarity has two acceptable values: 0 == upright, 1 == overturned (points down)
- Note that azimuth is the same as dip direction
surfe.AddPlanarConstraintwAzimuthDipPolarity(const double &x,const double &y,const double &z,
const double &azimuth,const double &dip,const int &polarity);
- Supply a matrix or array
surfe.SetPlanarConstraints(const MatrixXd &planar_constraints);
Accepted structure for MatrixXd:
Nx6 matrix, N = number of planar points
Columns: x, y, z, nx, ny, nz
Tangent constraints
- 3D points that have vector
attributed to them. The relationship between this vector and the scalar field is
. In other words, the vector
is orthogonal (90°) with respect to the gradient of the scalar field
. This constraint does not have a lot of effect on changing the modeled geometry since there is a large amount of freedom with fitting this constraint. However, if a lot of these are specified in addition to supplying two tangent constraints at the same point (different vectors!) then this constraint can be useful especially for foliation orientations (there is no polarity).
此约束对更改建模的几何图形没有太大影响,因为拟合此约束有很大的自由度。但是,如果除了在同一点上提供两个切线约束(不同的向量!)这样,这个约束特别适用于foliation orientations
(没有极性)。 - Added by
surfe.AddTangentConstraint(const double &x,const double &y,const double &z,
const double &tx,const double &ty,const double &tz);
- Added by
surfe.SetTangentConstraints(const MatrixXd &tangent_constraints);
Accepted structure for MatrixXd:
Nx6 matrix, N = number of tangent points
Columns: x, y, z, tx, ty, tz
Inequality constraints
- 3D points sampling lithologies or inside/outside points relative to an interface/s. E.g. Does not belong to an interface. Very useful constraint since observing interface between lithologies is extremely rare and this is the most abundant type of data especially for Geological surveys.
三维点采样岩性或相对于界面的内部/外部点。例如,不属于界面。非常有用的约束条件,因为观察岩性之间的界面非常罕见,这是最丰富的数据类型,特别是地质调查 - (x, y, z, level)
- Note: the level property has to be compatible and make sense with respect to the level property of the interface level property. For example, you can’t have an inequality point having the same level value as an interface point.注意:
level
属性必须与interface level
属性的leve
l属性兼容并且有意义。例如,不能有与界面点具有相同级别值的不等式点。
- Note: the level property has to be compatible and make sense with respect to the level property of the interface level property. For example, you can’t have an inequality point having the same level value as an interface point.注意:
- Added by
surfe.AddInequalityConstraint(const double &x,const double &y,const double &z,const double &level);
- Added by
surfe.SetInequalityConstraints(const MatrixXd &inequality_constraints);
Accepted structure for MatrixXd:
Nx4 matrix, N = number of inequality points
Columns: x, y, z, level
Outputs
Scalar field
Obtain the value of the scalar field at a 3D point
获取三维点处标量场的值
Get by
surfe.EvaluateInterpolantAtPoint(const double &x,const double &y,const double &z);
Returns a double
Gradient field
Obtain the gradient of the scalar fied at a 3D point
获取三维点上标量的梯度
Get by
surfe.EvaluateVectorInterpolantAtPoint(const double &x,const double &y,const double &z);
Returns a 3D vector
Spatial Parameters
Obtain spatial metrics of the inputted data constraints
Get by
surfe.GetDataBoundsAndResolution();
Returns a data structure called SpatialParameters
SpatialParameters data members: xmin, xmax, ymin, ymax, zmin, zmax, resolution
Note: resolution is the spatial resolution required to model the structural variability.
Interface Reference Points
For every interface being model there is a reference point for that interface (used for increments). Getting these reference points is critical when increments are used to model multiple conformal surfaces. This is because we don’t know what the scalar field values are associated to the interfaces. These values are need when applying marching cubes/tetrahedral algorithms to extract the iso surface from the modelled scalar field. The get these values you need to evaluate the interpolant at these reference points.
对于作为模型的每个界面,都有一个该界面的参考点(用于增量)。当增量用于建模多个共形曲面时,获取这些参考点至关重要。这是因为我们不知道标量字段值与界面关联的是什么。当应用marching cubes
(四面体算法)从建模标量场中提取等值曲面时,需要这些值。获取这些值需要在这些参考点计算插值。
Get by
surfe.GetInterfaceReferencePoints();
Returns an array/matrix
- Nx3 matrix, N = number of intefaces
Number of Interfaces
Gets the number of modeled interfaces by the interpolant.
Get by
surfe.GetNumberOfInterfaces();
Returns an integer
Has Interpolant been computed
This is a convenience method to determine whether or not the interpolant has been computed yet.
Get by
surfe.InterpolantComputed();
Returns a Boolean
Usage Example
#include <surfe_api.h>
#include <geo_builder.h>
#include <iostream>
using namespace std;
int main() {
Geo_Builder model;
model.InitializeGRBFInterpolantObject(4);
model.surfe->SetRestrictedRange(true, 10, 0);
model.SetInterfaceDataFile("well.csv");
model.surfe->SetRBFKernel("Thin Plate Spline");
try
{
model.LoadConstraintsFromFiles();
}
catch (const std::exception &e)
{
std::cout << "Surfe Exceptions: " << e.what() << " occurred. " << std::endl;
return EXIT_FAILURE;
}
try
{
model.surfe->ComputeInterpolant();
}
catch (const std::exception&e)
{
cout << "Exception: " << e.what() << endl;
return EXIT_FAILURE;
}
model.BuildRegularGrid(25,25);
model.WriteVTKIsoSurfaces("well_R3_RR_25_25.vtp");
model.WriteVTKEvaluationGrid("well_R3_RR_25_25.vti");
model.VisualizeVTKData();
return 0;
}