伴隨優化演算法格林函數

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個已知數據點，這樣井點數據的控製作用自然得到體現。同時，曲面樣條方法充分考慮了數據間的相對位置，其插值精度很高，在外推過程中，總是沿數據點的分布趨勢外推，因此曲面樣條法是具有一定外推能力的插值方法。同時曲面樣條方程得到的是一張連續光滑的曲面。

曲面樣條插值方法特別適合於地層層面的生成和地層厚度的插值。如果從曲面樣條法嚴格的過點性、良好的光滑性及外推性看，適用於那些光滑、趨勢性明顯、變化連續的儲層物性參數諸如油氣飽和度的插值。

曲面樣條法插值的精度很大程度上取決於數據點分布的均勻程度，稀疏區域主要由鄰近區域的外推得到，其插值結果可能偏差較大。同時，曲面樣條法的嚴格過點性使得它不能分別對待不同精度點的數據，因此它不具備數據的校正能力。