RPIM求解点源二维变分问题的最优形状参数

引用本文

李俊杰, 严家斌. RPIM求解点源二维变分问题的最优形状参数[J]. 物探与化探, 2017,39(6): 1233-1237.
LI Jun-Jie, YAN Jia-Bin. Optimal shape parameters of RPIM for resolving point source two-dimensional variational problem[J]. Geophysical and Geochemical Exploration, 2017,39(6): 1233-1237. 复制到剪切板

Doi:10.11720/wtyht.2015.6.22
Permissions

《物探与化探》编辑部

RPIM求解点源二维变分问题的最优形状参数

李俊杰¹, 严家斌²

1.浙江省水利水电勘测设计院,浙江杭州 310002

2.中南大学地球科学与信息物理学院有色资源与地质灾害探查湖南省重点实验室,湖南长沙 410083

作者简介: 李俊杰(1989-),江西上饶人,硕士,助理工程师,从事大地电磁无网格化正演研究。Email:lijunjiecsu@163.com

收稿日期: 2014-12-09

基金: 国家自然科学基金资助项目(40874055); 湖南省自然科学基金资助项目(14JJ2012)

摘要

径向基点插值法(RPIM)作为一种高精度的无网格方法,其形函数采用与径向基函数结合的插值方法构造,边界条件可直接加载。将RPIM用于点源二维变分问题的求解,介绍了RPIM的近似原理;推导了点源二维问题的RPIM总体矩阵表达式,简述了背景网格积分技术,研究了高斯点数目对RPIM计算精度的影响;最后通过数值试验得出了支持域无量纲尺寸 α最优选择区间与RPIM形状参数最优值。研究结果表明:RPIM求解点源二维变分问题具有较好的鲁棒性, α最优区间为1.0~1.2。

关键词: 径向基点插值法; 点源; 径向基函数; 点源二维变分问题

中图分类号:P631 文献标志码:A 文章编号:1000-8918(2015)06-1233-05

Optimal shape parameters of RPIM for resolving point source two-dimensional variational problem

LI Jun-Jie¹, YAN Jia-Bin²

1. Zhejiang Design Institute of Water Conservancy and Hydroelectric Power, Hangzhou 310002, China

2. Key Laboratory of Nonferrous Resources and Geological Hazard Detection, School of Geosciences and Info-Physics, Central South University, Changsha 410083, China

Abstract

Radial point interpolation method (RPIM) is a kind of high precision meshfree method. As its shape function is constructed by interpolation method in combination with radial basis function, the boundary conditions can be directly loaded. This paper utilizes RPIM to the calculation of point source two-dimensional electric field. Firstly, the approximate principle of RPIM is introduced in detail and the discrete system matrix expression is deduced corresponding to point source two-dimensional variational problem. Secondly, background grid integral technology is briefly introduced and the influence of different number of gauss points on calculation accuracy of RPIM is discussed. Lastly, the optimal range of support domain dimensionless size and the shape parameter optimal value of RPIM are obtained through numerical experiments. Studies show that RPIM has robustness for solving point source two-dimensional variational problem, and the optimal α range is 1.0 to 1.2.

Keyword: radial point interpolation method; point source; radial basis function; point-source two-dimensional variational problem

Show Figures

有限元法(finite element method , FEM $\begin{matrix} )^{[18]} \end{matrix}$ 计算复杂介质点源二维问题需采用非结构化三角网格剖分计算区域, 此种网格的生成较困难。无网格法^[9]作为FEM的补充和发展, 是近十几年来兴起于计算力学领域的一类新型数值算法, 其物性加载在只与坐标位置有关的高斯点上, 故节点规则分布下也能适应复杂模型的计算。无单元Galerkin法(element-free galerkin method, EFGM)与点插值法(point interpolation method, PIM)作为两种较成熟的全域弱式无网格方法, 已被成功应用于地震波场^[10]、雷达波场^[11]与大地电磁场响应 $\begin{matrix} ^{[1214]} \end{matrix}$ 的计算中。EFGM的缺点在于形函数通过拟合方法构造, 不满足Kronecker delta函数(N_i(X)=δ _ij)的性质, 边界条件加载不便; PIM形函数的构建在节点分布不合理时则易出现奇异性问题。

径向基点插值法(Radial point interpolation method, RPIM)相较PIM引入了径向基函数, 改善了PIM插值函数的奇异性, 且该算法支持域无量纲尺寸的选择区间更大, 从而能更好地处理各类工程与科学计算问题。文中以点源二维电场问题为例介绍了RPIM的近似原理, 从点源二维变分问题出发结合Galerkin法推导了总体矩阵的无网格化离散表达式, 最后通过数值计算得出了适用于点源二维正演问题的RPIM支持域无量纲尺寸、高斯点数量及径向基函数形状参数的最优值。

1 点源二维电场波数域变分问题

点源为三维源, 需先采用傅里叶变换将三维点源从空间域转化为波数域二维源, 当地下电性结构为二维时, 取走向为z轴, x轴与z轴垂直, y轴垂直向上, 求解域为矩形区域, 4个顶点依次以A、B、C、D逆时针编号, 波数域电位满足变分问题^[15]

式中: 是二维哈密顿算子, = e_x+ e_y; σ 为电导率; k为波数; I为供电电流; δ 为单位脉冲函数; P为点源的位置; r为点源到边界Γ 上的距离; K₀与K₁分别为第二类零阶与一阶修正Bessel函数; cos(r, n)为矢径r与边界外法向n的夹角的余弦。

2 变分问题的无网格解法

RPIM利用位于积分点支持域内的场节点构造形函数。支持域是人为定义的一个区域, 其形状主要有圆形与矩形两种, 对于节点均匀分布的情况, 支持域尺寸一般取横、纵向节点间距与支持域无量纲尺寸α 之积。α 用于控制实际支持域的大小, 是对无网格法计算精度影响很大的重要参数^[13]。图1为RPIM背景网格、支持域、积分点与场节点示意, 由于背景网格的积分常选用高斯积分法, 故积分点又称高斯积分点或高斯点。

	Figure Option View Download New Window
	图1 RPIM背景网格、支持域、积分点与场节点示意

2.1 径向基点插值近似

求解域Ω 中任意一点X处的场变量U(X)的径向基点插值近似表达式为

式中:R^T(X)为径向基函数(radial basis function, RBF), p(X)为二维空间坐标X^T=[x, y]的基函数; a与b是系数向量; n为RBF的个数, m为单项式的个数。常用RBF主要有MQ(multi-quadrics)函数、高斯指数(EXP)函数、薄板样条(thin plate spline, TPS)函数:

$\begin{matrix} \begin{matrix} MQ & R_{i} (x, y) = [r_{i}^{2} + (α_{c} d_{c})^{2}]^{q} \\ EXP & R_{i} (x, y) = \exp [- α_{c} {(\frac{r_{i}}{d_{c}})}^{2}] \\ TPS & R_{i} (x, y) = r_{i}^{η} \end{matrix}\} (3) \end{matrix}$

式中:d_c为与节点间距有关的特征长度, 当节点均匀分布时可取d_c= $\begin{matrix} \sqrt{d_{cx}^{2} + d_{cy}^{2}} \end{matrix}$ ; r_i= $\begin{matrix} \sqrt{(x - x_{i})^{2} + (y - y_{i})^{2}} \end{matrix}$ ; x与y为高斯点位置的坐标, x_i与y_i为支持域内节点的坐标; α _c、q与η 为径向基函数的形状参数, 在不同学科领域其最优值一般不同。

式(2)的p(X)可用Pascal三角形确定, 其线性基函数为p^T(X)=[1 x y], 式(2)简化为

$\begin{matrix} U_{s} = R_{0} a + P_{m} b, (4) \end{matrix}$

式中:U_s=[u₁u₂u₃ … u_n]^T; a=[a₁a₂a₃ … a_n]^T; b=[b₁b₂b₃ … b_n]^T。R₀与P_m的展开式分别为

$\begin{matrix} \begin{matrix} R_{0} = [\begin{matrix} R_{1} (r_{1}) & R_{2} (r_{1}) & R_{3} (r_{1}) & \dots & R_{n} (r_{1}) \\ R_{1} (r_{2}) & R_{2} (r_{2}) & R_{3} (r_{2}) & \dots & R_{n} (r_{2}) \\ R_{1} (r_{3}) & R_{2} (r_{3}) & R_{3} (r_{3}) & \dots & R_{n} (r_{3}) \\ ︙ & ︙ & ︙ & ︙ & ︙ \\ R_{1} (r_{n}) & R_{2} (r_{n}) & R_{3} (r_{n}) & \dots & R_{n} (r_{n}) \end{matrix}], (5) \\ P_{m} = [\begin{matrix} 1 & x_{1} & y_{1} & x_{1} y_{1} & \dots & p_{m} (x_{1}) \\ 1 & x_{2} & y_{2} & x_{2} y_{2} & \dots & p_{m} (x_{2}) \\ 1 & x_{3} & y_{3} & x_{3} y_{3} & \dots & p_{m} (x_{3}) \\ ︙ & ︙ & ︙ & ︙ & ︙ & ︙ \\ 1 & x_{n} & y_{n} & x_{n} y_{n} & \dots & p_{m} (x_{n}) \end{matrix}], (6) \end{matrix} \end{matrix}$

式(4)中有n+m个变量, 可使用下面m个约束条件添加m个方程

$\begin{matrix} \end{matrix}$

联立式(4)与式(7), 可得到矩阵方程

$\begin{matrix} U_{g} = [\begin{matrix} U_{s} \\ 0 \end{matrix}] = [\begin{matrix} R_{0} & P_{m} \\ {P^{T}}_{m} & 0 \end{matrix}] [\begin{matrix} a \\ b \end{matrix}] = G a_{0}, (8) \end{matrix}$

式中:U_g=[u₁u₂u₃ … u_n 0 0 … 0]^T; a₀=[a₁a₂ … a_n b₁b₂ … b_m]^T。因为矩阵R₀是对称的, 故矩阵G也是对称的。求解式(8)可得

$\begin{matrix} a_{0} = [\begin{matrix} a \\ b \end{matrix}] = G^{- 1} U_{g} 。 (9) \end{matrix}$

将式(4)写为矩阵形式

$\begin{matrix} \begin{matrix} U^{H} (X) & = R^{T} (X) a + p^{T} (X) b = [R^{T} (X) p^{T} (X)] [\begin{matrix} a \\ b \end{matrix}], (10) \end{matrix} \end{matrix}$

将式(9)代入式(10)得

$\begin{matrix} \begin{matrix} U^{H} (X) & = [R^{T} (X) p^{T} (X)] G^{- 1} U_{g} = {Φ^{T}}_{g} (X) U_{g}, \\ (11) \end{matrix} \end{matrix}$

式中 $\begin{matrix} {Φ^{T}}_{g} \end{matrix}$ (X)=[ϕ ₁(X) ϕ ₂(X) … ϕ _n(X) ϕ _n+₁(X) ϕ _n+₂(X) … ϕ _n+m(X)]。取 $\begin{matrix} {Φ^{T}}_{g} \end{matrix}$ (X)的前n项即得计算点X在支持域内的RPIM形函数Φ ^T(X):

$\begin{matrix} Φ^{T} (X) = [ϕ_{1} (X) ϕ_{2} (X) \dots ϕ_{n} (X)], (12) \end{matrix}$

相应的导数表达式为

$\begin{matrix} \begin{matrix} Φ^{(k)} (X) & = \{\frac{\partial^{k} [R^{T} (X)]}{\partial x^{k}} \frac{\partial^{k} [p^{T} (X)]}{\partial x^{k}}\} = \\ {ϕ_{1}^{(k)} (X) ϕ_{2}^{(k)} (X) \dots ϕ_{n}^{(k)} {(X)}}^{T} 。 (13) \end{matrix} \end{matrix}$

RPIM形函数满足Kronecker delta函数性质, 第一类边界条件加载便利; 然而, RBF中包含了对RPIM计算精度有直接影响的形状参数, 其最优值需要通过数值试验获得。

2.2 RPIM系统方程的构造

定义了无网格形函数, 便可以构造式(1)的离散系统方程, 将变分符号代入式(1), 有

式中:Φ 为形函数矩阵; n为支持域内的节点数; u为支持域内n个节点的场向量。将式(16)、式(15)与式(13)代入式(14), 引入总体编号体系, 式(14)最终变为:

式(17)右端的积分代表点源对二维电场的贡献, δ 函数积分存在式(18)的数值解

$\begin{matrix} \int_{Ω} δ (P) dΩ = \int_{Ω} δ (x_{0}) δ (y_{0}) dΩ = \{\begin{matrix} \frac{1}{2} & P \in Ω \\ 0 & P \notin Ω \end{matrix} (18) \end{matrix}$

式中:x₀与y₀为点源P所在位置的坐标。将式(18)代入式(17)得

$\begin{matrix} \begin{matrix} δ U^{T} (KU - P) = 0, (19) \\ K = \int_{Ω} [σ (\frac{\partial ϕ_{i}}{\partial x} \frac{\partial ϕ_{j}}{\partial x} + \frac{\partial ϕ_{i}}{\partial y} \frac{\partial ϕ_{j}}{\partial y}) + k^{2} σ ϕ_{i} ϕ_{j}] dΩ + \\ \int_{Γ_{\infty}} [\frac{kσ K_{1} (kr)}{K_{0} (kr)} \cos (r, n)] ϕ_{i} ϕ_{j} dΓ, \\ P = {[0 \dots 0 \frac{I}{2} 0 \dots 0]}^{T}; \end{matrix} \end{matrix}$

由于δ U^T是任意的, 故式(19)成立的条件为

$\begin{matrix} KU = P 。 (20) \end{matrix}$

式(20)即为无网格法构造的系统方程。由于求解域为矩形, 具体计算时应先将K包含的无穷边界积分项展开为

$\begin{matrix} \begin{matrix} \int_{Γ_{\infty}} [\frac{kσ K_{1} (kr)}{K_{0} (kr)} \cos (r, n)] ϕ_{i} ϕ_{j} dΓ = \\ \int_{Γ_{AB}} [\frac{kσ K_{1} (kr)}{K_{0} (kr)} (- \frac{| x_{0} - x |}{r})] ϕ_{i} ϕ_{j} dΓ + \\ \int_{Γ_{BC}} [\frac{kσ K_{1} (kr)}{K_{0} (kr)} (- \frac{| y |}{r})] ϕ_{i} ϕ_{j} dΓ + \\ \int_{Γ_{CD}} [\frac{kσ K_{1} (kr)}{K_{0} (kr)} (- \frac{| x_{0} - x |}{r})] ϕ_{i} ϕ_{j} dΓ, (21) \end{matrix} \end{matrix}$

式中r= $\begin{matrix} \sqrt{(x_{0} {- x)}^{2} + (y_{0} {- y)}^{2}} \end{matrix}$ 。求解出波数域的电位U后, 就可利用如表1所示的一组最优化波数^[16], 将波数域电位转化为空间域电位转换表达式为

表1 最优化波数

2.3 求解域背景网格积分

K的表达式中包含对求解域Ω 与求解域边界Γ 的积分, 为计算这些积分, 一般采用一组背景网格将区域离散(图1), 总体积分表示成离散单元积分之和的形式, 每个单元的积分利用高斯积分法求解:

式中:n_c与n_cΓ分别为背景单元及边界单元的个数; n_g与n_gΓ分别为背景单元及边界单元内的高斯点数目, ω _i为第i个高斯点X_Qi的权系数; $\begin{matrix} J_{ik}^{D} \end{matrix}$ 与 $\begin{matrix} J_{il}^{B} \end{matrix}$ 分别为背景单元k在积分点X_Qi处面积分及子边界l在积分点X_Qi处曲线积分的雅可比(Jacobian)矩阵。式(23)即为总体刚度矩阵的最终离散表达式。

3 RPIM最优形状参数试验

首先通过不含形状参数的EFGM来试验点源二维正演问题支持域无量纲尺寸的最优值, 计算对象为ρ =1 000 Ω · m的均匀介质模型。表2为不同支持域无量纲尺寸的EFGM数值计算结果, 由表2可知, 点源二维正演问题α 最优值为1.0~1.2。由于α 增大的同时增加了计算耗时^[13], 故之后的RPIM最优形状参数试验取α =1.0。

表2 不同支持域无量纲尺寸EFGM数值计算结果

表3~表6列出了基于MQ函数、EXP函数及TPS函数的RPIM在不同形状参数下的计算结果。如表3~表4所示, MQ函数形状参数具有较大的动态范围, α _c取值区间约[-6, 6], q取值区间约[-44, 9], 超过此区间, RPIM计算会产生奇异; 最优形状参数α _c=1.0, q=-5.8。表5、表6显示EXP函数最优形状参数α _c=6.3, TPS函数最优形状参数η =1.43或η =1.6。

表3 α _c=1.0时不同q值的RPIM(MQ)计算结果

表4 q=-5.8时不同α _c值的RPIM(MQ)计算结果

表5 RPIM(EXP)不同形状参数计算结果

表6 RPIM(TPS)不同形状参数计算结果

RPIM(MQ)处理点源二维正演问题具有很好的普适性, 其形状参数动态范围最大, RPIM(EXP)次之; 计算精度由高到低依次为RPIM(MQ)→ RPIM(TPS)→ RPIM(EXP)。

表7列出了最优形状参数下不同高斯点数量的RPIM计算结果。由表可见, 高斯点数量的增加对计算精度无明显提升效果, 2~4点高斯积分公式已能满足精度要求。

表7 不同高斯点数量RPIM计算结果

4 结论

1) 利用径向基函数近似求解域中的点源二维电位具有边界条件加载便利、精度高的优点。RBF形状参数对RPIM计算精度影响大, 通过数值试验获得了适用于直流点源问题MQ(multi-quadrics)函数、薄板样条(thin plate spline, TPS)函数、高斯指数(EXP)函数形状参数及支持域无量纲尺寸α 的建议值或区间:RPIM(MQ)形状参数建议取α _c=1.0, q=-5.8; RPIM(EXP)建议取α _c=6.3; RPIM(TPS)建议取η =1.43或η =1.6, 建议α 为1.0~1.2。

2) 对比了最优形状参数推荐值下RPIM模拟点源电阻率参量的精度, 相同高斯积分点RBF近似精度由高至低为MQ、TPS、EXP; 增加背景单元内的高斯点对提高RPIM求解精度作用有限, 反而会降低计算效率, 选用2~4点高斯积分公式即可。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	Key K, Weiss C. Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example[J]. Geophysics, 2006, 71(6): 291-299. [本文引用:1]
[2]	汤井田, 王飞燕. 基于非结构化网格的2. 5D直流电阻率模拟[J]. 物探化探计算技术, 2008, 30(5): 413-418. [本文引用:1]
[3]	刘长生, 汤井田, 任政勇, 等. 基于非结构化网格的三维大地电磁自适应矢量有限元模拟[J]. 中南大学学报: 自然科学版, 2010, 41(5): 1855-1859. [本文引用:1]
[4]	Ren Z Y, Tang J T. 3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method[J]. Geophysics, 2010, 75(1): 7-17. [本文引用:1]
[5]	张力员, 阮百尧, 刘海飞, 等. 三维全空间坑道直流聚焦超前探测数值模拟[J]. 物探与化探, 2011, 35(3): 419-422. [本文引用:1]
[6]	赵晓博, 朱自强, 李建慧. 基于非结构化网格的瞬变电磁2. 5维有限元正演模拟[J]. 物探化探计算技术, 2011, 33(5): 517-521. [本文引用:1]
[7]	严波, 刘颖, 叶益信. 基于对偶加权后验误差估计的2. 5维直流电阻率自适应有限元正演[J]. 物探与化探, 2014, 38(1): 145-150. [本文引用:1]
[8]	赵宁, 秦策. 基于非结构化网格的2. 5维时间域激发极化反演[J]. 科学技术与工程, 2014, 14(5): 17-23. [本文引用:1]
[9]	李俊杰, 严家斌. 无网格法进展及其在地球物理学中的应用[J]. 地球物理学进展, 2014, 29(1): 452-461. [本文引用:1]
[10]	贾晓峰, 胡天跃, 王润秋. 复杂介质中地震波模拟的无单元法[J]. 石油地球物理勘探, 2006, 41(1): 43-48. [本文引用:1]
[11]	冯德山, 王洪华, 戴前伟. 基于无单元Galerkin法探地雷达正演模拟[J]. 地球物理学报, 2013, 56(1): 298-308. [本文引用:1]
[12]	苏洲, 胡文宝. 二维大地电磁正演中的无网格算法[J]. 物探与化探, 2012, 36(6): 1024-1029. [本文引用:1]
[13]	严家斌, 李俊杰. 无网格法在大地电磁正演计算中的应用[J]. 中南大学学报: 自然科学版, 2014, 45(10): 3513-3520. [本文引用:2]
[14]	李俊杰, 严家斌. 无网格点插值法大地电磁二维正演数值模拟[J]. 石油物探, 2014, 53(5): 617-626. [本文引用:1]
[15]	刘云, 宋滔, 王赟. 电导率连续变化2. 5维直流电阻率法有限元数值模拟[J]. 地球物理学进展, 2014, 29(3): 1194-1200. [本文引用:1]
[16]	徐世浙. 地球物理中的有限单元法[M]. 北京: 科学出版社, 1994: 159-165. [本文引用:1]

2006

0.0

2008

0.0

2010

0.0

2010

0.0

2011

0.0

2011

0.0

2014

0.0

2014

0.0

2014

0.0

... 无网格法^[9]作为FEM的补充和发展,是近十几年来兴起于计算力学领域的一类新型数值算法,其物性加载在只与坐标位置有关的高斯点上,故节点规则分布下也能适应复杂模型的计算 ...

2006

0.0

... 无单元Galerkin法(element-free galerkin method,EFGM)与点插值法(point interpolation method,PIM)作为两种较成熟的全域弱式无网格方法,已被成功应用于地震波场^[10]、雷达波场^[11]与大地电磁场响应 1214的计算中 ...

2013

0.0

2012

0.0

2014

0.0

... 用于控制实际支持域的大小,是对无网格法计算精度影响很大的重要参数^[13] ...

... 增大的同时增加了计算耗时^[13],故之后的RPIM最优形状参数试验取#cod#x003b1 ...

2014

0.0

2014

0.0

... 1 点源二维电场波数域变分问题点源为三维源,需先采用傅里叶变换将三维点源从空间域转化为波数域二维源,当地下电性结构为二维时,取走向为z轴,x轴与z轴垂直,y轴垂直向上,求解域为矩形区域,4个顶点依次以A、B、C、D逆时针编号,波数域电位满足变分问题^[15] ...

0.0

... 求解出波数域的电位U后,就可利用如表1所示的一组最优化波数^[16],将波数域电位转化为空间域电位转换表达式为 ...