2014, Vol. 57
2. 中国科学院大学, 北京 100049;
3. 中国石油冀东油田公司, 唐山 063004
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. PetroChina Jidong Oilfield Company, Tangshan 063004, China
目前我国油气勘探已转移到西部地区，西部地区剧烈起伏的地形给地震勘探提出了严峻的挑战，经典反射地震勘探技术在复杂地表地区不再适用，地震波数值模拟中的网格法由于具有对地形和介质灵活划分的特点成为解决此类问题最有效的方法.网格法以结构网格或非结构网格对模拟对象进行离散，一般而言，基于非结构网格的模拟方法要比基于 结构网格的模拟方法需要更多计算量(董良国，2005； 王雪秋和孙建国，2008；朱多林和白超英，2011；孙小东等，2012).贴体网格是为适应复杂几何边界而出现的一种结构网格，该网格已在计算流体力学中广泛应用，目前也在地球物理模拟中推广应用(孙章庆等，2009；孙建国和蒋丽丽，2009；兰海强等，2011；Lan and Zhong， 2011；刘一峰和兰海强，2012).贴体网格生成方法主要有代数法、插值法、保角变换法和偏微分方程法(陶文铨，2001；蒋丽丽和孙建国，2008).这些方法生成贴体网格时，必须首先给定离散的边界条件.目前对面向地质条件的贴体网格生成技术研究不足，尤其在生成网格时对边界离散优化方法研究不足.边界的离散优化问题可描述为给定节点量的曲线逼近和曲面网格优化问题，本文研究内容集中在给定节点量的曲线逼近和曲面网格优化方法，以解决在地质条件下的贴体网格生成中的边界离散问题.
在数值模拟和可视化技术中，曲线和曲面分别以折线段和四边形去逼近.虽然学者对曲线逼近和曲面网格优化做了很多研究，大部分研究都是自由曲线曲面的造型技术(吴福鸣和潘日晶，2011；甘屹等，2002).虽然理论相对成熟，但并不适合给定几何形状和节点量的曲线逼近和曲面网格优化.文中，我们提出最大长度准则、最大面积准则、最小距离修正法和限定网格步长的优化算法，以实现曲线逼近和曲面网格优化，为贴体网格生成提供高质量的离散边界.
曲线作为平面和曲面的边界，在二维贴体网格生成时必须先根据节点量完成曲线离散，离散结果要尽可能的逼近曲线形状.最常逼近曲线的方法有样条法、圆弧法和直线段法(王振禄等，2007；孙彤和童小华，2010；郑永前等，2011).虽然贴体网格在物理域上是以曲线作为边界，但在网格生成时并不会关心相邻节点间曲线形状.因此，利用样条曲线段和圆弧段逼近曲线的方法在贴体网格生成中很少应用，目前生成网格时主要用直线段来逼近曲线.直线段逼近曲线的方法主要有：等间距法、等步长法、等误差法.等间距法无法自动实现复杂曲线形状的逼近，等误差法无法限定节点数量，故只有等步长法适合给定节点量的曲线逼近.等步长法具有实施简单、效率高的特征，但在实施过程中没有考虑曲线形状特征.本文依据最大长度准则实现定节点量的直线段逼近曲线，根据该准则离散优化的折线段能很好地反应原曲线的几何形状，该准则对机械加工、数值模拟、曲线显示都具有重要意义.
曲面作为三维几何体的边界，通常以三角形和四边形进行曲面离散.三维贴体网格边界以四边形结构网格形式离散.通常离散方法有插值算法、NURBS曲面逼近、偏微分方程法等(朱心雄，2001；杨旭东等，2002；吁日新和朱自强，2003；梅中义和范玉青，2003).这些方法均不含对网格的优化过程，目前曲面网格主要的优化方法有：Laplacian光顺方法、等参数修匀法(田红等，2009).Laplacian光顺方法将内部节点移至与其相邻共线结点坐标的平均值处，等参数修匀是将网格内部节点的坐标修改为相邻共面结点坐标的平均值.这两种优化方法都没有考虑优化对曲面网格形状特征的影响，文中，我们提出一种新的四边形曲面结构网格的优化算法，最大面积准则优化法.在不改变网格拓扑结构和节点数量的条件下使得网格能最多的描述几何体的边界形状.最大面积准则优化的网格不适用于数值计算，文中提出了限定网格步长的优化算法，可优化数值计算的网格.
最大准则的优化问题求解费时.为解决此问题我们提出改进型单粒子优化算法.单粒子优化算法的思想来源于粒子群优化算法，粒子群优化算法(Particle Swarm Optimization algorithm，PSO)，也可以称为微粒群算法，最早是由美国的心理学研究者Kennedy博士和从事计算智能研究的Eberhart于1995年提出的一种基于群智能的随机搜索算法(Kennedy and Eberhart， 1995).该算法最初是受到鸟类群体活动规律性的启发，进而利用群体智能建立的一个简化模型.由于该算法概念简明、实现方便、收敛速度快、参数设置少，是一种高效的搜索算法，近年来受到学术界的广泛重视.但在优化复杂高维多模函数时，PSO算法很容易陷入局部最优，并出现早熟收敛现象.为了改善算法性能，Eberhart和Kennedy于1995提出了全局粒子群算法和局部粒子群算法(Eberhart and Kennedy， 1995).Shi和Eberhart于1998年对PSO算法的速度项引入了惯性权重，并提出在进化过程中动态调整惯性权重以平衡收敛的全局性和收敛速度(Shi and Eberhart， 1998).Angeline于1998年将进化计算中的标准选择机制引入到粒子群优化算法中，提出了带选择机制的粒子群优化算法(Angeline，1998).Clerc和Kenndy于2002年提出了带收缩因子的粒子群优化 算法(Clerc and Kennedy， 2002).Mendes和Kennedy 于2004年提出了一种信息充分共享的粒子群优化算法(Mendes et al，2004； Kennedy and Mendes， 2006).纪震于2010年提出了智能单粒子优化算法(Intelligent Single Particle Optimizer，ISPO)，确保算法在优化复杂的具有大量局部最优化点的高维多模函数方面的优势(纪震等，2010).本文在ISPO算法的基础上，提出了改进的单粒子优化算法，该算法能高效求解给定节点量的曲线逼近和曲面优化问题. 2 限定节点量的曲线逼近和曲面网格优化 2.1 限定节点量的曲线逼近
直线逼近空间曲线时，应做到节点在曲线曲率变化大的位置密集，在曲率变化小的位置疏散，而等步长法在逼近空间曲线时忽略了该关键点.本文利用最大长度准则解决此问题.最大长度准则是在给定曲线形状和节点量的条件下，直线段长度和最大的离散结果能更好地描述曲线.
最大长度准则定义为：设Q={xi|i=1，2，…，N}为曲线上N个节点集合，该曲线的形状由f(x)给出，L为集合Q组成的折线段长度，若Lmax为任意集合Q组成折线段的长度集合中的最大值，则组成Lmax的集合Q所描述的曲线最逼真.下面通过实例展示最大长度准则和等步长法的对比结果.
我们采用y=30sin(3πx/100)，0≤x≤100的曲线来说明根据最大长度准则逼近曲线的性能.图 1至图 4是最大长度准则法和等步长法的对比结果.从图 1和图 2可以看出，以最大长度准则逼近曲线时，在描述曲线形状，逼近误差方面均优于等步长法.图 3说明以最大长度准则逼近曲线时用少量节点数就能接近原曲线的长度，而等间距法需要用大量节点才能接近原曲线长度，故最大长度准则逼近曲线的误差小于等步长法.图 4说明两种方法的节 点平均斜率跳跃值相近，最大长度准则结果略好于等步长法，平均斜率跳跃值为
为平均斜率跳跃值，n为直线段数量，i为直线段编号，ki为i直线段斜率.因此最大长度准则逼近曲线的性能优于等步长法，可以更好描述曲线形状，且逼近误差小.
 | 图 1 节点量为6的曲线逼近效果对比Fig. 1 Comparison of six-node curve approximation by two methods |
 | 图 2 节点量为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 两种逼近算法生成节点的平均斜率跳跃值的对比Fig. 4 Comparison of mean curvature jumping values generated by two approximation algorithms |
曲面网格优化方法中，Laplacian光顺法和等参数修匀法主要用于数值模拟计算网格的优化，这两种方法均没有考虑优化对曲面网格形状的影响.本文提出的最大面积优化准则可实现限定节点量曲面网格优化，使网格描述的曲面形状更加逼真.
最大面积优化准则定义为：设Q={xi，j|i=1，2，…，N；j=1，2，…，M}为曲面上N×M个节点集合，该曲面的形状由f(x)给出，S为集合Q组成四边形网格的面积，若Smax为任意集合Q组成四边形网格的面积集合中的最大值，则组成Smax的集合Q所描述的曲面最逼真.
本文采用z=30sin(πx/50)+30cos(πy/50)，0≤x≤200，0≤y≤200的曲面来说明最大面积 准则优化曲面网格时的性能.根据图 5至图 7可以看出，用Laplacian光顺法优化等参数修匀法生成的网格后网格步长趋于一致，而用最大面积准则优化的网格在法曲率半径变化大位置(谷和峰)网格密集.从图 8中可以看出，最大面积准则优化的网格面积最大，Laplacian光顺法优化的网格面积小于等参数修匀法生成的网格面积.因此，最大面积准则优化的网格在用于曲面显示时效果好于Laplacian光顺法和等参数修匀法，能更好地描述曲面.
 | 图 5 等参数修匀法生成原始曲面网格Fig. 5 Original surface grid generated by optimization-based smoothing |
 | 图 6 Laplacian光顺方法优化曲面网格的效果图Fig. 6 Surface grid optimized by Laplacian smoothing |
 | 图 7 最大面积准则优化曲面网格的效果图Fig. 7 Surface grid optimized by maximum area criteria |
 | 图 8 三种网格的面积对比图Fig. 8 Area comparison of three kinds of grids |
根据最大准则将曲线逼近和曲面网格优化问题抽象为数学优化问题，即为给定曲线(曲面)上节点数量，寻找一组满足最大准则的节点.为求解此优化问题，文中提出了改进的单粒子优化算法.ISPO算法的计算耗时主要花费在适应值计算上，求解一次适应值的耗时会随着粒子维数增加而增加，使得总求解耗时随粒子维数急剧增加(纪震等，2010).本文提出的改进单粒子优化算法吸收了ISPO算法的学习策略和多样性部分，同时改进了适应值的计算.在改进的单粒子优化算法中以子粒子与其邻近子粒子的关系建立局部适应值函数，通过求解局部适应值来完成粒子更新，从而达到提高求解速度的目的.ISPO算法根据实际问题的目标函数直接建立适应值函数，而改进的单粒子优化算法需要根据实际问题对目标函数进行分析以建立局部适应值函数.
3.1 改进的单粒子优化算法在改进的单粒子优化算法中，整个N维矢量被分为N个子维矢量，定义为N个子粒子，即把整个位置矢量分解成N个子粒子位置矢量，每一个子粒子位置矢量与它对应的速度矢量被分别表示成 j和 j，j=1，…，N；每个子粒子用局部适应值判断其影响域内局部解的优劣性，以循环迭代更新子粒子位置的方法求得整体粒子的最优解.其算法原理描述如下：
(1)子粒子.将N维粒子分解为N个子粒子，每个子粒子有自己的位置矢量和速度矢量.如用N-1段线段去逼近一段曲线，曲线产生包含端点在内N个节点.这N个节点组成一个粒子，每一个节点为一个子粒子，即为N个子粒子.
(2)更新过程.改进单粒子优化算法的更新过程是每个子粒子位置矢量按照先后顺序从 1到 N进行循环求解，用每个子粒子在当前局部条件下的最优位置更新相应子粒子位置矢量，迭代执行D次或者达到全局最优时结束迭代.更新子粒子时速度和位置调整公式为
参数L表示学习矢量，r 表示在区间[-0.5，0.5]上均匀分布的随机矢量； a，p，s和b分别表示多样性因子、下降因子、搜索因子和加速因子，其取值与具体求解问题有关，文中a=1，p=0.5，s=2和b=1；f()为局部适应值函数，文中为子粒子与相邻子粒子的距离之和.算法流程如下：
(1)初始化子粒子的位置矢量z 0j、a、p、s和b参数，并置k=1；
(2)设置第一个计算子粒子的位置，即j=0；
(3)计算Vkj、zjk、f(zkj)和f(z k-1j)的值，并比较f(z kj)与f(z k-1j)的大小；
(4)若f(z kj)大于f(z k-1j)，则 kj=k-1j+kj，L = kj，i=0，并转到(3)循环计算；否则kj=k-1j，L = L /s，i=i+1，并向下执行；
(5)若i小于3，则转到(3)计算；否则判断j与N的大小；
(6)若j不大于N，则j=j+1，L =0，并转到(3)循环计算；否则判断k与D的大小；
(7)若k大于D，则结束计算；否则k=k+1，并转到(2)循环计算.
在子粒子矢量更新过程中，子粒子速度矢量决定着位置矢量，每个子粒子速度矢量包括两个部分：学习部分b× L和多样性部分(a/kp)× r .多样性部分(a/kp)会随着迭代次数的增加而下降，这样可以使得优化的收敛速度加快.学习部分b× L是根据子粒子上一次的速度更新情况智能地调整子粒子速度矢量，使子粒子速度具有更大的多样性，避免陷入局部最优.
3.2 改进的单粒子优化算法在曲线逼近中的应用及性能分析在曲线逼近时，最大长度准则优化问题的目标函数描述为
 | 图 9 两种算法优化耗时的对比图Fig. 9 Comparison of time-consuming of two optimizers |
本文借鉴网格能量法的思想，间接实现最大面积准则对曲面网格的优化，即采用最大长度准则和最小距离修正法优化曲面网格线.最小距离修正法为节点与相邻节点的距离之和最小，以保证网格中不会产生图 10所示的畸形网格，使网格保持光顺.
 | 图 10 畸形网格Fig. 10 Distorted grid |
从图 10和图 11可以看出，经最小距离修正法处理后的网格正交性改善.类比改进的单粒子优化算法在曲线逼近中的应用，可实现改进的单粒子优化算法对曲面网格优化问题的求解.
 | 图 11 修改后的网格Fig. 11 Modified grid |
最大面积准则优化曲面网格算法流程如下：
(1)根据曲面网格的特性将网格曲线分为两族，同一族的曲线互不相交，分别定义为X方向的曲线和Y方向的曲线；
(2)根据最大长度准则优化Y方向的曲线；
(3)根据最大长度准则优化X方向的曲线；
(4)根据最小距离修正法修正曲面网格；
(5)重复执行(2)至(4)步，直到满足收敛要求.
4 最大准则优化算法及限定网格步长的优化算法在地质体贴体网格中应用最大准则优化曲面网格能提高网格质量，可更好地描述曲线曲面形状.最大准则优化的网格可更好地描述地质体的谷、峰、地势平坦和坡度变化较小的区域.图 13所示网格是根据图 12所示等高线图生成的三维贴体网格，从图 12和13可以看出：A类标记的位置为地势平坦或者坡度变化较小区域，网格间距较大；B类标记的位置为凹陷的谷底区域，网格较密；C类标记的位置为凸起的峰顶区域，网格较密.因此，网格在坡度变化较大区域密集而在坡度变化较小区域稀疏，能更好地描述出地质体特征.
 | 图 12 地质体的等高线图Fig. 12 Contour lines of a geology body |
 | 图 13 最大准则优化的三维贴体网格边界Fig. 13 Boundaries of a three dimensional body-fitted grid optimized by maximum criteria |
数值计算中要求网格具有合理的步长范围和良好的正交性.为将地质体网格(图 13)用于数值计算，本文提出了限定网格步长的优化算法.该算法引入最大步长系数amax和最小步长系数amin，使生成网格步长均匀.图 14为限定网格步长的优化算法优化 图 13网格的结果图(其中amax为1.2，amin为0.8)，
 | 图 14 限定网格步长优化算法优化的三维贴体网格边界Fig. 14 Boundaries of a three-dimensional body-fitted grid optimized by the limited grid step length algorithm |
从图看出优化后的网格步长得到控制，消除了细条 状和面积过大的网格；图 15为文本算法与Laplacian 光顺算法优化图 13网格后网格线夹角大小分布图(网格夹角数为14400)，从图可知，两种算法优化网 格后网格线夹角分布范围一致，而本文算法在85°~95° 之间分布高，在其它范围内的分布低，说明本文算法优化的网格具有更好地正交性.根据限定网格步长优化算法优化逼近曲线的最大长度准则折线段的算法流程如下：
(1)初始化网格曲线长度最大步长系数amax和最小步长系数amin；
(2)初始化节点标号i=1，并设置判断修正结束因子k=0；
(3)计算网格曲线的长度，求得步长的平均值Lave；
(4)计算最大步长Lmax=amax*Lave和最小步长Lmin=amin*Lave；
(5)计算节点i-1和i的直线距离Li-1，i；
(6)若i大于1且小于N-1，则跳转到(8)步，否则继续往下执行；
(7)若i等于N，则优化曲线段的节点下标为：j1=N-2，j2=N-1，j3=N，否则j1=0，j2=1，j3=2，并跳转到(12)步；
(8)计算直线距离Li-2，i-1和Li，i+1；
(9)若Li-2，i-1大于Li，i+1则继续往下执行，否则跳到(11)步；
(10)若Li-1，i小于Lmin，则j1=i-2，j2=i-1，j3=i，否则j1=i-1，j2=i，j3=i+1，并跳转到(12)步；
(11)若Li-1，i大于Lmax，则j1=i-2，j2=i-1，j3=i，否则j1=i-1，j2=i，j3=i+1；
(12)若Li-1，i小于Lmin，则在曲线段j1j3上寻找节点j2的坐标，使得Li-1，i等于Lmin，并设置k=1；
(13)若Li-1，i大于Lmax，则在曲线段j1j3上寻找节点j2的坐标，使得Li-1，i等于Lmax，并设置k=1；
(14)若i小于节点数N，则i=i+1，并转到(3)，否则向下执行；
(15)若k等于1，则设置i=1和k=0，并转到(3)，否则结束计算.
 | 图 15 网格线正交性分布图Fig. 15 Orthogonality distribution columns of grid lines |
通过对生成贴体网格中边界离散技术的研究，提出了最大准则、限定网格步长的优化算法和改进的单粒子优化算法，使边界离散技术得以提高，可生成高质量的贴体网格.
本文提出的最大长度准则逼近曲线的性能高于等步长法，能在确定节点量的条件下最大限度地描述出原曲线形状，为二维贴体网格和空间曲面网格生成打好基础.根据最大面积准则优化的曲面网格能很好地保持原曲面形状，且网格光顺.为求解最大准则的优化问题，文中提出了改进的单粒子优化算法，该算法求解曲线逼近和曲面网格优化问题的效率远高于ISPO算法.为解决最大准则优化曲面网格步长无法控制的问题，本文引入限制网格步长的优化算法，通过该算法优化曲面网格，生成的新网格在保证网格步长不超出要求范围内的同时具有更好的正交性.通过各种对比分析说明，本文提出的最大长度准则、最大面积准则、最小距离修正法、限定网格步长的优化算法和改进的单粒子优化算法性能都高于目前的优化算法，为生成高质量的贴体网格及其在地球物理模拟中应用提供技术支撑.
本文主要研究了贴体网格生成中边界的离散技术，而在生成贴体网格时没有考虑对内部网格步长控制和内部约束条件.受多旋回构造运动的影响，通常地下的地质条件复杂.受岩石岩性和流体特性的影响，地层呈非均质性.因此未来工作主要包括含有内部约束条件的贴体网格生成技术及贴体网格在地质地球物理中的应用.
| [1] | Angeline P J. 1998. Using selection to improve particle swarm optimization. //IEEE World Congress on Computational Intelligence. USA: Anchorage Alaska, 84-98. |
| [2] | 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. |
| [3] | 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. |
| [4] | 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. |
| [5] | 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. |
| [6] | Ji Z, Zhou J R, Liao H L, et al. 2010. A novel intelligent single particle optimizer. Chinese Journal of Computers (in Chinese), 33(3): 556-561. |
| [7] | Jiang L L, Sun J G. 2008. Source terms of elliptic system in grid generation. Global Geology (in Chinese), 27(3): 298-305. |
| [8] | Kennedy J, Eberhart R C. 1995. Particle swarm optimization. //Proceedings of the IEEE international Conference on Neural Networks. Perth, Australia: IEEE, 4: 1942-1948. |
| [9] | 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. |
| [10] | 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. |
| [11] | 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. |
| [12] | 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. |
| [13] | Mei Z Y, Fan Y Q. 2003. An algorithm for subdividing and approximating NURBS surface into quadrilateral mesh. Journal of Engineering Graphics (in Chinese), 24(3): 105-110. |
| [14] | Mendes R, Kennedy J, Neves J. 2004. The fully informed particle swarm: simpler, maybe better. IEEE Transactions on Evolutionary Computation, 8(3): 204-210. |
| [15] | 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. |
| [16] | 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. |
| [17] | 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. |
| [18] | 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. |
| [19] | Sun Z Q, Sun J G, Zhang D L. 2009. 2D DC Electric field numerical modeling including surface topography using coordinate transformation method. Journal of Jilin University: Earth Science Edition (in Chinese), 39(3): 528-534. |
| [20] | Tao W Q. 2001. Numerical Heat Transfer (Second Edition) (in Chinese). Xi'an: Xi'an Jiaotong University Press. |
| [21] | 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. |
| [22] | 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. |
| [23] | 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. |
| [24] | 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. |
| [25] | Yang X D, Qiao Z D, Zhu B. 2002. An easily generated and high-quality grid generation method for aircraft configuration. Journal of Northwestern Polytechnical University (in Chinese), 20(1): 49-53. |
| [26] | 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. |
| [27] | 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. |
| [28] | Zhu D L, Bai C Y. 2011. Review on the seismic wavefield forward modelling. Progress in Geophysics (in Chinese), 26(5): 1588-1599. |
| [29] | Zhu X X. 2001. Free curve and surface modeling techniques (in Chinese). Beijing: Science Press. |
| [30] | 董良国. 2005. 复杂地表条件下地震波传播数值模拟. 勘探地球物理进展, 28(3): 187-194 . |
| [31] | 甘屹, 齐从谦, 陈亚洲. 2002. 基于遗传算法的曲线曲面光顺. 同济大学学报, 30(3): 322-325 . |
| [32] | 纪震, 周家锐, 廖惠连等. 2010. 智能单粒子优化算法. 计算机学报, 33(3): 556-561 . |
| [33] | 蒋丽丽, 孙建国. 2008. 基于Poisson方程的曲网格生成技术. 世界地质, 27(3): 298-305 . |
| [34] | 兰海强, 刘佳, 白志明. 2011. VTI介质起伏地表地震波场模拟. 地球物理学报, 54(8): 2072-2084 . |
| [35] | 刘一峰, 兰海强. 2012. 曲线坐标系程函方程的求解方法研究. 地球物理学报, 55(6): 2014-2026 . |
| [36] | 梅中义, 范玉青. 2003. NURBS曲面的四边形网格的分割与逼近. 工程图学学报, 24(3): 105-110 . |
| [37] | 孙建国, 蒋丽丽. 2009. 用于起伏地表条件下地球物理场数值模拟的正交曲网格生成技术. 石油地球物理勘探, 44(4): 494-500 . |
| [38] | 孙彤, 童小华. 2010. 基于折线逼近的曲线位置与模型误差综合建模. 测绘工程, 19(3): 26-30 . |
| [39] | 孙小东, 李振春, 王小六. 2012. 三角网格有限差分法叠前逆时偏移方法研究. 地球物理学进展, 27(5): 2077-2083 . |
| [40] | 孙章庆, 孙建国, 张东良. 2009. 二维起伏地表条件下坐标变换法直流电场数值模拟. 吉林大学学报: 地球科学版, 39(3): 528-534 . |
| [41] | 陶文铨. 2001. 数值传热学 (第二版). 西安: 西安交通大学出版社. |
| [42] | 田红, 徐能雄, 梅钢. 2009. 含复杂插值曲面实体六面体网格优化方法. 计算机工程与设计, 30(2): 286-290 . |
| [43] | 王雪秋, 孙建国. 2008. 地震波有限差分数值模拟框架下的起伏地表处理方法综述. 地球物理学进展, 23(1): 40-48 . |
| [44] | 王振禄, 刘鹏玉, 蔡慧林. 2007. 直线逼近非圆曲线的优化算法. 机械设计与制造, (9): 176-178 . |
| [45] | 吴福鸣, 潘日晶. 2011. 基于能量优化与非均匀B样条小波的曲面光顺. 计算机工程与应用, 47(18): 176-178, 185 . |
| [46] | 杨旭东, 乔志德, 朱兵. 2002. 一种改进的四边界插值网格生成方法及其应用. 西北工业大学学报, 20(1): 49-53 . |
| [47] | 吁日新, 朱自强. 2003. B样条应用于超限插值的网格生成方法. 北京航空航天大学学报, 29(1): 27-30 . |
| [48] | 郑永前, 王云鹏, 王科委. 2011. 基于粒子群优化算法的自由曲线双圆弧逼近. 中国机械工程, 22(21): 2530-2535 . |
| [49] | 朱多林, 白超英. 2011. 基于波动方程理论的地震波场数值模拟方法综述. 地球物理学进展, 26(5): 1588-1599 . |
| [50] | 朱心雄. 2001. 自由曲线曲面造型技术. 北京: 科学出版社. |