賈艷艷，邢學(xué)軍，史基安，孫國強(qiáng)

1中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，蘭州 730000

2中國科學(xué)院大學(xué)，北京 100049

3中國石油冀東油田公司，唐山 063004

1 引言

目前我國油氣勘探已轉(zhuǎn)移到西部地區(qū)，西部地區(qū)劇烈起伏的地形給地震勘探提出了嚴(yán)峻的挑戰(zhàn)，經(jīng)典反射地震勘探技術(shù)在復(fù)雜地表地區(qū)不再適用，地震波數(shù)值模擬中的網(wǎng)格法由于具有對地形和介質(zhì)靈活劃分的特點(diǎn)成為解決此類問題最有效的方法.網(wǎng)格法以結(jié)構(gòu)網(wǎng)格或非結(jié)構(gòu)網(wǎng)格對模擬對象進(jìn)行離散，一般而言，基于非結(jié)構(gòu)網(wǎng)格的模擬方法要比基于結(jié)構(gòu)網(wǎng)格的模擬方法需要更多計(jì)算量（董良國，2005；王雪秋和孫建國，2008；朱多林和白超英，2011；孫小東等，2012）.貼體網(wǎng)格是為適應(yīng)復(fù)雜幾何邊界而出現(xiàn)的一種結(jié)構(gòu)網(wǎng)格，該網(wǎng)格已在計(jì)算流體力學(xué)中廣泛應(yīng)用，目前也在地球物理模擬中推廣應(yīng)用（孫章慶等，2009；孫建國和蔣麗麗，2009；蘭海強(qiáng)等，2011；Lan and Zhong，2011；劉一峰和蘭海強(qiáng)，2012）.貼體網(wǎng)格生成方法主要有代數(shù)法、插值法、保角變換法和偏微分方程法（陶文銓，2001；蔣麗麗和孫建國，2008）.這些方法生成貼體網(wǎng)格時(shí)，必須首先給定離散的邊界條件.目前對面向地質(zhì)條件的貼體網(wǎng)格生成技術(shù)研究不足，尤其在生成網(wǎng)格時(shí)對邊界離散優(yōu)化方法研究不足.邊界的離散優(yōu)化問題可描述為給定節(jié)點(diǎn)量的曲線逼近和曲面網(wǎng)格優(yōu)化問題，本文研究內(nèi)容集中在給定節(jié)點(diǎn)量的曲線逼近和曲面網(wǎng)格優(yōu)化方法，以解決在地質(zhì)條件下的貼體網(wǎng)格生成中的邊界離散問題.

在數(shù)值模擬和可視化技術(shù)中，曲線和曲面分別以折線段和四邊形去逼近.雖然學(xué)者對曲線逼近和曲面網(wǎng)格優(yōu)化做了很多研究，大部分研究都是自由曲線曲面的造型技術(shù)（吳福鳴和潘日晶，2011；甘屹等，2002）.雖然理論相對成熟，但并不適合給定幾何形狀和節(jié)點(diǎn)量的曲線逼近和曲面網(wǎng)格優(yōu)化.文中，我們提出最大長度準(zhǔn)則、最大面積準(zhǔn)則、最小距離修正法和限定網(wǎng)格步長的優(yōu)化算法，以實(shí)現(xiàn)曲線逼近和曲面網(wǎng)格優(yōu)化，為貼體網(wǎng)格生成提供高質(zhì)量的離散邊界.

曲線作為平面和曲面的邊界，在二維貼體網(wǎng)格生成時(shí)必須先根據(jù)節(jié)點(diǎn)量完成曲線離散，離散結(jié)果要盡可能的逼近曲線形狀.最常逼近曲線的方法有樣條法、圓弧法和直線段法（王振祿等，2007；孫彤和童小華，2010；鄭永前等，2011）.雖然貼體網(wǎng)格在物理域上是以曲線作為邊界，但在網(wǎng)格生成時(shí)并不會關(guān)心相鄰節(jié)點(diǎn)間曲線形狀.因此，利用樣條曲線段和圓弧段逼近曲線的方法在貼體網(wǎng)格生成中很少應(yīng)用，目前生成網(wǎng)格時(shí)主要用直線段來逼近曲線.直線段逼近曲線的方法主要有：等間距法、等步長法、等誤差法.等間距法無法自動(dòng)實(shí)現(xiàn)復(fù)雜曲線形狀的逼近，等誤差法無法限定節(jié)點(diǎn)數(shù)量，故只有等步長法適合給定節(jié)點(diǎn)量的曲線逼近.等步長法具有實(shí)施簡單、效率高的特征，但在實(shí)施過程中沒有考慮曲線形狀特征.本文依據(jù)最大長度準(zhǔn)則實(shí)現(xiàn)定節(jié)點(diǎn)量的直線段逼近曲線，根據(jù)該準(zhǔn)則離散優(yōu)化的折線段能很好地反應(yīng)原曲線的幾何形狀，該準(zhǔn)則對機(jī)械加工、數(shù)值模擬、曲線顯示都具有重要意義.

曲面作為三維幾何體的邊界，通常以三角形和四邊形進(jìn)行曲面離散.三維貼體網(wǎng)格邊界以四邊形結(jié)構(gòu)網(wǎng)格形式離散.通常離散方法有插值算法、NURBS曲面逼近、偏微分方程法等（朱心雄，2001；楊旭東等，2002；吁日新和朱自強(qiáng)，2003；梅中義和范玉青，2003）.這些方法均不含對網(wǎng)格的優(yōu)化過程，目前曲面網(wǎng)格主要的優(yōu)化方法有：Laplacian光順方法、等參數(shù)修勻法（田紅等，2009）.Laplacian光順方法將內(nèi)部節(jié)點(diǎn)移至與其相鄰共線結(jié)點(diǎn)坐標(biāo)的平均值處，等參數(shù)修勻是將網(wǎng)格內(nèi)部節(jié)點(diǎn)的坐標(biāo)修改為相鄰共面結(jié)點(diǎn)坐標(biāo)的平均值.這兩種優(yōu)化方法都沒有考慮優(yōu)化對曲面網(wǎng)格形狀特征的影響，文中，我們提出一種新的四邊形曲面結(jié)構(gòu)網(wǎng)格的優(yōu)化算法，最大面積準(zhǔn)則優(yōu)化法.在不改變網(wǎng)格拓?fù)浣Y(jié)構(gòu)和節(jié)點(diǎn)數(shù)量的條件下使得網(wǎng)格能最多的描述幾何體的邊界形狀.最大面積準(zhǔn)則優(yōu)化的網(wǎng)格不適用于數(shù)值計(jì)算，文中提出了限定網(wǎng)格步長的優(yōu)化算法，可優(yōu)化數(shù)值計(jì)算的網(wǎng)格.

最大準(zhǔn)則的優(yōu)化問題求解費(fèi)時(shí).為解決此問題我們提出改進(jìn)型單粒子優(yōu)化算法.單粒子優(yōu)化算法的思想來源于粒子群優(yōu)化算法，粒子群優(yōu)化算法（Particle Swarm Optimization algorithm，PSO），也可以稱為微粒群算法，最早是由美國的心理學(xué)研究者Kennedy博士和從事計(jì)算智能研究的Eberhart于1995年提出的一種基于群智能的隨機(jī)搜索算法（Kennedy and Eberhart，1995）.該算法最初是受到鳥類群體活動(dòng)規(guī)律性的啟發(fā)，進(jìn)而利用群體智能建立的一個(gè)簡化模型.由于該算法概念簡明、實(shí)現(xiàn)方便、收斂速度快、參數(shù)設(shè)置少，是一種高效的搜索算法，近年來受到學(xué)術(shù)界的廣泛重視.但在優(yōu)化復(fù)雜高維多模函數(shù)時(shí)，PSO算法很容易陷入局部最優(yōu)，并出現(xiàn)早熟收斂現(xiàn)象.為了改善算法性能，Eberhart和Kennedy于1995提出了全局粒子群算法和局部粒子群算法（Eberhart and Kennedy，1995）.Shi和Eberhart于1998年對PSO算法的速度項(xiàng)引入了慣性權(quán)重，并提出在進(jìn)化過程中動(dòng)態(tài)調(diào)整慣性權(quán)重以平衡收斂的全局性和收斂速度（Shi and Eberhart，1998）.Angeline于1998年將進(jìn)化計(jì)算中的標(biāo)準(zhǔn)選擇機(jī)制引入到粒子群優(yōu)化算法中，提出了帶選擇機(jī)制的粒子群優(yōu)化算法（Angeline，1998）.Clerc和Kenndy于2002年提出了帶收縮因子的粒子群優(yōu)化算法（Clerc and Kennedy，2002）.Mendes和Kennedy于2004年提出了一種信息充分共享的粒子群優(yōu)化算法（Mendes etal，2004；Kennedy and Mendes，2006）.紀(jì)震于2010年提出了智能單粒子優(yōu)化算法（Intelligent Single Particle Optimizer，ISPO），確保算法在優(yōu)化復(fù)雜的具有大量局部最優(yōu)化點(diǎn)的高維多模函數(shù)方面的優(yōu)勢（紀(jì)震等，2010）.本文在ISPO算法的基礎(chǔ)上，提出了改進(jìn)的單粒子優(yōu)化算法，該算法能高效求解給定節(jié)點(diǎn)量的曲線逼近和曲面優(yōu)化問題.

2 限定節(jié)點(diǎn)量的曲線逼近和曲面網(wǎng)格優(yōu)化

2.1 限定節(jié)點(diǎn)量的曲線逼近

直線逼近空間曲線時(shí)，應(yīng)做到節(jié)點(diǎn)在曲線曲率變化大的位置密集，在曲率變化小的位置疏散，而等步長法在逼近空間曲線時(shí)忽略了該關(guān)鍵點(diǎn).本文利用最大長度準(zhǔn)則解決此問題.最大長度準(zhǔn)則是在給定曲線形狀和節(jié)點(diǎn)量的條件下，直線段長度和最大的離散結(jié)果能更好地描述曲線.

最大長度準(zhǔn)則定義為：設(shè)Q＝｛xi｜i＝1，2，…，N｝為曲線上N個(gè)節(jié)點(diǎn)集合，該曲線的形狀由f（x）給出，L為集合Q組成的折線段長度，若Lmax為任意集合Q組成折線段的長度集合中的最大值，則組成Lmax的集合Q所描述的曲線最逼真.下面通過實(shí)例展示最大長度準(zhǔn)則和等步長法的對比結(jié)果.

我們采用y＝30sin（3πx／100），0≤x≤100的曲線來說明根據(jù)最大長度準(zhǔn)則逼近曲線的性能.圖1至圖4是最大長度準(zhǔn)則法和等步長法的對比結(jié)果.從圖1和圖2可以看出，以最大長度準(zhǔn)則逼近曲線時(shí)，在描述曲線形狀，逼近誤差方面均優(yōu)于等步長法.圖3說明以最大長度準(zhǔn)則逼近曲線時(shí)用少量節(jié)點(diǎn)數(shù)就能接近原曲線的長度，而等間距法需要用大量節(jié)點(diǎn)才能接近原曲線長度，故最大長度準(zhǔn)則逼近曲線的誤差小于等步長法.圖4說明兩種方法的節(jié)點(diǎn)平均斜率跳躍值相近，最大長度準(zhǔn)則結(jié)果略好于

圖1 節(jié)點(diǎn)量為6的曲線逼近效果對比Fig.1 Comparison of six-node curve approximation by two methods

圖2 節(jié)點(diǎn)量為21的曲線逼近效果對比Fig.2 Comparison of twenty-one node curve approximation by two methods

圖3 兩種逼近算法生成直線的長度與原曲線長度的對比Fig.3 Comparison of curve lengths generated by two approximation algorithms and the original curve length

圖4 兩種逼近算法生成節(jié)點(diǎn)的平均斜率跳躍值的對比Fig.4 Comparison of mean curvature jumping values generated by two approximation algorithms

等步長法，平均斜率跳躍值為

2.2 限定節(jié)點(diǎn)量的曲面網(wǎng)格優(yōu)化

曲面網(wǎng)格優(yōu)化方法中，Laplacian光順法和等參數(shù)修勻法主要用于數(shù)值模擬計(jì)算網(wǎng)格的優(yōu)化，這兩種方法均沒有考慮優(yōu)化對曲面網(wǎng)格形狀的影響.本文提出的最大面積優(yōu)化準(zhǔn)則可實(shí)現(xiàn)限定節(jié)點(diǎn)量曲面網(wǎng)格優(yōu)化，使網(wǎng)格描述的曲面形狀更加逼真.

最大面積優(yōu)化準(zhǔn)則定義為：設(shè)Q＝｛xi，j｜i＝1，2，…，N；j＝1，2，…，M｝為曲面上N×M個(gè)節(jié)點(diǎn)集合，該曲面的形狀由f（x）給出，S為集合Q組成四邊形網(wǎng)格的面積，若Smax為任意集合Q組成四邊形網(wǎng)格的面積集合中的最大值，則組成Smax的集合Q所描述的曲面最逼真.

圖5 等參數(shù)修勻法生成原始曲面網(wǎng)格Fig.5 Original surface grid generated by optimization-based smoothing

圖6 Laplacian光順方法優(yōu)化曲面網(wǎng)格的效果圖Fig.6 Surface grid optimized by Laplacian smoothing

圖7 最大面積準(zhǔn)則優(yōu)化曲面網(wǎng)格的效果圖Fig.7 Surface grid optimized by maximum area criteria

本文采用z＝30sin（πx／50）+30cos（πy／50），0≤x≤200，0≤y≤200的曲面來說明最大面積準(zhǔn)則優(yōu)化曲面網(wǎng)格時(shí)的性能.根據(jù)圖5至圖7可以看出，用Laplacian光順法優(yōu)化等參數(shù)修勻法生成的網(wǎng)格后網(wǎng)格步長趨于一致，而用最大面積準(zhǔn)則優(yōu)化的網(wǎng)格在法曲率半徑變化大位置（谷和峰）網(wǎng)格密集.從圖8中可以看出，最大面積準(zhǔn)則優(yōu)化的網(wǎng)格面積最大，Laplacian光順法優(yōu)化的網(wǎng)格面積小于等參數(shù)修勻法生成的網(wǎng)格面積.因此，最大面積準(zhǔn)則優(yōu)化的網(wǎng)格在用于曲面顯示時(shí)效果好于Laplacian光順法和等參數(shù)修勻法，能更好地描述曲面.

圖8 三種網(wǎng)格的面積對比圖Fig.8 Area comparison of three kinds of grids

3 改進(jìn)的單粒子優(yōu)化算法在曲線逼近和曲面網(wǎng)格優(yōu)化中的應(yīng)用

根據(jù)最大準(zhǔn)則將曲線逼近和曲面網(wǎng)格優(yōu)化問題抽象為數(shù)學(xué)優(yōu)化問題，即為給定曲線（曲面）上節(jié)點(diǎn)數(shù)量，尋找一組滿足最大準(zhǔn)則的節(jié)點(diǎn).為求解此優(yōu)化問題，文中提出了改進(jìn)的單粒子優(yōu)化算法.ISPO算法的計(jì)算耗時(shí)主要花費(fèi)在適應(yīng)值計(jì)算上，求解一次適應(yīng)值的耗時(shí)會隨著粒子維數(shù)增加而增加，使得總求解耗時(shí)隨粒子維數(shù)急劇增加（紀(jì)震等，2010）.本文提出的改進(jìn)單粒子優(yōu)化算法吸收了ISPO算法的學(xué)習(xí)策略和多樣性部分，同時(shí)改進(jìn)了適應(yīng)值的計(jì)算.在改進(jìn)的單粒子優(yōu)化算法中以子粒子與其鄰近子粒子的關(guān)系建立局部適應(yīng)值函數(shù)，通過求解局部適應(yīng)值來完成粒子更新，從而達(dá)到提高求解速度的目的.ISPO算法根據(jù)實(shí)際問題的目標(biāo)函數(shù)直接建立適應(yīng)值函數(shù)，而改進(jìn)的單粒子優(yōu)化算法需要根據(jù)實(shí)際問題對目標(biāo)函數(shù)進(jìn)行分析以建立局部適應(yīng)值函數(shù).

3.1 改進(jìn)的單粒子優(yōu)化算法

在改進(jìn)的單粒子優(yōu)化算法中，整個(gè)N維矢量被分為N個(gè)子維矢量，定義為N個(gè)子粒子，即把整個(gè)位置矢量分解成N個(gè)子粒子位置矢量，每一個(gè)子粒子位置矢量與它對應(yīng)的速度矢量被分別表示成~zj和~vj，j＝1，…，N；每個(gè)子粒子用局部適應(yīng)值判斷其影響域內(nèi)局部解的優(yōu)劣性，以循環(huán)迭代更新子粒子位置的方法求得整體粒子的最優(yōu)解.其算法原理描述如下：

（1）子粒子.將N維粒子分解為N個(gè)子粒子，每個(gè)子粒子有自己的位置矢量和速度矢量.如用N-1段線段去逼近一段曲線，曲線產(chǎn)生包含端點(diǎn)在內(nèi)N個(gè)節(jié)點(diǎn).這N個(gè)節(jié)點(diǎn)組成一個(gè)粒子，每一個(gè)節(jié)點(diǎn)為一個(gè)子粒子，即為N個(gè)子粒子.

（2）更新過程.改進(jìn)單粒子優(yōu)化算法的更新過程是每個(gè)子粒子位置矢量按照先后順序從~z1到~zN進(jìn)行循環(huán)求解，用每個(gè)子粒子在當(dāng)前局部條件下的最優(yōu)位置更新相應(yīng)子粒子位置矢量，迭代執(zhí)行D次或者達(dá)到全局最優(yōu)時(shí)結(jié)束迭代.更新子粒子時(shí)速度和位置調(diào)整公式為

算法流程如下：

（2）設(shè)置第一個(gè)計(jì)算子粒子的位置，即j＝0；

（5）若i小于3，則轉(zhuǎn)到（3）計(jì)算；否則判斷j與N的大小；

（6）若j不大于N，則j＝j(luò)+1，L＝0，并轉(zhuǎn)到（3）循環(huán)計(jì)算；否則判斷k與D的大小；

（7）若k大于D，則結(jié)束計(jì)算；否則k＝k+1，并轉(zhuǎn)到（2）循環(huán)計(jì)算.

在子粒子矢量更新過程中，子粒子速度矢量決定著位置矢量，每個(gè)子粒子速度矢量包括兩個(gè)部分：學(xué)習(xí)部分b×L和多樣性部分（a／kp）×r.多樣性部分（a／kp）會隨著迭代次數(shù)的增加而下降，這樣可以使得優(yōu)化的收斂速度加快.學(xué)習(xí)部分b×L是根據(jù)子粒子上一次的速度更新情況智能地調(diào)整子粒子速度矢量，使子粒子速度具有更大的多樣性，避免陷入局部最優(yōu).

3.2 改進(jìn)的單粒子優(yōu)化算法在曲線逼近中的應(yīng)用及性能分析

在曲線逼近時(shí)，最大長度準(zhǔn)則優(yōu)化問題的目標(biāo)函數(shù)描述為

其中N為節(jié)點(diǎn)數(shù)，Li，i+1為節(jié)點(diǎn)i和i+1之間的直線距離.根據(jù)ISPO算法的適應(yīng)值函數(shù)定義，式（5）為ISPO算法的適應(yīng)值函數(shù).從式（5）中不難發(fā)現(xiàn)，當(dāng)節(jié)點(diǎn)i的位置發(fā)生變化時(shí)，只對距離Li-1，i和Li，i+1有直接影響.故曲線逼近優(yōu)化問題的目標(biāo)函數(shù)也可以近似描述為

max（Li-1，i+Li，i+1）（i＝12，…，N-1）， （6）文中以式（6）作為改進(jìn)的單粒子優(yōu)化算法求解曲線逼近優(yōu)化問題的局部適應(yīng)值函數(shù)，每一個(gè)節(jié)點(diǎn)為改進(jìn)的單粒子優(yōu)化算法的子粒子.在文中采用y＝30sin（πx／100），0≤x≤1000函數(shù)曲線分析兩種優(yōu)化算法的性能，分析結(jié)果如圖9所示，在節(jié)點(diǎn)為200時(shí)，以改進(jìn)的單粒子優(yōu)化算法的計(jì)算效率已是ISPO算法的30倍左右.

圖9 兩種算法優(yōu)化耗時(shí)的對比圖Fig.9 Comparison of time-consuming of two optimizers

3.3 改進(jìn)的單粒子優(yōu)化算法在曲面網(wǎng)格中的應(yīng)用

本文借鑒網(wǎng)格能量法的思想，間接實(shí)現(xiàn)最大面積準(zhǔn)則對曲面網(wǎng)格的優(yōu)化，即采用最大長度準(zhǔn)則和最小距離修正法優(yōu)化曲面網(wǎng)格線.最小距離修正法為節(jié)點(diǎn)與相鄰節(jié)點(diǎn)的距離之和最小，以保證網(wǎng)格中不會產(chǎn)生圖10所示的畸形網(wǎng)格，使網(wǎng)格保持光順.

圖10 畸形網(wǎng)格Fig.10 Distorted grid

圖11 修改后的網(wǎng)格Fig.11 Modified grid

從圖10和圖11可以看出，經(jīng)最小距離修正法處理后的網(wǎng)格正交性改善.類比改進(jìn)的單粒子優(yōu)化算法在曲線逼近中的應(yīng)用，可實(shí)現(xiàn)改進(jìn)的單粒子優(yōu)化算法對曲面網(wǎng)格優(yōu)化問題的求解.

最大面積準(zhǔn)則優(yōu)化曲面網(wǎng)格算法流程如下：

（1）根據(jù)曲面網(wǎng)格的特性將網(wǎng)格曲線分為兩族，同一族的曲線互不相交，分別定義為X方向的曲線和Y方向的曲線；

（2）根據(jù)最大長度準(zhǔn)則優(yōu)化Y方向的曲線；

（3）根據(jù)最大長度準(zhǔn)則優(yōu)化X方向的曲線；

（4）根據(jù)最小距離修正法修正曲面網(wǎng)格；

（5）重復(fù)執(zhí)行（2）至（4）步，直到滿足收斂要求.

4 最大準(zhǔn)則優(yōu)化算法及限定網(wǎng)格步長的優(yōu)化算法在地質(zhì)體貼體網(wǎng)格中應(yīng)用

最大準(zhǔn)則優(yōu)化曲面網(wǎng)格能提高網(wǎng)格質(zhì)量，可更好地描述曲線曲面形狀.最大準(zhǔn)則優(yōu)化的網(wǎng)格可更好地描述地質(zhì)體的谷、峰、地勢平坦和坡度變化較小的區(qū)域.圖13所示網(wǎng)格是根據(jù)圖12所示等高線圖生成的三維貼體網(wǎng)格，從圖12和13可以看出：A類標(biāo)記的位置為地勢平坦或者坡度變化較小區(qū)域，網(wǎng)格間距較大；B類標(biāo)記的位置為凹陷的谷底區(qū)域，網(wǎng)格較密；C類標(biāo)記的位置為凸起的峰頂區(qū)域，網(wǎng)格較密.因此，網(wǎng)格在坡度變化較大區(qū)域密集而在坡度變化較小區(qū)域稀疏，能更好地描述出地質(zhì)體特征.

圖12 地質(zhì)體的等高線圖Fig.12 Contour lines of a geology body

圖13 最大準(zhǔn)則優(yōu)化的三維貼體網(wǎng)格邊界Fig.13 Boundaries of a three dimensional body-fitted grid optimized by maximum criteria

數(shù)值計(jì)算中要求網(wǎng)格具有合理的步長范圍和良好的正交性.為將地質(zhì)體網(wǎng)格（圖13）用于數(shù)值計(jì)算，本文提出了限定網(wǎng)格步長的優(yōu)化算法.該算法引入最大步長系數(shù)amax和最小步長系數(shù)amin，使生成網(wǎng)格步長均勻.圖14為限定網(wǎng)格步長的優(yōu)化算法優(yōu)化圖1 3網(wǎng)格的結(jié)果圖（其中amax為1.2，amin為0.8），從圖看出優(yōu)化后的網(wǎng)格步長得到控制，消除了細(xì)條狀和面積過大的網(wǎng)格；圖15為文本算法與Laplacian光順?biāo)惴▋?yōu)化圖13網(wǎng)格后網(wǎng)格線夾角大小分布圖（網(wǎng)格夾角數(shù)為14400），從圖可知，兩種算法優(yōu)化網(wǎng)格后網(wǎng)格線夾角分布范圍一致，而本文算法在85°～95°之間分布高，在其它范圍內(nèi)的分布低，說明本文算法優(yōu)化的網(wǎng)格具有更好地正交性.根據(jù)限定網(wǎng)格步長優(yōu)化算法優(yōu)化逼近曲線的最大長度準(zhǔn)則折線段的算法流程如下：

（1）初始化網(wǎng)格曲線長度最大步長系數(shù)amax和最小步長系數(shù)amin；

（2）初始化節(jié)點(diǎn)標(biāo)號i＝1，并設(shè)置判斷修正結(jié)束因子k＝0；

圖14 限定網(wǎng)格步長優(yōu)化算法優(yōu)化的三維貼體網(wǎng)格邊界Fig.14 Boundaries of a three-dimensional body-fitted grid optimized by the limited grid step length algorithm

圖15 網(wǎng)格線正交性分布圖Fig.15 Orthogonality distribution columns of grid lines

（3）計(jì)算網(wǎng)格曲線的長度，求得步長的平均值Lave；

（6）若i大于1且小于N-1，則跳轉(zhuǎn)到（8）步，否則繼續(xù)往下執(zhí)行；

（7）若i等于N，則優(yōu)化曲線段的節(jié)點(diǎn)下標(biāo)為：j1＝N-2，j2＝N-1，j3＝N，否則j1＝0，j2＝1，j3＝2，并跳轉(zhuǎn)到（12）步；

（14）若i小于節(jié)點(diǎn)數(shù)N，則i＝i+1，并轉(zhuǎn)到（3），否則向下執(zhí)行；

（15）若k等于1，則設(shè)置i＝1和k＝0，并轉(zhuǎn)到（3），否則結(jié)束計(jì)算.

5 結(jié)論

通過對生成貼體網(wǎng)格中邊界離散技術(shù)的研究，提出了最大準(zhǔn)則、限定網(wǎng)格步長的優(yōu)化算法和改進(jìn)的單粒子優(yōu)化算法，使邊界離散技術(shù)得以提高，可生成高質(zhì)量的貼體網(wǎng)格.

本文提出的最大長度準(zhǔn)則逼近曲線的性能高于等步長法，能在確定節(jié)點(diǎn)量的條件下最大限度地描述出原曲線形狀，為二維貼體網(wǎng)格和空間曲面網(wǎng)格生成打好基礎(chǔ).根據(jù)最大面積準(zhǔn)則優(yōu)化的曲面網(wǎng)格能很好地保持原曲面形狀，且網(wǎng)格光順.為求解最大準(zhǔn)則的優(yōu)化問題，文中提出了改進(jìn)的單粒子優(yōu)化算法，該算法求解曲線逼近和曲面網(wǎng)格優(yōu)化問題的效率遠(yuǎn)高于ISPO算法.為解決最大準(zhǔn)則優(yōu)化曲面網(wǎng)格步長無法控制的問題，本文引入限制網(wǎng)格步長的優(yōu)化算法，通過該算法優(yōu)化曲面網(wǎng)格，生成的新網(wǎng)格在保證網(wǎng)格步長不超出要求范圍內(nèi)的同時(shí)具有更好的正交性.通過各種對比分析說明，本文提出的最大長度準(zhǔn)則、最大面積準(zhǔn)則、最小距離修正法、限定網(wǎng)格步長的優(yōu)化算法和改進(jìn)的單粒子優(yōu)化算法性能都高于目前的優(yōu)化算法，為生成高質(zhì)量的貼體網(wǎng)格及其在地球物理模擬中應(yīng)用提供技術(shù)支撐.

本文主要研究了貼體網(wǎng)格生成中邊界的離散技術(shù)，而在生成貼體網(wǎng)格時(shí)沒有考慮對內(nèi)部網(wǎng)格步長控制和內(nèi)部約束條件.受多旋回構(gòu)造運(yùn)動(dòng)的影響，通常地下的地質(zhì)條件復(fù)雜.受巖石巖性和流體特性的影響，地層呈非均質(zhì)性.因此未來工作主要包括含有內(nèi)部約束條件的貼體網(wǎng)格生成技術(shù)及貼體網(wǎng)格在地質(zhì)地球物理中的應(yīng)用.

Angeline P J.1998.Using selection to improve particle swarm optimization.∥IEEE World Congress on Computational Intelligence.USA：Anchorage Alaska，84-98.

Clerc M，Kennedy J.2002.The particle swarm explosion，stability，and convergence in a multidimensional complex space.IEEE Transactions on Evolutionary Computation，6（1）：58-73.

Dong G L.2005.Numerical simulation of seismic wave propagation under complex near surface conditions.Progress in Exploration Geophysics（in Chinese），28（3）：187-194.

Eberhart R C，Kennedy J.1995.A new optimizer using particle swarm theory.∥Proceedings of the Sixth International Symposium on Micro Machine and Human Science，Nagoya，Japan，39-43.

Gan Y，Qi C Q，Chen Y Z.2002.Smoothing of curves and surfaces based on genetic algorithm.Journal of Tongji University（in Chinese），30（3）：322-325.

Ji Z，Zhou J R，Liao H L，etal.2010.A novel intelligent single particle optimizer.Chinese Journal of Computers（in Chinese），33（3）：556-561.

Jiang L L，Sun J G.2008.Source terms of elliptic system in grid generation.Global Geology（in Chinese），27（3）：298-305.

Kennedy J，Eberhart R C.1995.Particle swarm optimization.∥Proceedings of the IEEE international Conference on Neural Networks.Perth，Australia：IEEE，4：1942-1948.

Kennedy J，Mendes R.2006.Neighborhood topologies in fully informed and best-of-neighborhood particle swarms.IEEE Transactions on Systems，Man，and Cybernetics，Part C：Applications and Reviews，36（4）：515-519.

Lan H Q，Liu J，Bai Z M.2011.Wave-field simulation in VTI media with irregular free surface.Chinese J.Geophys.（in Chinese），54（8）：2072-2084.

Lan H Q，Zhong J Z.2011.Three-dimensional wave-field simulation in heterogeneous transversely isotropic medium with irregular free surface.Bulletin of the Seismological Society of America，101（3）：1354-1370.

Liu Y F，Lan H Q.2012.Study on the numerical solutions of the eikonal equation in curvilinear coordinate system.Chinese J.Geophys.（in Chinese），55（6）：2014-2026.

Mei Z Y，F(xiàn)an Y Q.2003.An algorithm for subdividing and approximating NURBS surface into quadrilateral mesh.Journal of Engineering Graphics（in Chinese），24（3）：105-110.

Mendes R，Kennedy J，Neves J.2004.The fully informed particle swarm：simpler，maybe better.IEEE Transactions on Evolutionary Computation，8（3）：204-210.

Shi Y H，Eberhart R C.1998.A modified particle swarm optimizer.∥Proceedings of the IEEE International Conference on Evolutionary Computation.Piscataway，NJ，USA：IEEE，69-73.

Sun J G，Jiang L L.2009.Orthogonal curvilinear grid generation technique used for numeric simulation of geophysical fields in undulating surface condition.Oil Geophysical Prospecting（in Chinese），44（4）：494-500.

Sun T，Tong X H.2010.A comprehensive curve error model based on broken-line approximation in curve position.Engineering of Surveying and Mapping（in Chinese），19（3）：26-30.

Sun X D，Li Z C，Wang X L.2012.Pre-stack reverse-time migration using a finite difference method based on triangular grids.Progress in Geophysics（in Chinese），27（5）：2077-2083.

Sun Z Q，Sun J G，Zhang D L.2009.2DDC Electric field numerical modeling including surface topography using coordinate transformation method.Journal of Jilin University：Earth Science Edition（in Chinese），39（3）：528-534.

Tao W Q.2001.Numerical Heat Transfer（Second Edition）（in Chinese）.Xi′an：Xi′an Jiaotong University Press.

Tian H，Xu N X，Mei G.2009.Method for optimizing hexahedral mesh of object with complex interpolation surfaces.Computer Engineering and Design（in Chinese），30（2）：286-290.

Wang X Q，Sun J G.2008.The state-of-the-art in numerical modeling including surface topography with finite-difference method.Progress in Geophysics（in Chinese），23（1）：40-48.

Wang Z L，Liu P Y，Cai H L.2007.Analysis of an optimizing calculation of no-round curve by beeline approximation.Machinery Design ＆Manufacture（in Chinese），（9）：176-178.

Wu F M，Pan R J.2011.Surface fairing based on energy optimization and nonuniform B-spline wavelets.Computer Engineering and Applications（in Chinese），47（18）：176-178，185.

Yang X D，Qiao Z D，Zhu B.2002.An easily generated and highquality grid generation method for aircraft configuration.Journal of Northwestern Polytechnical University（in Chinese），20（1）：49-53.

Yu R X，Zhu Z Q.2003.Transfinite interpolation grid generation with the application of B spline.Journal of Beijing University of Aeronautics and Astronautics，29（1）：27-30.

Zhong Y Q，Wang Y P，Wang K W.2011.Biarc approximation of freeform curves based on particle swarm optimization algorithm.China Mechanical Engineering（in Chinese），22（21）：2530-2535.

Zhu D L，Bai C Y.2011.Review on the seismic wavefield forward modelling.Progress in Geophysics（in Chinese），26（5）：1588-1599.

Zhu X X.2001.Free curve and surface modeling techniques（in Chinese）.Beijing：Science Press.

附中文參考文獻(xiàn)

董良國.2005.復(fù)雜地表?xiàng)l件下地震波傳播數(shù)值模擬.勘探地球物理進(jìn)展，28（3）：187-194.

甘屹，齊從謙，陳亞洲.2002.基于遺傳算法的曲線曲面光順.同濟(jì)大學(xué)學(xué)報(bào)，30（3）：322-325.

紀(jì)震，周家銳，廖惠連等.2010.智能單粒子優(yōu)化算法.計(jì)算機(jī)學(xué)報(bào)，33（3）：556-561.

蔣麗麗，孫建國.2008.基于Poisson方程的曲網(wǎng)格生成技術(shù).世界地質(zhì)，27（3）：298-305.

蘭海強(qiáng)，劉佳，白志明.2011.VTI介質(zhì)起伏地表地震波場模擬.地球物理學(xué)報(bào)，54（8）：2072-2084.

劉一峰，蘭海強(qiáng).2012.曲線坐標(biāo)系程函方程的求解方法研究.地球物理學(xué)報(bào)，55（6）：2014-2026.

梅中義，范玉青.2003.NURBS曲面的四邊形網(wǎng)格的分割與逼近.工程圖學(xué)學(xué)報(bào)，24（3）：105-110.

孫建國，蔣麗麗.2009.用于起伏地表?xiàng)l件下地球物理場數(shù)值模擬的正交曲網(wǎng)格生成技術(shù).石油地球物理勘探，44（4）：494-500.

孫彤，童小華.2010.基于折線逼近的曲線位置與模型誤差綜合建模.測繪工程，19（3）：26-30.

孫小東，李振春，王小六.2012.三角網(wǎng)格有限差分法疊前逆時(shí)偏移方法研究.地球物理學(xué)進(jìn)展，27（5）：2077-2083.

孫章慶，孫建國，張東良.2009.二維起伏地表?xiàng)l件下坐標(biāo)變換法直流電場數(shù)值模擬.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，39（3）：528-534.

陶文銓.2001.數(shù)值傳熱學(xué)（第二版）.西安：西安交通大學(xué)出版社.

田紅，徐能雄，梅鋼.2009.含復(fù)雜插值曲面實(shí)體六面體網(wǎng)格優(yōu)化方法.計(jì)算機(jī)工程與設(shè)計(jì)，30（2）：286-290.

王雪秋，孫建國.2008.地震波有限差分?jǐn)?shù)值模擬框架下的起伏地表處理方法綜述.地球物理學(xué)進(jìn)展，23（1）：40-48.

王振祿，劉鵬玉，蔡慧林.2007.直線逼近非圓曲線的優(yōu)化算法.機(jī)械設(shè)計(jì)與制造，（9）：176-178.

吳福鳴，潘日晶.2011.基于能量優(yōu)化與非均勻B樣條小波的曲面光順.計(jì)算機(jī)工程與應(yīng)用，47（18）：176-178，185.

楊旭東，喬志德，朱兵.2002.一種改進(jìn)的四邊界插值網(wǎng)格生成方法及其應(yīng)用.西北工業(yè)大學(xué)學(xué)報(bào)，20（1）：49-53.

吁日新，朱自強(qiáng).2003.B樣條應(yīng)用于超限插值的網(wǎng)格生成方法.北京航空航天大學(xué)學(xué)報(bào)，29（1）：27-30.

鄭永前，王云鵬，王科委.2011.基于粒子群優(yōu)化算法的自由曲線雙圓弧逼近.中國機(jī)械工程，22（21）：2530-2535.

朱多林，白超英.2011.基于波動(dòng)方程理論的地震波場數(shù)值模擬方法綜述.地球物理學(xué)進(jìn)展，26（5）：1588-1599.

朱心雄.2001.自由曲線曲面造型技術(shù).北京：科學(xué)出版社.