伴随优化算法格林函数

1. 李乐伟的关于格林函数

格林函数是积分方程方法的重要核心，包括矩量法和边界元方法。如果没有格林函数，上述两个方法就根本不能够被解决。因而是电磁理论与计算电磁学。候选人自1984年的硕士学位开始研究颂陆至今已经20年。 在场论野谨顷中，空间任意坐标下的偶极子，格林函数代表矢量场。如果 应用电路理论去等效，格林函数就好比多端口网络系统中的频响函数，如果知道了源，那么电场和磁场就可以应用格林函数来确定。由此，函数直接提供了源与它的辐射场的关系。
然而，在各种材料和不同的结构中建立格林函数公式就是相当复杂和深奥的数学,此外格林函数的快速计算形式把计算电磁学推向一个高度甚至更为重要。
3．解决了用于快速计算电大尺寸电磁散射问题的快速多极子算法、自适应积分方程方法和快速傅立叶变化方法，自主开发出三款商用软件包。
4．对负折射材料以及其他复合材料中的电磁场极化特性、散射特晌液性以及吸收特性进行了深 入研究，由于其在负折射材料领域的重大贡献和突出成就，作为唯一的外国人被美国空军遴选为相关重大研究项目的9名首席科学家之一。
5．在IEEE AP、IEEE MTT、Radio Science等微波界着名SCI国际期刊上发表论文近310多篇， 被SCI收录400多篇（包括专着章节和部分国际会议论文），被SCI引用1500多次。
6．已出版英文专着2本:
Spheroidal Wave Functions in Electromagnetic Theory, Wiley: New York, 320 pages, Wiley-Interscience Series, Edited by Kai Chang, ISBN No.: 0-471-03170-4, November, 2001. 《电磁学中的椭球波函数》。
Dyadic Greens Functions in Inhomogeneous Media, John Wiley & Sons, New York, 350, 2006 《不均匀介质中的并矢格林函数》。
完成了21本书中共45个章节的撰写，见电磁学研究进展系列丛书，PIER series (EMW Publishing, Boston, USA) 。
由于他的突出贡献，2004年当选为IEEE fellow。

2. 电磁场与电磁波领域，需要哪些数学知识作为基础

基本的预备知识是麦克斯组中的四个方程所涉及的物理学定律。库伦定律，磁场和电场下的高斯定律，法拉第电磁感应定律，以及毕奥萨伐尔扒祥毁春备环流定律（后来扩展成了安培环流定律。）物理学方面的基础知识就是这些。主要是数学方面的，需要掌握微积分相关知识，以及矢量相关的几个简单的运算法则。掌握一定的解微分方程，偏微分方程的方法（很对公式推导需要）。基本上地磁场和电磁波理论方面的学习就需要这些了。学电磁场与电磁波理论的书籍，肯定要有高等数学的基础，和数学物理做铺垫的！要不然里面的麦克斯韦方程积分形式，微分形式你就会很难看懂！所以你的选择很正确，把高等数学学了，再看一下数学物理的方法，还是可宴闹以自学电磁场与电磁波的，这么课还有一个特点就是公式太多，记一下慢慢消化。

3. 第四讲：弹性体模拟-弹性力学基础

这一讲的主题是弹性体模拟中的弹性力学基础，首先来看一个问题，就是当一个球跟多个球同时发生碰撞，如下图所示（称之为Bernoulli Problem）：

或者当敬历第一个球跟第二个球碰撞的瞬间，第二个球同时与第三个球碰撞，如下图所示（称之为Newton's Cradle）：

这些问题的求解都可以归结为线性互补问题Linear Complementary Problem（简称LCP）的求解，线性互补问题的定义可以参考 物理引擎之约束求解(一)——线性互补问题 。假设有两个刚体（分别标注为1号跟2号），由于每个刚体都有六个自由度（translation+rotation），那么这两个刚体的位置烂稿早就可以用下述公式来表示：

其中

同样 也有对应的变量表达，而由于刚体是不可穿插形变的，因此q还会有一些额外的约束：

这里的i指的是第i个碰撞点， 可以表示快要碰撞时候两个物体在第i个碰撞点处的signed distance。

如果用上一讲中的impulse-based方法来解方程，那么就意味着我们需要在碰撞的瞬间为两者施加一个冲量，而要想产生冲量，也就意味着两者产生了穿插，前面的条件将不再满足，即：

而如果我们希望在下一个timestep回复到正确的状态，那么上面这个函数的导数就需要是为正的：

这个表达式的物理意义是，等式后面表达式中的前面部分 表达的是碰撞点处的表面的法线（signed distance的梯度就是signed distance下降最快的方向，即表面法线方向），后面部分是刚体的速度，两者相乘大于等于0表示的是这两个向量的夹角小于九十度，也就是说在这个速度下，会使得两者的signed distance不断变大。

在impulse的作用下，速度会从 变成 ，即速度从碰撞变成背离。在Houdini中有个叫做restitution coefficient（恢复系数）的参数 ，这个参数是做什么的呢？用一个例子来说明，一个球在重力的作用下下坠，碰到地表，如果这个参数为0，那么这个小球在碰撞后，就会贴在地表上，如果是1的话，碰撞之后就会回弹到起始的高度，也就是说，这个参数表达的是经过碰撞之后能量的残留比例。

将恢复系数跟表面法线结合起来看，我们可以得到如下公式（这个条件我们称之为一号条件）：

如下图所示：

入射方向为 ，出射方向为 ， 可以简单看成是是出射方向的速度的长度与入射方向的速度长度的比例， 是法线方向，那么上面这个公式自然就是成立的，毕竟能量是守恒的，出射速度的长度不可能超过入射速度的长度。

上面这个公式是针对单个碰撞点（contact point）而言的，当我们有多个碰撞点的时候， 就变成了一个矩阵，这个矩阵我们称之为 （这个矩阵的物理含义为，将刚体的速度映射到碰撞点上的沿着法线方向的速度？不是特别理解，如果单个碰撞点的情况，对应的是碰撞点的法线，那么这里就对应的是多个碰撞点法线组成的矩阵，既然是多个碰撞点，也就没有一个整体的法线概念了），这里的转置后面会解释。

另外，这里还有另外一个公式：

这个公式中 表示的是刚体的质量， 是一个矩阵，这个矩阵就是上面梯度向量组成的矩阵的转置（推导这里就不说了，抛开多个碰撞的情况不说，单个碰撞点的情况下，冲量的作用方向肯定是垂直于碰撞表面的，也就是说经过冲量作用后的加速度应该也是与碰撞表面的法线方向保持一致的，从主观上来说，这个结果是可以理解的）， 则是前面说的碰撞时的冲量impulse，如果有多个碰撞点的话，这个变量也就是一个向量（几个碰撞点就是几维的）。矩阵 跟冲量的乘法得到的是一个力，而作用力除以质量就得到了加速度，通过这个加速度（在时间的作用下）就能完成速度的一个变化。

这个公式中，冲量（每个分量）必须要大于等于0（当刚体碰撞时的速度在表面法线方向上的速饥雀度非0，且速度法线方向的分量与法线方向相反，那么就取＞0，否则就取等于0），否则在碰撞后，刚体就直接装进另一个刚体中，从而发生穿插，这个条件我们称之为二号条件：

一号条件跟二号条件组成linear complementary condition（即如果 （对应的是碰撞点碰撞之后没有施加冲量，也就是说碰撞时候的速度与表面法线正好垂直，这种情况下，表面法线与速度正交，一号条件不等式前面的结果就是0）的时候，就不需要考虑一号条件；而如果 就必须要考虑一号条件，我们将这种情况称之为互补），最终我们需要求解的问题就可以表示为，在一二号条件：

约束下，求解：

即LCP，这个问题如果我们能够求解出冲量向量，那么我们就能得到需要的解，但是在多个碰撞点的情况下，这个会比较复杂。

Bullet物理引擎用的是一种叫做Gauss Siedal Method，一个个碰撞点去考虑，比如先考虑一号碰撞点的impulse是多少，之后再看一号impulse对二号碰撞点的速度的改变，之后再看二号碰撞点的impulse为多少并加上一号碰撞点的影响，同时算出对三号碰撞点的影响，以此类推，直到最后一个碰撞点对一号碰撞点的影响，不断迭代，经历多轮迭代最终达到一个平衡，这种方法是串行计算，对GPU不友好，性能较差。

另一种方法则是将LCP（这是一种timestep的方法）转换成一个优化问题，这种思想是物理模拟中的一种十分重要的思想，在很多问题的求解上都有应用。

先来说说拉格朗日乘子（ lagrange multiplier ），这是数学上的一种优化策略中的术语，这个优化策略用以在给定的一些等式约束下，求得某个函数的局部极值（极大极小值），这个算法的基本思想可以叙述为，对于一个需要求解极小值（极大值可以转换为极小值）的函数f(x)，已知需要满足g(x) = 0条件（或者说，已知极值是在g(x) = 0条件下满足），那么拉格朗日函数可以表示为：

我们知道，如果不考虑条件函数g(x) = 0的话，f(x)的极值可以直接通过 来求得，而在加上条件之后，问题就会变得复杂一点，上面的拉格朗日函数是自变量x跟 （这个就叫拉格朗日乘子）的函数，而原函数f(x)的极点则肯定会出现在拉格朗日函数的saddle point（极点）上，所以只需要对拉格朗日函数求偏导，并令各个偏导结果为0进行求解，就能得到取得极值的坐标点。

而如果对上面的情况进行泛化处理，比如将条件从等式变成不等式，那么我们就需要使用 KKT（Karush Kuhn Tucker） 条件方法，同样的，KKT会将原始的函数f(x)，与等式约束 以及不等式约束 写到一个式子里：

f(x)取得极值的条件为：

联合这些等式，就能得到最终的极值点的坐标。这里值得一提的是，上面第三个条件中，由于h(x)是小于等于0的，因此要想条件成立，要么h(x) = 0，要么 ，也就是说，这个条件可以拆成两个条件：

且这两个条件是complementary（互补）的，如果移除g(x)的干扰，这个条件就跟前面LCP问题的求解的形式完全一致了，也就是说，LCP的求解就变成了KKT（优化问题）的求解，而这种问题的求解方法就十分丰富了，这里就不做展开，否则可以直接讲一个学期。

但是我们发现，用LCP来求解Bernoulli Problem可以得到正确的结果，但是用来求解Newton Cradle问题的时候结果却是错误的，这是因为理论上LCP问题的求解是不唯一的，也就是说我们得到的解有多个，其中选取的那个可能并不是正确的解。

这里需要注意的是，上面的求解是没有考虑碰撞时的摩擦力的，实际上如果添加了摩擦力的话，问题会变得更为复杂。但是有意思的是，即使加上摩擦力，依然可以表达为LCP问题（不过更为复杂），转化成的优化问题，其限制条件将会变成二次函数，因此其求解也会变得更为复杂。

最后做个总结，在目前情况下的物理引擎，对刚体碰撞的模拟依然是存在较大的精度问题，比如说大多数都没有办法解决上面的Bernoulli跟Newton Cradle问题，即使是不考虑摩擦力的情况，目前都没有物理引擎能给出完善的解决方案，更何况是更为精确的考虑摩擦力的情况。

Continuum Mechanics中描述了形变的相关理论，包括弹性形变与塑性形变，我们这里着重考察弹性形变。

最简单的弹性形变就是前面说的弹簧+质点的粒子系统，当质点不断移动从而压缩弹簧的话，弹簧就会不断的积累势能，当手松开之后，势能可能就会转换为动能，此时弹簧产生的力叫做conservative force。

但是弹性形变不只是局限于弹簧系统，普通的物件也会发生此类形变，但是这里有亮点需要明确：

可以知道，弹性体的形变是一个连续的变化，即形变前相邻的点在形变后依然是相邻的。我们可以利用微分思想来对这个连续的规律进行考虑，假设形变前物体上某一点X（处在material space）经过形变后变成了点x（处在deformed space），那么形变就可以表示成如下的公式：

对这个映射关系采用泰勒展开：

由于 是一个三维的变量（三维空间中的点），因此上面公式中的偏微分项 就是一个3x3的矩阵，这个矩阵我们叫做deformation gradient（形变梯度矩阵），通常用F来表示。

从这个公式我们可知道，这个矩阵的第一列进行乘法运算时对应的是 的第一个元素，换句话说，第一列就是在material space中沿着x方向移动一个距离时，在deformed space中就需要沿着某个方向移动一个相同大小的距离，这个方向就用第一列向量来表示，同理，其他两列就分别对应y/z方向上移动时对应的deformed space中的移动方向。

接下来我们看看要如何用数学公式表示物体的弹性形变长度，假设在material space中某个点移动的位移为 ，那么移动的长度的平方就可以表示为 ，而在deformed space中对应点的移动位移就可以用 来表示，其长度则为 ，那么具体的形变量就可以用后者开平方后减去前者的开平方，不过这里为了避免因为开平方导致的高额计算消耗，可以先直接用平方之差来表示形变的测度：

上面等式中最后的 称之为Green（格林） Strain（形变） Tensor（本质上是个矩阵），这个矩阵可以用来描述微分意义下物体形变的规律，可对于推导或者计算出这个 的点X而言，在其周边较小的一块区域（微分算子）里的所有点都可以通过这个矩阵来计算其形变量，但是，如果超出这个区域范围， 可能会不同，也就是说， 准确来说是物件上某点的一个函数，可以写成

接下来，就是如何将前面的 变成势能，势能可以表示为如下形式的函数：

我们知道，这个函数有如下的性质：

根据上面两个性质，我们可以推断，这个函数可以大致用下图来表示：

至少是一个两阶以上的模型，应用泰勒展开，可以得到：

而由于 跟 都是0，那么上面公式就变成了：

由于 是个矩阵，那么这个公式就可以改写成如下的形式，这个公式我们称之为constitution model：

到目前为止，所有的公式都是通过数学推导给出的，没有掺杂物理相关的概念与理论，接下来，为了能够得到上面的公式的具体形式，就需要给出系数 ，而物理学中针对这个系数提出了众多的模型，甚至还有人考虑通过神经网络训练得到这个系数，这就造成了不同的势能物理模型。

要知道这个公式是如何在物理中进行应用的，我们就放到下一讲中进行介绍了。

4. 反演技术是什么

（1）首先介绍下什么是反演技术。
①水平多层土壤的视在电阻率正演：是指已知水平多层土壤参数和四极法极距布置参数，计算不同四极法电极布置情况下的视在电阻率。水平多层土壤的视在电阻率正演模型是水平多层土壤的视在电阻率反演的基础，反演其实就是由优化方法控制的正演迭代。
②接地理论中定义格林函数为单位点源对应的空间电位函数，而使用格林函数求解接地问题的方法称为格林函数法。
（2）土壤反演技术
①土壤参数反演问题属于小规模的优化问题(变量一般只有几个)，但是土壤参数反演问题的非线性程度极高，基本上所有的优化方法都要花费大量的计算时间进行求解空间的搜索，同时也容易陷入局部解。
②迭代初值对于土壤反演来说是一个十分重要的参数。虽然没有免费午餐定律说明所有算法对所有优化问题的总体求解效率是相等的，但是对于土壤反演这种特殊问题，使用不同方法进行对比和分析后仍可以发掘和提出更多更有效方法。
（3）传统优化方法
①鉴于现场土壤电阻率结构千差万别，现场测量数据的反演解释是整个反演过程最困难和抽象的部分。优化方法结合水平多层土壤视在电阻率的计算模型是反演过程的核心，
②优化函数库的优化方法可分为以下两种方法：
Ⅰ、直接搜索法：只需要使用目标函数值的方法。
Ⅱ、梯度法：可能需要目标函数的一阶或者二阶偏导数信息的表达方法。

5. HFSS算法及应用场景介绍

安氏

前言

相信每一位使用过HFSS的工程师都有一个疑问或者曾经有一个疑问：我怎么才能使用HFSS计算的又快又准？对使用者而言，每个工程师遇到的工程问题不一样，工程经验不能够直接复制；对软件而言，随着HFSS版本的更新，HFSS算法越来越多，针对不同的应用场景对应不同的算法。因此，只有实际工程问题切合合适的算法，才能做到速度和精度的平衡。工程师在了解软件算法的基础上，便能够针对自己的需求进行很好的算法选择。

由于当今世界计算机的飞速发展，让计算电磁学这门学科也有了很大的发展，如图1所示，从大的方面来看，我们将计算电磁学分为精确的全波算法和高频近似算法，在每一类下面又分了很多种算法，结合到HFSS软件，通过ANSYS公司40余年来坚持不懈的研发和战略性的收购，到目前为止，HFSS有FEM、IE（MoM）、DGTD、PO、SBR+等算法，本文会针对每种算法和应用场景逐一介绍，相信你看完这篇文章应该对HFSS算法和应用场景会有更深的认识。

算法介绍

全波算法-有限元算法（ FEM）

有限元算法是ANSYS HFSS的核心算法，已有二十多年的商用历史，也是目前业界最成熟稳定的三维电磁场求解器，有限元算法的优点是具有极好的结构适应性和材料适应性，充分考虑材料特性：趋肤效应、介质损耗、频变材料；是精确求解复杂材料复杂结构问题的最佳利器，有限元算法采用四面体网格，对仿真物体能够很好的进行还原。

FEM算法的支配方程见下图：

HFSS有限元算法在网格划分方面能够支持自适应网格剖分、网格加密、曲线型网格，在求解时支持切向矢量基函数、混合阶基函数和直接法、迭代法、区域分解法的强大的矩阵求解技术。

在应用领域，HFSS主要针对复杂结构进行求解，尤其是对于一些内部问题的求解，比高速信号完整性分析，阵列天线设计，腔体问题及电磁兼容等应用场景，非常适合有限元算法求解。

有限元算法结合ANSYS公司的HPC模块，ANSYS HFSS有限元算法可以进行电大尺寸物体的计算，大幅度提升仿真工程师的工作效率。针对宽带问题，FEM推出了宽带自适应网格剖分，大大提升了仿真精度。

全波算法-积分方程算法（ IE）

积分方程算法基于麦克斯维方程的积分形式，同时也基于格林函数，所以可自动满足辐射边界条件，对于简单模型及材料的辐射问题，具有很大的优势，但原始的积分方程法计算量太大，很难用于实际的数值计算中，针对此问题， HFSS 中的 IE算法提供了两种加速算法，一种是 ACA 加速，一种是 MLFMM，分布针对不同的应用类型。 ACA 方法基于数值层面的加速技术，具有更好的普适性，但效率相比 MLFMM 稍差， MLFMM 算法基于网格层面的加速，对金属材料，松散结构，具有更高的效率。

IE算法的支配方程见下图：

IE算法是三维矩量法积分方程技术，支持三角形网格剖分。IE算法不需要像FEM算法一样定义辐射边界条件，在HFSS中主要用于高效求解电大尺寸、开放结构问题。与HFSS FEM算法一样，支持自适应网格技术，也可以高精度、高效率解决客户问题，同时支持将FEM的场源链接到IE中进行求解。HFSS-IE算法对金属结构具有很高的适应性，其主要应用领域天线设计、天线布局、 RCS、 EMI/EMC仿真等方向。

高频近似算法-PO算法

FEM算法和IE算法是精确的全波算法，在超大电尺寸问题上，使用精确全波算法会造成效率的降低。针对超大电尺寸问题，ANSYS推出PO（物理光学法）算法，PO 算法属于高频算法，非常适合求解此类问题，在适合其求解的问题中，具有非常好的效率优势。

PO算法主要原理为射线照射区域产生感应电流，而且在阴影区域设置为零电流，不考虑射线追迹或多次反射，以入射波作为激励源，将平面波或链接FEM（IE）的场数据作为馈源。但由于不考虑射线的多次反射和绕射等现象，一般针对物理尺寸超大，结构均匀的物体电磁场计算，在满足精度的要求，相比全波算法效率明显提高。比如大平台上的天线布局，大型反射面天线等等。

高频近似算法-SBR+算法

PO算法可以解决超大电尺寸问题的计算，但由于未考虑到多次反射等物理物体，主要用于结构均匀物理的电磁场计算。针对复杂结构且超大电尺寸问题，ANSYS通过收购Delcross公司（Savant软件）引入了SBR+算法， SBR+是在SBR算法（天线发射出射线，在表面“绘制” PO电流）的基础上考虑了爬行波射线（沿着表面追迹射线）、物理绕射理论PTD（修正边缘处的PO电流）、一致性绕射理论UTD（沿着边缘发射衍射射线，绘制阴影区域的电流），因此SBR+算法是高频射线方法，具有非常高效的速度，同时具有非常好的精度，在大型平台的天线布局中效果非常好。

SBR+支持从FEM、IE中导入远场辐射方向图或者电流源，也支持导入相应的测试数据，SBR+算法主要用于天线安装分析，支持多核、GPU等并行求解方式并且大多数任务可在低于8 GB内存下完成。

混合算法（ FEBI， IE-Region，PO-Region，SBR+ Region）

前面对频率内的各种算法做了介绍并说明了各种算法应用的场景，很多时候碰到的工程问题既包括复杂结构物理也包括超大尺寸物理，如新能源汽车上的天线布局问题，对仿真而言，最好的精度是用全波算法求解，最快的速度是采用近似算求解，针对该问题，ANSYS公司将FEM算法、 IE 算法、PO 算法、SBR+算法等融合起来，推出混合算法。在一个应用案例中，采用不同算法的优点而回避不同算法的缺点，可极大限度的提高算法的效率，以及成为频域内解决大型复杂问题的必备算法。

HFSS中FEM与IE可以通过IE Region与FEBI边界进行混合求解，FEM与PO、SBR+算法可以通过添加PO Region及SBR+ Region进行混合，混合算法的使用扩大了HFSS的使用范围。

时域算法-transient算法

HFSS时域求解是基于间断伽略金法（discontinuous Galerkin method, DGTD）的三维全波电磁场仿真求解器，采用基于四面体有限元技术，能得到和HFSS频域求解器一样的自适应网格剖分精度，该技术使得HFSS的求精精度成为电磁场行业标准。这项技术完善了HFSS的频域求解器技术，帮助工程师对更加深入详细了解其所设计器件的电磁性能。

Transient算法支配方程见下图：

采用HFSS-Transient算法，工程师可利用短脉冲激励对静电放电、电磁干扰、雷击和等应用问题开展研究，还包括时域反射阻抗以及短时激励下的瞬态场显示也可以借助它来完成。

谐振分析-Eigenmode算法

谐振特性是每个结构都存在固有的电磁谐振，谐振的模式、频率和品质因子，与其结构尺寸相关，这些谐振既可能是干扰源的放大器，也可能是敏感电路的噪声接收器。谐振会导致信号完整性、电源完整性和电磁兼容问题，因而了解谐振对加强设计可靠性很有帮助。

Eigenmode算法支配方程见下图：

在HFSS中，使用eigenmode算法可计算三维结构谐振模式，并可呈现图形化空间的谐振电压波动，分析结构的固有谐振特性。依据谐振分析的结果，指导机箱内设备布局和PCB层叠布局，改善电磁兼容特性。

总结

HFSS里面有各种不同的算法，有全波算法、近似算法以及时域算法，工程师可以格局需要选择不同算法（最高的精度和最高的效率）。首先针对频域算法，使用范围见图14，通常FEM算法和IE算法非常适合于中小尺寸问题，对大型问题，FEM/IE运行时间/内存需求非常巨大； PO方法适合解决超大电尺寸问题，但对问题复杂度有限制，通常通常不能提供客户所期望的精度，但对于均匀物体是一个很好的选择；SBR+算法适合解决超大电尺寸问题，对复杂结构也能够提供很好的精度和速度；针对既有电小尺寸复杂结构计算问题，又有电大尺寸布局计算问题，混合算法是一个很好的选择。Transient算法适合解决与时间相关的电磁场问题，如ESD、TDR等；Eigenmode算法专门针对谐振仿真。

想要更多，点击此处，关注技术邻官网

6. 二维网格数值插值技术

大部分油气藏的数据是散乱分布的，因此称为散乱数据。散乱数据指的是在二维平面上或三维空间中，无规则的、随机分布的数据。利用散乱数据建模就要对散乱数据进行插值或拟合。

设在二维平面上有n个点 （x_i，y_i） （i ＝1，…，n），并有Z_i ＝f（x_i，y_i）。插值问题就是要构造一个连续的函数F （x，y），使其在 （x_k，y_k） （后＝1，…，n） 点的函数值为Z_k，即Z_k ＝f（x_k，y_k） （k＝1 ，…，n）。

早在20世纪60年代，散乱数据的插值问题就已引起人们的注意。近50年来，已经有多种算法被提了出来。但是，由于应用问题千差万别，数据量大小不同，对连续性的要求也不同等等，没有一种算法适用于所有的场合。而且大多数算法只能适用于具有中、小规模数据量的散乱点插值问题。大规模散乱数据 （例如，10000个点以上） 的插值问题还正在研究之中。

据散乱数据的复杂程度，其可分为单自变量、双自变量及多自变量3种类型。下面将主要讨论双自变量散乱数据的插值问题。

（一） 插值的一般概念

插值的概念最早可追溯到 “控制论之父” 诺伯特·维纳的不朽着作 《平稳时间序列的外推、插值和光滑及其工程应用》 （Wiener，1949）。随着计算机技术的发展，插值的概念已广泛地应用于数据处理、自动控制、数值分析、地球物理及数学地质等领域。

1. 插值方法

计算机插值方法一般可分为两大类：拟合函数方法和加权平均方法。它们的原理都是来自手工方法。Crain （1970） 把这两种方法得到的曲面分别称为数学曲面和数值曲面。Alfeld ＆ Barnhill （1984） 分别称它们为分片方法和点方法，而Cuyt （1987） 则把这两种方法分别称为系数问题和数值问题。

利用拟合函数方法进行插值，就是利用二元多项式来表示一个插值曲面，插值问题化为确定这个二元多项式的系数的问题。一般来说，往往可通过求解一个线性代数方程组来获取这些系数。这个方程组的系数可由观测数据来确定，它们代表了观测数据的影响，而其方程的次数则表示了控制多项式拟合程度。方程次数越高，拟合的程度就越高。当这个多项式的系数确定以后，将空间某一点处的坐标代入该多项式，即可求得该点处的值。

这个方法的特点是可以制服畸变的原始数据或带有噪声的原始数据。所以，用一个函数进行拟合是一个光滑的过程，一些局部的细节可能消失。所得的插值曲面的复杂程度取决于多项式的次数，即所求解的线性方程的数目。

加权平均插值方法把求插值的问题化为求取观测数据的加权平均。每个观测数据点对应的加权系数恰恰反映了该数据点对插值点的影响大小。为了求取一个插值，必须要计算出一组加权系数。

加权平均方法的一个主要优点是，可以获取在观测数据点附近变量的小尺度趋势。而利用一个适当次数的多项式是无法获取曲面的这种局部细节的。

从原则上讲，这两种方法的差别就在于：加权平均方法强调了曲面的局部细节，而拟合函数方法则概括了曲面的整体性质。从计算时间上看，前者花费的时间比后者要多得多。

2. 插值效果评判

从理论上讲，一个插值方法的效果如何应通过插值结果和客观存在的原始曲面的比较，按以下3条标准来进行判断：

（1） 原始曲面和插值结果之间差异的最大值为最小。

（2） 原始曲面和插值结果之差的平方和为最小。

（3） 在每个观测数据点处，插值结果本身的数值及其1阶到k阶导数和原始曲面的相等。

由于原始曲面本身是未知的，所以在以上3个标准中，第一个和第二个标准是无法检验的，仅有第三个标准是在一定的模型假设之下可以进行检验。在一些实际应用中，当原始曲面可用解析函数来表达时，仅利用第三个标准来检验插值的效果也是可行的。然而，在地质建模中，观测数据点往往不够多，且还有一定的观测误差，不可能断定原始曲面是否可用解析函数表达。这时，插值技术的合理性必须从直观的几何和人们的经验等方面进行评价。

利用计算机进行插值所遇到的困难，主要来自观测数据点数目不足和观测误差。如果观测数据充分多且精确，那么几乎所有的插值方法都会给出良好的效果。另一方面，对于圆形或狭长的隆起，凹陷和鞍点等变量的空间变化几何特征，在数据点分布较稀的情况下，用任何插值方法都是难以推断出它们的存在的。所以，插值方法需要考虑曲面的局部斜率的影响。

下面主要介绍加权平均方法和拟合函数方法中最常用的插值算法。这些算法能解决大部分油气藏建模问题。

（二） 与距离成反比的加权法

距离成反比加权插值方法是基于如下的模型：每个数据点都有局部影响，这个影响随着数据点和插值点距离的增加而减弱，且在一定的范围以外，可以忽略不计；这个影响是以该数据点为中心；而在任一点处的插值恰是各数据点影响之和。

1. 与距离成反比加权插值公式

这一方法首先是由气象学及地质学工作者提出来的，后来由于D. Shepard的工作被称为Shepard方法。其基本思想是将插值函数F （x，y）定义为各数据点函数值f i的加权平均，即：

油气田开发地质学

在 （x_k，y_k） 点处函数值可写成：

油气田开发地质学

式中： 表示由第i个（x_i，y_i）点到插值点（x_k，y_k）的距离；W_i（x_k，y_k）——权函数；μ——功率因素，通过改变值来调整权函数与距离的关系，与距离成反比加权和与距离成平方反比加权分别是μ＝1和μ＝2时的特殊情形。

与距离成反比加权插值方法是最早使用的计算机插值方法，至今仍被广泛地应用着。在大多数商业性的等值线图绘制软件包中被用来形成网格化数据。这种方法较为直观：一个数据点对于插值点的影响模型化为与这两点之间的距离成反比。

2. 与距离成反比加权插值改进

距离成反比加权算法中的功率因素μ应该取为μ≥0，否则表明距离越远的点作用越突出，这违背了普通常识。考虑以下最极端的情况是：

油气田开发地质学

假设有n个数据点，一个插值点 （x_k，y_k） 位于第i个数据点附近，相应的权函数可写成：

油气田开发地质学

式中：d_j （x_k，y_k），d_i （x_k，y_k）——插值点 （x_k，y_k）到各数据点的距离。

当插值点 （x_k，y_k） 和第i个数据点很靠近时，可以认为其他数据点对W_D的影响是一个常数，即C为常数。当该插值点和第i个数据点的距离趋于零时，对于不同的μ，W_D会有不同的性质。首先看W_D对D的导数，有：

油气田开发地质学

然后再有：

油气田开发地质学

可见，当μ＝1时，随着D趋于0，W′_D近似为不随D变化的常数，这可用图6-8中左端的图形来表示。

当μ>1时，随着D趋于零，W′_D趋于零。这说明在该数据点附近，加权系数的变化为零，即可用图6-8中间的图形来表示。

当μ<1时，随着D趋于零，W′_D趋于无穷大。这说明加权系数在该数据点附近还有一个尖点，如图6-8右端的图形所示。

图6-8 与距离反比加权的权数随参数μ的变化

（1）功率因素μ越小时，近距离点和远距离点的作用越接近，生成的平面网格数据越平滑。随着μ增大，平面网格数据的光滑性越差。同时，μ直接影响网格数据的极值和均值，μ越小网格数据的均值越接近原始数据的均值，但极值相差越大。因此，当要求插值结果尽可能接近原始数据的均值时，功率因素不能选择过大。例如对于开发早期的油气藏，因为仅有少数探井控制，此时网格化得到的各类物性参数应该在总体上符合井点的统计结果，均值是比极值更有价值的参数，因而通常将功率因素取为1。相反，μ取值越大，网格数据越能恢复原始数据的极值，但也容易使均值误差增大。原因是当μ取较大的值时，近距离点的作用越突出，原始数据点分布的不均匀性使得部分点在网格节点上发挥了更大的作用，而另外一些点的作用则受到屏蔽。因此，当要求突出数据的局部特征，体现储层的非均质性特征时，功率因素应选择得大些，一般取为2。

（2） 利用与距离成反比加权法进行插值时，当增加、删除或改变一个点时，权函数W_i （x_k，y_k） 均需重新计算，因而该方法是一个全局插值算法。

为了克服Shepard方法的上述缺陷，Franke及Nielson提出了MQS （Modified QuadraticShepard） 方法，它仍然是一个与距离成反比的加权方法。对它的改进如下：

插值点 （x_k，y_k） 到已知数据点的距离d_i作适当修改，使其只能在局部范围内起作用，以改变Shepard方法的全局插值性质。这时重新定义距离函数：

油气田开发地质学

式中：r_w为一个常数。而

油气田开发地质学

因此，当 （x_k，y_k） 点与某一点的距离大于r_w时，权值就为零。

（3） 与距离成反比加权法的权函数W_i（x_k，y_k）始终满足W_i（x_k，y_k）≤1/n，因此插值结果不会大于或小于原始数据的最大值、最小值。当已知数据点过少时 （这种情况在早期地质研究中是最为常见），使用具有外推能力的曲面样条或趋势面分析，得到的结果往往背离实际。其原因是这两种方法在远点不具有控制能力，它将沿趋势无限发展下去。特别是在仅有少数井资料可用的情况下，与距离成反比加权法应是优先选择的方法。

但是，隐含在原始数据中的尖峰会被淹没而无法显示出来。因为距离成反比加权进行插值时，每个数据点所发挥的作用基本上是中心对称的，因此对山脊和山谷等非各向同性的几何形状的显示不利。为了克服这种状况，需要考虑变量的局部变化趋势，为此需要对梯度进行估计。

当已知数据点用与距离成反比加权方法形成数据插值曲面时，该曲面被称数据曲面。与距离成反比加权方法也可用于各个数据点处的切平面，把任一点处的插值值取成各切平面在该点处取值的一个加权平均，所形成的曲面称为与距离成反比加权梯度插值曲面，简称梯度曲面。如果在数据点以外的一个点是变量的局部高点，那么梯度曲面在该点的值容易大于变量的真实值，即呈现 “过估计” 的状态。如果数据曲面在该点的值小于变量的真实值，则呈现 “欠估计” 的状态。因此，数据曲面可以通过和梯度曲面的相互结合来克服本身的缺陷。

可以用数据曲面和梯度曲面之差乘以一个系数作为一个修正量，对数据曲面进行修正，这样，可以用下述的曲面来代替单纯的数据曲面：

油气田开发地质学

式中：L （x，y） ——距离反比加权插值；W_i——和第i个数据点的距离反比加权系数；τ_i——曲面在该点处的粗糙度指数；S_i（x，y）——第i个数据点 （x_i，y_i） 处的切平面在点（x，y） 处的值；H（W_i，τ_i）——混合函数。

如此得到的曲面称为混合曲面。它通过所有的数据点，具有连续的坡度，其变化在空间的分布更均匀。数据曲面和梯度曲面是混合曲面的两种极端情况。由于数据曲面和梯度曲面之差在各数据点处为零，还因为混合函数的变化范围为0～1，且当混合函数等于0或1时，其一阶导数为零，故混合曲面和梯度曲面相切于各数据点处，且其高阶导数在各数据点处亦为零。

（4） 与距离成反比加权法仅考虑了插值点与数据点之间距离的影响，没有考虑到各数据点之间的关系，物性参数分布的趋势性没有得到充分的体现。为此，Franke及Nielson进行了改进。

用节点函数Q_i （x，y） 代替f_i，Q_i （x，y） 是一个插值于 （x_i，y_i） 点的二次多项式，即有Q_i （x_i，y_i） ＝f，i＝1，…，n。Q_i可由下式表示：

Q_i（x，y）＝f_i＋a₁（x-x_i）＋a₂（y-y_i）＋a₃（x-x_i）2＋a₄（x-x_i）（y-y_i）＋a₅（y-y_i）²

式中：a₁，a₂，…，a₅是按下式最小二乘法得出的优化解：

油气田开发地质学

式中：f_i，f_j分别为 （x_i，y_i）和 （x_j，y_j）点的函数值，而ρ_j可按下式选取：

油气田开发地质学

其中r_q为一常数，而

油气田开发地质学

求出Q_i （x，y） 后，插值函数可表示为：

油气田开发地质学

上述方法消除了Shepard方法中的一些缺陷，因而在散乱点插值中得到广泛的应用。但是，为了求得Q_i （x，y） （i＝1，…，n），需要多次求解线性方程组，计算量大，因此，一般只用于中、小规模散乱点的插值运算。

（三） 多项式趋势面法

由计算机产生的曲面一般不会总是和原始的观测数据一致。如果两者的差别在给定的尺度之下不是很明显，那么产生的曲面可被认为是插值曲面，否则就被认为是近似曲面。如果观测数据含有明显的观测误差，近似曲面就显得更合理。这时，和插值曲面相比，近似曲面由于数据的各种误差所产生的扰动不太容易看得清，但是近似曲面空间变化的一些主要性质还是能清晰地被体现出来的。

近似曲面和每个数据点之间的差称为残差，可视为每一个数据点上的一种误差表示。然而，计算出来的这种残差是意味着对未知的观测误差的一种度量，还是意味着一种允许的插值误差，或者意味两者都是，这要依变量的空间性质和观测数据的获取方法而定。

确定近似曲面的方法可分为3种。第一种方法是以残余的平方和最小为条件，确定多项式的系数，以获取曲面。第二种方法是利用观测数据误差的附加信息，并满足最小曲率的原则以确定曲面。最后一种方法是利用观测误差和插值误差的附加信息，以满足最小平方差或最小曲率为条件确定曲面。以下主要讨论多项式构造趋势面法。

多项式构造趋势面是目前最常用的方法，一次多项式表示的趋势面是空间的一个平面，二次趋势面是抛物面，椭球面或双曲面，三次及三次以上的趋势面是形态复杂的空间曲面，随着趋势面的次数增高，曲面的形态就越复杂。

如果有一组总共n个观测数据，其观测点的平面坐标为 （x_i，y_i），地质变量的观测值为f_i（x_i，y_i）。对于这组观测数据的多项式趋势面方程表示成如下形式：

油气田开发地质学

式中： ——第i个观测点的趋势值；a₁，a₂，…，a₅——待定系数，它们的个数m与所选用趋势面方程的多项式次数n存在下列关系式：

m＝［n（n＋3）＋2］/2

为使趋势面最大限度地逼近原始观测数据，可采用最小二乘法使每个观测点的观测值与趋势值之差 （残差） 的平方和最小，即：

油气田开发地质学

得到需求解的m阶正规方程组。当系数矩阵满秩时，趋势面方程也就被唯一确定了。

由于趋势面分析不具有过点性，使得局部井点上误差可能很大。同时参数场在三维空间中的分布过于复杂，无论从理论上还是实验中都无法确证某类参数场能较好地符合某确定次数的曲面，多项式次数过高，会导致趋势面发生频繁振动，多项式次数过低，得到的趋势面又过于光滑，丧失许多细节。再者，趋势面方程在外推过程中容易使参数场发生畸变，产生无意义的结果。因而现代地质建模研究中已经很少将趋势面分析单独作为插值方法使用。但对于下列两种情况趋势面分析仍然能达到较为理想的效果，一是对小范围内具有显着趋势性分布的数据点；二是数据点分布过于密集，而且可能存在若干异常数据点时，趋势面分析会自动削弱异常数据点的影响。为提高趋势面分析的精度，可采用残差来校正趋势面分析的结果。具体过程如下：

（1） 由趋势面分析得到任一网格节点 （x_k，y_k）处的趋势值

油气田开发地质学

（2）计算n个已知点 （x_i，y_i）处的残差△f_i（x_i，y_i）＝f_i（x_i，y_i）

油气田开发地质学

（3） 调用某种插值方法将残差分配到每个网格节点上，对网格节点 （x_k，y_k） 有△f_k（x_k，y_k）；

（4） 网格节点 （x_k，y_k）经校正后的最终结果为

油气田开发地质学

特别地，如果残差分配的插值算法也是趋势面分析，就形成所谓的多级趋势面分析。此时，网格节点上的值是多次趋势面分析的结果，多级趋势面分析在不断减少观测点插值误差的同时，整张参数场曲面仍然保持其连续性和光滑性，原因是该算法同样符合线性迭加原理。

（四） 径向基函数插值法

径向基函数的名字来源于这样一种情况，即基函数是由单个变量的函数构成的。一个点 （x，y） 的这种基函数的形式往往是h_k（x，y）＝h（d_k），这里的d_k表示由点 （x，y） 至第k个数据点的距离。一般说来，这种方法不具有多项式精度，但只要稍加改进，即可获得具有多项式精度的插值公式：

油气田开发地质学

式中：q_k（x，y）是一个多项式基，其阶次小于m。

上式中的系数a_k和b_k应满足下面的联立方程组：

油气田开发地质学

油气田开发地质学

第一式中的n个方程式满足了插值要求，而第二式中的m个方程式则保证了多项式精度。两式中共有m＋n个未知数，同时存在m＋n个方程式，联立求解，即可得出待定系数。

下面，介绍两种主要的径向基函数插值法。

1. Multiquadric方法

Multiquadric方法是由R. L. Hardy在1971年提出来的。它是最早被提出并且应用得最为成功的一种径向基函数插值法。它采用的插值函数，即 （x，y）处的值F （x，y）：

油气田开发地质学

式中： 为基函数；e_i——非负常数；a_i——加权系数，满足如下方程组：

MV_a＝V_z

式中：V_a=(a₁,a₂,…,a_n)^T,V_z=[f(x₁,y₁),f(x₂,y₂),…,f(x_n,y_n)]^T,f(x_i,y_i)是(x_i,y_i)处的数据点的值。

油气田开发地质学

由于M和V_z都不依赖于插值点的坐标 （x，y），所以a_i也不依赖于 （x，y）。然而，基函数C（X-X_i）则是以Xi为参数的 （x，y）的函数。所以说，插值曲面F（x，y）是n个基函数C（X-X_i）所构成的n个空间曲面配置而成的。

Arther提出了如下形式的基函数：C（d）＝1-d²/e²，其中d是数据点到插值点之间距离，而e则是一个常数。显然，基函数C（d）是d的一个衰减函数。当d＝0时，C（d）取得最大值1，而当d≤e时，有0≤C（d）≤1。这时，基函数呈现为椭圆抛物面，而加权系数a_i（i＝1，2，…，n）所满足的线性代数方程组的矩阵M应作相应的改动，其对角线元素应改成1。Hardy（1971）引入如下的基函数：C（d）＝（d²＋e²）^1/2，该基函数呈现为椭圆双曲面。Hardy还建议将e²取成0.815乘以数据点间距离的平均值。

2. 薄板样条法

样条 （Spline）本来是绘图员用来绘制光滑曲线的工具，是一种用木材或金属等弹性材料做成的细条。在绘图时，沿着通过图纸各已知点的样条，便可绘出一条光滑曲线。数学上所说的样条 （多项式样条） 实质上是分段多项式曲线的光滑连接。当函数为分段的m次多项式，在分段点上有直至m-1阶连续导数，那么该函数则称为m次样条函数，简称为样条。一般来说，研究和应用得比较多的是三次样条。零次和一次样条函数分别是台阶状函数和折线状函数。以上所述的是关于一维样条函数，对于二维样条函数也可作为类似的考虑。

样条函数的主要功能是进行插值，其主要优点在于，能在插值多项式的次数尽可能低的条件下，使插值曲线或插值曲面取得较高的光滑度，且只需要利用函数本身的值，而不需要提供函数的各阶导数的值。

三次样条曲面包含有三种不同的类型：双三次样条、伪三次样条及薄片样条。这些样条曲面以m和s为其两个参数，使得希氏空间Hs中元素的m阶导数的范数所构成的一个泛函达到最小。此外，这个泛函具有旋转不变性。

对于薄片样条曲面，m＝2，s＝0。这一方法是由R.L. Harder及R. N. Desmarais在1972年提出来的，后来由J. Duchon及J. Meinguet等人予以发展。薄板样条法得名于如下事实，即用此方法求出的散乱点的插值函数使下面这一泛函表达式具有最小值：

油气田开发地质学

在这里，I（F）表示受限于插值点的无限弹性薄板的弯曲能量。因此，这一方法的实质从力学观点看是使插值函数所代表的弹性薄板受限于插值点，并且具有最小的弯曲能量。这是一个泛函求极值的问题。这一变分问题的解即为我们所需要的插值函数，具有径向基函数插值法的一般形式。

R.L. Harder及R. N. Desmarais提出解析形式如下：

油气田开发地质学

且有 其中t（x，y）和t_i（x_i，y_i）是二维空间中的点，而f_i是t_i处的观测值。还有，K（t_i，t） 这里，r_i代表点t和t_i之间的距离，

由解析表达式及其约束条件，可给出用以确定系数的线性代数方程组：

油气田开发地质学

其中，K=[k_ij]_n×n,k_ij=K(t_i,t_j),k_ii=0,F^T[f₁,f₂,…,f_n],α^T=[b,a₁,a₂],A^T［λ₁，λ₂，…，λ_n］，且有：

油气田开发地质学

求解上述n＋3阶方程组则得到待定系数b，a₁，a₂，λ_i，然后可插值出平面任一位置的函数值F（x，y）。

上述方程与Enriguez等人给出的方程相类似，其微小的差异就是基函数中的自然对数（In） 变成了常用对数 （lg）。

在构造薄片样条曲面的过程中，Franke （1982）提出了以r²lgr作为基函数，Sandwell（1987）则提出了以双调和格林函数r²（lgr-1）作为基函数。另外，Ayeni （1979）也对不同的基函数进行了讨论。

使用平面上n个已知点进行曲面样条插值时，实质上是求解一个n＋3阶线性方程组以确定n＋3个系数。为了保证解的存在性和唯一性，系数矩阵应该是满秩的。对下列3种情况必须避免：

（1） 在给定的n个已知点中存在着距离过近的两点 （x_i，y_i） 与 （x_j，y_j） 极端的情形是同一个数据点的重复输入，此时系数矩阵的第i行与第j行对应各元素非常接近，导致线性方程组的系数矩阵是奇异的。因此在数据预处理过程中必须消除沉余数据。

（2） 给定的已知数据点过少，此时系数矩阵的后3行线性相关，矩阵是不满秩的。

（3） n个已知点数据分布在一条直线上，显然由这样的n个点不能唯一决定一张曲面，系数矩阵表现为不满秩。

针对情况 （1），通常在做曲面样条插值时，首先对数据进行预处理，通过给定一个适当的距离下限r_min来滤掉那些相距过近的点，研究发现r_min＝（△x＋△y）/8是一个较合理选择 （△x和△y分别表示x和y方向的步长）。情况 （2） 和 （3） 实际上意味着不能进行曲面样条插值，除非通过数据均整来改变数据分布状态。

曲面样条插值方法是一种严格的过点插值法，即由生成的样条曲面必定通过给定的n个已知数据点，这样井点数据的控制作用自然得到体现。同时，曲面样条方法充分考虑了数据间的相对位置，其插值精度很高，在外推过程中，总是沿数据点的分布趋势外推，因此曲面样条法是具有一定外推能力的插值方法。同时曲面样条方程得到的是一张连续光滑的曲面。

曲面样条插值方法特别适合于地层层面的生成和地层厚度的插值。如果从曲面样条法严格的过点性、良好的光滑性及外推性看，适用于那些光滑、趋势性明显、变化连续的储层物性参数诸如油气饱和度的插值。

曲面样条法插值的精度很大程度上取决于数据点分布的均匀程度，稀疏区域主要由邻近区域的外推得到，其插值结果可能偏差较大。同时，曲面样条法的严格过点性使得它不能分别对待不同精度点的数据，因此它不具备数据的校正能力。