多连通域上泊松方程柯西问题的一种新型无网格方法

Cauchy Problem of the Poisson Equation in a Multi-Connected Domain by a New Meshless Method

摘要: 本文提出了一种求解多连通域中泊松方程柯西问题的无网格数值方法。结合拉普拉斯方程的基本解和径向基函数得到了数值解。由于系数矩阵是不适定的，因此采用正则化方法来求解所得到的线性方程组。通过对正则化参数的适当选取和对柯西数据的先验假设，得到了上述问题的正则化解，并且利用数值例子验证了该方法的有效性和准确性。

Abstract: This paper presents a meshless numerical scheme to solve the Cauchy problem of the Poisson equation in a multi-connected domain. Fundamental solutions of Laplace’s equations and radial basis functions (RBFs) are used to obtain a numerical solution. Because the coefficient matrix is illposed, the Tikhonov regularization method is applied to solve the resulting system of linear equations. By the suitable choices of a regularization parameter and a priori assumption to the Cauchy data, the regularized solution to above problem is obtained. Several numerical examples are given to verify the efficiency and accuracy of the proposed method.

1. 引言

本文讨论了多连通域中泊松方程的一个柯西问题。设 $\Omega$ 是 $R^2$ 中一个有界的多连通区域，有两个光滑边界 ${\Gamma }_{\text{1}}$ 和 ${\Gamma }_{\text{0}}$。假设 ${\Gamma }_{\text{1}}$ 是可以测量柯西数据的外部可达边界， ${\Gamma }_{\text{0}}$ 为内部不可达边界。函数 $u\in C^{\text{2}}\left(\Omega \right)\cap C^{\text{1}}\left(\bar{\Omega }\right)$ 满足，

$-\Delta u=f\left(x,y\right),\text{in}\Omega$ (1)

以及边界条件

${u|}_{{\Gamma }_1}=g_1\left(x,y\right)$ (2)

${\frac{\partial u}{\partial n}|}_{{\Gamma }_1}=g_2\left(x,y\right)$ (3)

其中 $g_1$ 和 $g_2$ 分别是在边界 ${\Gamma }_{\text{1}}$ 上指定的Dirichlet和Neumann数据； $\frac{\partial u}{\partial n}$ 表示 $u$ 在 ${\Gamma }_{\text{1}}$ 上的外法向导数。

众所周知，在Hadamard [1] 意义下，泊松方程的柯西问题是不适定的 [2] [3]。也就是说，给定数据中非常小的扰动可能会导致解的巨大偏差。该问题在等离子体物理学、心电图学、腐蚀无损评价等许多不同的领域都有应用 [4] [5] [6] [7]。由于其缺陷，柯西问题的数值逼近变得非常困难和具有挑战性。通常，我们使用正则化技术来获得稳定的数值解，如Tikhonov正则化 [8] 和拟逆方法 [9]。

对于环域中具有Dirichlet边界条件的泊松方程的直接问题，Lin等人研究了一种基于傅里叶展开法和超球面光谱方法的简单而有效的光谱方法 [10]。魏等人用基本解法结合Tikhonov正则化法求解了圆的邻域上拉普拉斯方程的柯西问题 [11]。然而，据我们所知，关于泊松方程柯西问题的文献很少，本文的目的是用MFS-RBFs来解决这个问题。

最近，作者提出了一种新的，改进的无网格方法来解决逆热源问题 [12]。他们结合了MFS和RBF方法给出了逆源问题的数值解。受他们思想的启发，我们用这种方法来求解泊松方程。

通常对于实际问题，由于 $g_1$ 和 $g_2$ 是通过测量得到的，因此不可避免地会有误差 $\delta$，而且我们只知道测量的数据 $g_1^{\delta }$， $g_2^{\delta }$。假设

$\left|g_1^{\delta }\left(x,y\right)-g_{\text{1}}\left(x,y\right)\right|\le \delta \left|g_2^{\delta }\left(x,y\right)-g_2\left(x,y\right)\right|\le \delta 对所有\left(x,y\right)\in {\Gamma }_1$ (4)

本文的组织结构如下：在第1节和第2节中，我们分别介绍了拉普拉斯方程的基本解方法和径向基函数法。第3节中应用基本解–径向基函数(MFS-RBFs)方法解决问题。第4节介绍了Tikhonov正则化方法。最后，我们在第5节中给出了数值例子。

2. 基本解方法

基本解法是求解偏微分方程的常用方法，是一种无网格方法。其主要思想是为一个给定的偏微分方程找到一种具有特殊形式和奇异性的解，然后假设这些解可以得到所需的定解问题的解。然而，基本解方法有一个很大的缺点，即使是对于正定问题，所得到的线性方程系统也是病态的 [13]。因此，需要一种特殊的正则化方法来求解这个代数系统。

在本节中，我们给出了一种用周和魏在 [11] [14] 中提出的MFS近似拉普拉斯方程柯西问题解的方法。在(1)中，当 $f\left(x,y\right)\equiv 0$ 时，原方程变成拉普拉斯方程

$\Delta u=0\text{in}\Omega$ (5)

设 $R_{\text{1}}$ 和 $R_{\text{2}}$ 为实数，使 $0<R_2<R_1$，而N为自然数。我们在以下点上安排内外源点：

$Q_{\text{1},j}=\left(R_1\cos \left({\theta }_j\right),R_1\sin \left({\theta }_j\right)\right),j=1,2,\cdots ,N,$ 为外源点，

$Q_{\text{2},j}=\left(R_{\text{2}}\cos \left({\theta }_j\right),R_{\text{2}}\sin \left({\theta }_j\right)\right),j=1,2,\cdots ,N,$ 为内源点，其中 ${\theta }_j=\frac{2\text{π}j}{N}$。

我们知道拉普拉斯方程的基本解 [15] 定义为

$G\left(P,Q\right)=-\frac{1}{2\pi }\log \left|P-Q\right|$， (6)

其中 $P$ 和 $Q$ 是 $R^2$ 中的点， $\left|P-Q\right|$ 表示 $R^2$ 中的欧几里得距离。

根据MFS的理念和基本解(6)，我们假设柯西问题(5)，(2)，(3)的近似解 $u_N$ 可以用以下线性组合来表示

$u_N\left(x\right)=q+{\displaystyle \underset{s=1}{\overset{2}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\lambda }_{s,j}}}G\left(P,Q_{s,j}\right)$， (7)

其中 $q$ 是一个实常数， $\left\{Q_{s,j}\right\}$ 是源点，而 $\left\{{\lambda }_{s,j}\right\}$ 是要确定的未知系数，从而满足以下不变条件

${\displaystyle \underset{s=1}{\overset{2}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\lambda }_{s,j}}}=0$ (8)

在可接触边界 ${\Gamma }_{\text{1}}$ 上，选择N个配置点 $\left\{P_i\right\}=\left(r\cos {\theta }_i,r\sin {\theta }_i\right)$， $i=1,2,\cdots ,N$ 来拟合Dirichlet数据和Neumann数据。有关对配置点和源点的设置的说明，请参见图1。通过配置边界条件(2)~(3)，容易得到线性方程组，

$\{\begin{array}{l}{u}_N\left(P_i\right)=g_1\left(P_i\right),i=1,2,\cdots ,N,\\ \frac{\partial u_N}{\partial n}\left(P_i\right)=g_2\left(P_i\right),i=1,2,\cdots ,N.\end{array}$ (9)

Figure 1. Fork(×) are collocation points for boundary condition; Stars(*) represent source points; Solid line (-) denotes the outer boundary; Dashdot line (-.) represents the inner boundary

图1. (×)是边界条件的配置点；(*)表示源点；实线(-)表示外边界Γ₁；点线(-.)表示内边界Γ₀

从不变条件(8)和配置条件(9)可得，参数集 $q$ 和 ${\lambda }_{s,j}$ 是以下线性方程的解

$A\lambda =b$ (10)

其中 $\lambda ={\left(q{\lambda }_1{\lambda }_2\right)}^{\text{T}},b={\left(0g_1{g}_2\right)}^{\text{T}}$，以及 $A$ 由以下矩阵定义

$\left(\begin{array}{ccc}\text{0}& 1^{\text{T}}& 1^{\text{T}}\\ 1& L_1& L_2\\ 0& M_1& M_2\end{array}\right)$， (11)

其中

$L_s=\left(G\left(P_i,Q_{s,j}\right)\right)\in R_{N\times N},s=1,2,$

$M_s=\left(\frac{\partial G}{\partial n}\left(P_i,Q_{s,j}\right)\right)\in R_{N\times N},s=1,2,$

$1={\left(1,1,\cdots 1\right)}^{\text{T}}\in R^N,$

$0={\left(0,0,\cdots ,0\right)}^{\text{T}}\in R^N,$

${\lambda }_s={\left({\lambda }_{s,0},{\lambda }_{s,1},\cdots ,{\lambda }_{s,N}\right)}^{\text{T}}\in R^N,s=1,2,$

$g_1={\left(g_1\left(P_1\right),g_1\left(P_2\right),\cdots ,g_1\left(P_N\right)\right)}^{\text{T}}\in R^N,$

$g_2={\left(g_2\left(P_1\right),g_2\left(P_2\right),\cdots ,g_2\left(P_N\right)\right)}^{\text{T}}\in R^N.$

3. 径向基函数方法

在本节中，我们考虑了离散数据 [12] 的插值的RBF方法。如果 $X^{\text{*}}$ 是一个不动点， $X$ 是 $R^{\text{2}}$ 中的一个任意点，而径向函数 ${\psi }^{*}$ 通过 ${\psi }^{\text{*}}={\psi }^{*}\left(r\right)$ 定义，其中 $r={\Vert X-X^{*}\Vert }_2$，径向函数 ${\psi }^{*}$ 仅由 $X$ 和 $X^{\text{*}}$ 之间的距离决定。这个性质意味着 ${\psi }^{*}$ 是径向对称的。表1中给出了一些常用的无限光滑的RBFs。这些函数依赖于一个自由参数C，称为形状参数，它使用了RBFs的近似理论。

Table 1. Some well-known radial basis functions

表1. 一些常用的径向基函数

设 ${\left\{\left(x_k,y_k\right)\right\}}_{k=1}^N$ 是 $R^{\text{2}}$ 中的域 $\Omega$ 中一组给定的不同点的集合。RBFs的主要思想是使用相同类型的RBFs的线性组合进行如下插值：

$\psi \left(x,y\right)={\displaystyle \underset{k=1}{\overset{N}{\sum }}{\lambda }_k{\psi }_k\left(x,y\right),}$ (12)

其中 ${\psi }_k\left(x,y\right)=\psi \left(\Vert \left(x,y\right)-\left(x_k,y_k\right)\Vert \right)$， ${\lambda }_k$ 是未知系数。假设给定值 $f_k=f\left(x_k,y_k\right),k=1,2,\cdots ,N$。当未知系数 ${\lambda }_k$ 的值已知时，有 $\psi \left(x_k,y_k\right)=f_k$，且满足以下线性方程

$A\lambda =b$， (13)

这里 $\lambda ={\left({\lambda }_1,\cdots ,{\lambda }_N\right)}^T$， $b={\left(f_1,\cdots ,f_N\right)}^{\text{T}}$ 和 $A=\left(a_{ij}\right)$，其中 $i,j=1,2,\cdots ,N,$ 对应的 $a_{ij}={\psi }_j\left(x{}_i,y_i\right)$。

由Schoenberg定理 [16] 可知，在一般情况下，对于高斯，逆多二次，多二次以及逆二次径向基函数中不同的插值点，矩阵是正定的。

4. 基本解–径向基函数方法

本节采用基本解和径向基函数法求解问题(1)~(3)的数值方案。通过MFS和RBFs的线性组合，直接给出解的线性表达式，即用RBFs的线性组合近似出方程(1)的非齐次项。

在区域 $\Omega$ 中，我们假设

$Q_{s,j}=\left\{\left(x_{s,j},y_{s,j}\right)\right\},s=1,2,j=1,2,\cdots ,N,$ $P_i=\left\{\left(x_i,y_i\right)\right\}\in {\Gamma }_1,i=1,2,\cdots ,N,$ 分别为源点和配置点。

则问题(1)~(3)的解可以被写成如下形式：

$u^{*}\left(x,y\right)=q+{\displaystyle \underset{s=1}{\overset{2}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\lambda }_{s,j}{\varphi }_{s,j}\left(x,y\right)}}+{\displaystyle \underset{s=1}{\overset{2}{\sum }}{\displaystyle \underset{j=N+1}{\overset{2N}{\sum }}{\lambda }_{s,j}{\psi }_{s,j}\left(x,y\right),}}$ (14)

这里 ${\varphi }_{s,j}\left(x,y\right)={\varphi }_{s,j}\left(x-x_{s,j},y-y_{s,j}\right)=G\left(P,Q_{s,j}\right)$，其中 $G\left(P,Q\right)$ 由(6)给出。而且 ${\psi }_{s,j}\left(x,y\right)=\psi \left({\Vert \left(x,y\right)-\left(x_{s,j},y_{s,j}\right)\Vert }_2\right)$ 是由表1给出的径向基函数中的一种。在本节中，我们选用高斯径向基函数。我们用近似解 $u^{*}\left(x,y\right)$ 去满足给定的方程。所以系数 ${\lambda }_{s,j}$ 通过以下矩阵方程确定：

$A\lambda =b,$ (15)

其中 $\lambda ={\left(q,{\lambda }_1,\cdots ,{\lambda }_{4N}\right)}^{\text{T}}$， $b={\left(f,g_1,g_2\right)}^{\text{T}}$。矩阵 $A$ 可以表示为如下形式：

$\left(\begin{array}{ccccc}\text{0}& \text{0}& \text{0}& \Delta {\psi }_{1,j}\left(x_i,y_i\right)& \Delta {\psi }_{2,k}\left(x_i,y_i\right)\\ 1& {\varphi }_{1,j}\left(x_i,y_i\right)& {\phi }_{2,j}\left(x_i,y_i\right)& {\psi }_{1,k}\left(x_i,y_i\right)& {\psi }_{2,k}\left(x_i,y_i\right)\\ 0& \frac{\partial {\varphi }_{1,j}}{\partial n}\left(x_i,y_i\right)& \frac{\partial {\varphi }_{2,j}}{\partial n}\left(x_i,y_i\right)& \frac{\partial {\psi }_{1,k}}{\partial n}\left(x_i,y_i\right)& \frac{\partial {\psi }_{2,k}}{\partial n}\left(x_i,y_i\right)\end{array}\right)$

其中 $i,j=1,2,\cdots ,N$， $k=N+1,\cdots ,2N$。

在实际应用中，由于测量数据 $g_1$， $g_2$ 和 $f$ 不能准确地测量并且通常会被噪音干扰，并且数值算法关于噪音数据的稳定性对于获得具有物理意义的解至关重要。故我们用噪音数据来证明所提出的数值方案的稳定性，并得到了以下噪声数据:

$b_i^{\delta }=b_i\left(1+\delta \cdot \left(2\ast rand\left(i\right)-1\right)\right),i=1,\cdots ,N$， (16)

其中 $b_i$ 是精确值，函数 $rand\left(i\right)$ 是在[0, 1]上均匀分布的随机数。

5. 正则化方法

由于泊松方程(1)~(3)的柯西问题是高度不适定的，方程(15)中的矩阵A是病态的，且矩阵A的条件数随着配置点数的增加而急剧增加。另一方面，由于径向基函数依赖于形状参数C也会影响A的条件数。因此，我们常用一些正则化方法来解决这种不适定问题 [17] [18]。本文中采用Tikhonov正则化方法 [8] 求解带有右侧误差 $b^{\delta }$ 的矩阵方程(13)，即解决以下最小二乘问题

$\underset{\lambda }{\min }{\Vert A\lambda -b^{\delta }\Vert }_2^2+{\alpha }^2{\Vert \lambda \Vert }_2^2,$ (17)

其中 ${\Vert \cdot \Vert }_{\text{2}}$ 是 $R^{\text{2}}$ 中的欧几里得范数， $\alpha >0$ 为正则化参数。容易知道(17)中的 ${\lambda }_{\alpha }^{\delta }$ 的最小值满足

$\left(\alpha I+A^{\text{T}}A\right){\lambda }_{\alpha }^{\delta }=A^{\text{T}}{b}^{\delta },$ (18)

其中 ${\lambda }_{\alpha }^{\delta }={\left(q_{\alpha }^{\delta },{\left({\lambda }_{\alpha }^{\delta }\right)}_{1,1},\cdots ,{\left({\lambda }_{\alpha }^{\delta }\right)}_{2,N},{\left({\lambda }_{\alpha }^{\delta }\right)}_{1,N+1},\cdots ,{\left({\lambda }_{\alpha }^{\delta }\right)}_{\text{2},\text{2}N}\right)}^{\text{T}}$。在计算过程中，我们使用Hansen的软件包 [19] 来解决(15)，而该问题(1)~(3)的近似解 $u_{\alpha }^{\delta }$ 表示如下：

$u_{\alpha }^{\delta }\left(x,y\right)=q_{\alpha }^{\delta }+{\displaystyle \underset{s=1}{\overset{2}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left({\lambda }_{s,j}\right)}_{\alpha }^{\delta }{\varphi }_{s,j}\left(x,y\right)}}+{\displaystyle \underset{s=1}{\overset{2}{\sum }}{\displaystyle \underset{j=N+1}{\overset{2N}{\sum }}{\left({\lambda }_{s,j}\right)}_{\alpha }^{\delta }{\psi }_{s,j}\left(x,y\right),}}$ (19)

因此，对 $\left(x_0,y_0\right)\in {\Gamma }_0$，内部不可达边界上的Dirichlet数据为

${u_{\alpha }^{\delta }|}_{{\Gamma }_0}=q_{\alpha }^{\delta }+{\displaystyle \underset{s=1}{\overset{2}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left({\lambda }_{s,j}\right)}_{\alpha }^{\delta }{\varphi }_{s,j}\left(x_0,y_0\right)}}+{\displaystyle \underset{s=1}{\overset{2}{\sum }}{\displaystyle \underset{j=N+1}{\overset{2N}{\sum }}{\left({\lambda }_{s,j}\right)}_{\alpha }^{\delta }{\psi }_{s,j}\left(x_0,y_0\right),}}$ (20)

6. 数值实验

在本节中，我们举了三个数值例子来验证该方法的可行性和有效性。计算程序都在Matlab中给出，为了估计数值近似的误差，我们选择了一些测试点，用以下公式计算相对均方根误差(RRMSE)。

$RRMSE\left(u\right)=\frac{\sqrt{{\displaystyle {\sum }_{i=1}^N{\left({\tilde{u}}_i-u_i\right)}^2}}}{\sqrt{{\displaystyle {\sum }_{i=1}^N{u}_i^2}}}$

其中N是在内边界 ${\Gamma }_{\text{0}}$ 上均匀分布的测试点的总数， ${\tilde{u}}_i$ 和 $u_i$ 是在 ${\Gamma }_{\text{0}}$ 上一点处的解的近似值和精确值。

例1. 已知精确解以及内外边界都为圆。

$u\left(x,y\right)=-\frac{y}{{\left(x-0.1\right)}^2+y^2}+\cos (\; x\; )$

${\Gamma }_1=\left\{\left(x,y\right)|x^2+y^2=1\right\}$

${\Gamma }_0=\left\{\left(x,y\right)|x^2+y^2=0.8\right\}$

Figure 2. The approximate solutions ( ${u_{\alpha }^{\delta }|}_{{\Gamma }_0}$ ) to different points on ( ${\Gamma }_0$ ) in Example 1

图2. 例一中( ${\Gamma }_0$ )上不同点的近似解 ${u_{\alpha }^{\delta }|}_{{\Gamma }_0}$

例2. 精确解和内外边界分别为

$u\left(x,y\right)=\left(1-\exp \left(-y\right)\right)\sin (\; x\; )$

${\Gamma }_1=\left\{\left(x,y\right)|x^2+y^2=\text{3}\right\}$

${\Gamma }_{\text{0}}=\left\{\left(x,y\right)|x^2+y^2=\text{2}\text{.6}\right\}$

Figure 3. The approximate solutions ( ${u_{\alpha }^{\delta }|}_{{\Gamma }_0}$ ) to different points on ( ${\Gamma }_0$ ) in Example 2

图3. 例二中( ${\Gamma }_0$ )上不同点的近似解 ${u_{\alpha }^{\delta }|}_{{\Gamma }_0}$

例3. 在此例中，我们考虑了精确解u未知的情况。设外边界是一个圆 ${\Gamma }_1=\left\{\left(x,y\right)|x^2+y^2=\text{2}\right\}$，内边界通过 $\left({\Gamma }_0\ast \cos \left(t\right),{\Gamma }_0\ast \sin \left(t\right)\right)$ 被参数化， ${\Gamma }_{\text{0}}$ 为

${\Gamma }_{\text{0}}=\frac{\text{5}}{\text{6}}\ast \sqrt{\cos {\left(t\right)}^2+0.25\ast \sin {\left(t\right)}^{\text{2}}},$

$t\in \left[0,2\text{π}\right].$

在 [14] 中给出的。通过求解以下直接问题，得到了外边界 ${\Gamma }_{\text{1}}$ 上的Neumann数据。

$\Delta u=0,\text{in}D$

${u|}_{{\Gamma }_1}=2+x,$

${u|}_{{\Gamma }_0}=\exp \left(x\right).$

在这个例子中，我们可以看到非齐次项是0，我们使用本文中提到的方法来解决它，即通过径向基函数来逼近0。

Figure 4. The approximate solutions ( ${u_{\alpha }^{\delta }|}_{{\Gamma }_0}$ ) to different points on ( ${\Gamma }_0$ ) in Example 3

图4. 例三中( ${\Gamma }_0$ )上不同点的近似解 ${u_{\alpha }^{\delta }|}_{{\Gamma }_0}$

图2~4展示了例1~3中 ${u_{\alpha }^{\delta }|}_{{\Gamma }_0}$ 的精确解和近似解之间的比较，以及使用MFS-RBFs添加到数据中的不同误差。结果表明，随着误差水平的增加，该近似函数具有可接受的精度。在图5中我们展示了例1~3的关于 $\delta =1\%$ 的形状参数C的相对误差RRMSE(u)。在图6中，我们分别绘制了在内边界 ${\Gamma }_{\text{0}}$ 上例1和2在不同误差水平下的解 $u$ 的Neumann数据。图7给出了例1~3的数值解相对于配置点的相对误差。可以看出，数值误差逐渐稳定，表明不需要使用大量的配置点来得到精确的解。然而，从图8中来看，当源点的数量在一定的范围内时，数值结果是振荡的。

(a) Relative error change of shape parameter C in Example 1 (b) Relative error change of shape parameter C in Example 2(c) Relative error change of shape parameter C in Example 3

Figure 5. Relative error change of shape parameter C in Examples 1~3

图5. 例1~3中δ = 1%时形状参数C的相对误差变化

(a) $R_1=2,R_2=0.5,C=0.001,N=65$ (b) $R_1=\text{4},R_2=\text{2},C=0.1,N=\text{350}$

Figure 6. The ${\frac{\partial u}{\partial n}|}_{{\Gamma }_0}$ for Examples 1 and 2 with various noise levels

图6. 例1和例2不同误差水平下的 ${\frac{\partial u}{\partial n}|}_{{\Gamma }_0}$

(a) $R_1=2,R_2=0.5,C=0.1,N=65,\delta =0.01$ (b) $R_1=\text{4},R_2=\text{1},C=0.1,N=\text{170},\delta =0.01$ (c) $R_1=\text{10},R_2=\frac{\text{1}}{\text{6}},C=0.\text{0002},N=\text{100},\delta =0.01$

Figure 7. The relative errors of solutions with respect to the number of collocation points for Example 1~3

图7. 例1~3的关于解的配置点的数目的相对误差

(a) $R_1=2,R_2=0.5,C=0.1,\alpha =\text{0}\text{.001},\delta =0.01$ (b) $R_1=\text{4},R_2=\text{1},C=0.1,\alpha =\text{0}\text{.01},\delta =0.01$

Figure 8. The relative errors of solutions with respect to the number of source points for Example 1 and 2

图8. 例1和例2的解相对于源点数目的相对误差

7. 结论

本文将MFS与RBF结合起来，重建了多连通域上二维泊松方程的内部Dirichlet和Neumann数据。用具有精确解和非精确解的数值例子验证了该方法是准确和稳定的。

基金项目

本文所述的工作由国家自然科学基金项目数学天元基金(No. 11326234)，甘肃省自然科学基金项目(No. 145RJZA099)，甘肃省高校科研项目(2014A-012)，西北师范大学青年教师科研能力提升项目(2020-08)支持。

NOTES

^*通讯作者。