文章编号：1004-0609(2013)10-2897-08

基于混沌振荡PSO-BP算法的电阻率层析成像非线性反演

戴前伟^{1, 2}，江沸菠^{1, 2}

(1. 中南大学 有色金属成矿预测教育部重点实验室，长沙 410083；

2. 中南大学 地球科学与信息物理学院，长沙 410083)

摘  要：粒子群优化算法是一种启发式的全局优化算法，将其与BP神经网络结合，能够有效地改善BP神经网络在进行电阻率层析反演中的收敛速度和求解质量。提出一种基于混沌振荡的粒子群算法，使用混沌振荡曲线来自适应调整惯性权重w以提高PSO算法的全局寻优能力，并使用其训练和优化BP神经网络的权值和阈值。比较不同隐含层节点数目和惯性权重w值对反演结果的影响，并给出混沌振荡PSO-BP算法非线性反演的具体实现方案。对均匀半空间中异常体理论模型进行反演，实验结果表明：混沌振荡PSO-BP不依赖初始模型，在稳定性和准确性上优于BP反演和标准PSO-BP反演，成像质量优于最小二乘法反演的。

关键词：电阻率层析成像；非线性反演；粒子群优化；反向传播网络；混沌序列

中图分类号：P631        　   　     文献标志码：A

Nonlinear inversion for electrical resistivity tomography based on chaotic oscillation PSO-BP algorithm

DAI Qian-wei^{1, 2}, JIANG Fei-bo^{1, 2}

(1. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education, Central South University, Changsha 410083, China;

2. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)

Abstract: The particle swarm optimization (PSO) is a heuristic global optimization method, which can effectively improve the convergence speed and the results quality with the BP neural network in resistivity tomography 2-D nonlinear inversion. A chaotic oscillation PSO algorithm was presented, and the chaos oscillation curve was used to adjust the inertia weight adaptively and improve the global optimum capability of PSO. And this algorithm was used to train and optimize the weights and threshold values of the BP neural network. The impacts of different numbers of the hidden layer nodes and types of the inertia weight to the inversion result were compared, and an implementation of chaotic oscillation PSO-BP algorithm was given. The half space abnormity synthetic model was inversed. The results show that the chaotic oscillation PSO-BP algorithm that is independent of the initial model has better performance than BP and standard PSO-BP algorithm in stability and accuracy, and has higher imaging quality than least square inversion.

Key words: electrical resistivity tomography; nonlinear inversion; particle swarm optimization; back propagation network; chaotic sequences

1987年，SHIMA等^[1]首次提出了“电阻率层析成像”(Resistivity tomography)一词，并提出了反演解释的方法。此后各国学者从不同角度对这一问题从理论、实验到应用开展了研究，电阻率层析成像ERT技术得到了快速发展，由于该方法场源易于实现，野外实施简便，对现场无破坏作用，分辨率高(与普通电阻率法比较)，勘探成本低(与地震波层析成像比较)，因此，穿透深度大(与电磁波层析成像比较)，特别适合于浅层构造和工程探测等研究领域^[2]。

目前，ERT资料处理、解释以线性反演为主，例如α中心法^[3]、最小二乘法^[4]、佐迪法^[5]、积分法^[6]等。由于ERT反演本质上是一个复杂的非线性问题，将非线性问题在初始模型附近进行线性化，存在着容易陷入局部收敛、依赖初始模型选择、偏导数矩阵求解困难等问题^[7]，因此，ERT的非线性反演方法的研究越来越受到人们的关注。

在ERT的非线性反演研究方面：卢元林等^[2]将模拟退火算法较早地引入二维电法资料的反演，MADAN等^[8]、等^[9]使用遗传算法实现了电测深曲线的遗传算法反演，LIU等^[10]使用改进的遗传算法进行了三维的电阻率反演，等^[11]将粒子群优化算法引入一维电阻率反演的过程中，取得了较好的效果。与以上这些非线性反演方法相比，神经网络由于能够模拟人脑神经元的信息学习和处理，直接从样本数据中构造数学模型，所以在电法的资料的反演解释中应用得最为广泛。GAD等^[12]使用人工神经网络进行了一维和二维电阻率反演的初步研究，但是其二维反演的网格剖分稀疏，反演的参数不多、规模有限；MAITI等^[13]提出一种贝叶斯神经网络的训练方法，但是其应用仅限于一维电测深数据；AHMAD等^[14]和徐海浪等^[7]使用人工神经网络对不同装置的高密度电阻率数据进行了训练并实现了反演，得到了较好的结果，但是两者的样本划分方式差别很大，不利于神经网络算法在反演中的应用推广。同时，两者均未能有效解决BP神经网络在反演时存在的局部极值问题。综上所述，虽然神经网络在电法资料的反演解释中得到了广泛的应用，但是其搜索空间随着模型参数的增加而急剧增加，在求解复杂反演问题时，其算法的收敛性和准确性成为了研究的重点。

为了解决单一非线性反演方法在解释复杂ERT资料时存在的局部极值和收敛缓慢等问题，本文作者提出一种基于混沌振荡的PSO算法，将混沌振荡PSO算法和BP神经网络相结合，使用PSO算法来优化神经网络的权值和阈值，改善BP神经网络反演的性能，针对不同的隐含层节点数目和惯性权重w值对混沌振荡PSO-BP算法的性能进行比较和分析，得到优化的混沌振荡PSO-BP算法模型，并将其性能与BP算法和标准PSO-BP算法进行对比，模型实验表明，该方法稳定有效且不易陷入局部最优，具有较好的性能。

1  PSO-BP算法的基本原理和应用

1.1  PSO算法的基本原理

粒子群优化算法(Particle swarm optimization, PSO)是由KENNEDY等^[15]于1995年提出的一种模拟鸟类群体觅食行为的智能算法。它利用鸟群中个体对信息的共享机制，使整个群体的运动在问题求解空间中产生无序到有序的演化过程，从而获得最优解。

标准的PSO算法采用“速度-位置”搜索模型，根据粒子对搜索空间适应度值的大小对粒子的优劣进行评价，其中适应度函数根据优化的目标设定。PSO算法首先在一个d维的搜索空间中初始化为一群规模为N的随机粒子，其中每个粒子i (i=1, 2, …, N)在空间中的位置用表示，速度用表示，然后通过迭代来寻找最优解。在每一次迭代过程中, 粒子通过追逐两个极值来更新自己的位置：一个是粒子自身所找到的当前最优解, 这个解称为个体极值；另一个是整个群体当前找到的最优解，这个解称为全局极值。粒子群的进化方程描述如下：

              (1)

                    (2)

式中：t为当前迭代次数；表示当前粒子i的个体极值；表示当前的全局极值；c ₁和c ₂为加速因子，通常在0~2之间取值；rand₁和rand₂为取值在0~1之间的独立随机数；w为惯性权重。迭代过程中，个体极值和全局极值不断更新，最后得到的全局极值就是算法得到的最优解。

1.2  基于混沌惯性权重的PSO算法

w惯性权重是影响粒子群算法的一个重要参数，合适的选择可以使粒子具有均衡的探索和开发能力。早期的粒子群算法采用固定权重，后来很多研究者对其进行了优化，采用了时变的惯性权重模型，史峰   等^[16]给出了常见的几种w惯性权重的计算方法。

                    (3)

                   (4)

          (5)

式中：t为当前迭代次数；T_max为最大迭代次数；w_s和w_e分别为惯性权重的初值和终值。

以上的惯性权重模型均为光滑递减模型，师学明等^[17]提出了一种阻尼振荡的w惯性权重模型，并证明了振荡递减的下降形态能够帮助PSO算法加快收敛速度，跳出局部极值。本文作者参照振荡递减的思想，结合混沌Logistic方程提出一种振荡递减的w惯性权重策略，其具体实现如下：

  t=0, 1, 2, …, n       (6)

                (7)

式(6)为Logistic方程，其中μ是控制参数，t是迭代次数。当、=4时，Logistic方程处于完全混沌状态。

式(7)为混沌振荡的w惯性权重方程，其中为混沌Logistic方程的输出，经过数值实验，混沌初值的选择对最终适应度的影响不大，混沌振荡PSO算法最终均能收敛到最优解附近。当混沌初值为时不同惯性权重w的下降形态如图1所示。

图1  惯性权重的迭代曲线

Fig. 1  Different curves of inertia weight

同时粒子群算法的进化方程更新成：

               (8)

式中：和为采用不同混沌初值时的混沌Logistic方程输出。

1.3  BP神经网络

BP神经网络是一种多层前馈神经网络，其主要特点为信号前向传播，误差反向传播。在前向传递中，输入信号从输入层经隐含层逐层处理至输出层。每一层的神经元状态只影响下一层的神经元状态。如果输出层得不到期望输出，则转入反向传播，根据误差调整网络的权值和阈值，从而使BP神经网络预测输出不断逼近期望输出，其中三层的BP网络结构如图2所示^[18]。

图2  三层BP网络结构

Fig. 2  BP neural network with three layers structure

图2中x₁, x₂, …, x_n是BP神经网络的输入值；y₁, y₂, …, y_m是BP神经网络的预测输出值；w_ij、w_jk是BP神经网络的权值。

网络初始化以后，根据输入向量x，输入层和隐含层的权值w_ij 以及隐含层阈值，计算隐含层节点的输出为

                         (9)

其中f为隐含层激励函数，常使用s型激活函数：

                              (10)

根据隐含层输出H，连接权值w_jk和阈值b，可计算BP神经网络的预测输出为

                          (11)

根据网络预测输出O_k和期望输出T_k之间的预测误差，可求出权值和阈值的更新公式：

            (12)

                          (13)

              (14)

                               (15)

式中：为学习速率。

1.4  混沌振荡PSO优化的BP网络算法

虽然BP 网络具有精确寻优能力，但由于网络模型初始权值和阈值是随机产生的，因此，对网络的收敛性和学习效率都有一定的影响，且容易陷入局部极小值，特别是在电阻率层析成像中输入节点较多时该问题尤为突出。

为了减小这种影响，采用前文所述的混沌振荡PSO算法优化BP算法，其具体的实现步骤如下：

1) 初始化BP神经网络，设定网络的输入层、隐含层和输出层的神经元个数，传递函数，学习函数和训练函数。

2) 根据BP神经网络的结构确定粒子的维数，使得粒子群中每个粒子的维度分量都对应神经网络的一个权值或阈值。然后在粒子的位置区间[X_min, X_max]和速度区间[v_min, v_max]内随机初始化粒子群，并设置惯性权重w的初始值，加速因子c₁、c₂、种群规模和迭代次数。

3) 计算每个粒子的适应度：对每个输入粒子调用步骤(1)中设定的BP神经网络计算出网络的输出值，并根据训练样本的期望输出计算出当前粒子的适应度，直至计算出每个粒子的适应度。粒子i适应度的计算公式为

                       (16)

式中：n为训练集的样本个数；m为神经网络输出神经元的个数；Y_ij为第i个样本的第j个理想输出值；y_ij为第i个样本的第j个实际输出值。

4) 根据适应度更新粒子的个体极值和全局极值，将每个粒子i的适应度f_i分别与个体极值p_id和全局极值p_gd进行比较，如果f_i较小，则更新相应的个体极值p_id和全局极值p_gd，并记录当前最好粒子的位置。

5) 根据进化方程更新粒子的速度和位置：根据本研究提出的混沌粒子群算法，使用新的进化方程更新粒子的速度和位置，并按照粒子的位置区间[X_min, X_max]和速度区间[v_min, v_max]规范超出限定范围的粒子，然后跳至步骤3)开始新一轮的迭代寻优。

6) 当达到迭代次数或目标精度时，PSO算法终止，所得到的全局最优解p_gd为BP神经网络的权值和阈值，使用该值对BP神经网络的测试样本进行求解，得到网络的预测输出。

2  二维电阻率层析模型反演

本文作者针对电阻率层析成像进行了二维PSO-BP非线性反演，其具体的数据采集情况如下：测量采用Wenner-Schlumberger装置，测量电极为37个，极距为1 m，一条测线上共采集306个数据点。

2.1  BP神经网络结构的设计

针对电阻率反演的BP神经网络建模，目前的文献中主要有两种方法，AHMAD等^[14]使用采集视电阻率的水平位置、垂直位置和视电阻率值为输入节点，对应位置的真电阻率值为输出节点，将每次测量的所有数据点设为一个数据集进行训练。其特点是神经网络的结构简单(输入节点数量为3，输出节点数量为1)、训练迅速。但存在两个问题：其一视电阻率和真电阻率一一对应，无法充分反映视电阻率是电场作用范围内地下电性不均匀体的综合反映这一视电阻率的本质特征；其二视电阻率的垂直位置一般是投影位置，与真实位置存在误差。徐海浪等^[7]考虑到上述缺点，将所有测量的视电阻率作为输入节点，所有模型参数作为输出节点，但输入输出节点数量巨大(输入节点990个，输出节点1 062个)，因此，采用3隐含层的神经网络，而且3隐含层的节点分别为200、160、200。如此大规模的神经网络建立不仅需要通过大量的测试来确定隐含层的最优节点数、而且训练和测试需要更多的数据和时间。本文作者综合考虑以上两种建模方式的优缺点，采取以下建模方式：将采集并预处理后的视电阻率数据以垂直方向分组，每组输入数据包含多列(奇数列)相邻垂直视电阻率数据，输出数据为输入数据组中间列数据对应位置的模型电阻率参数，同时允许不同输入数据组在相邻处存在数据重叠。采用以上方法，经过测试，本研究的BP神经网络输入节点数为54，输出节点为10，每次测量可获得多组数据集，同时大大简化BP神经网络的结构。网络的训练数据通过改变不同异常体的位置和数量来获得，共获取40次测量数据，640组数据集。同时，提供32组测试数据，测试数据均未参加网络训练。

为了完成神经网络的设计，还需要确定隐含层的节点数目，隐含层神经元个数对BP神经网络的性能有较大的影响。若隐含层神经元个数较少，则网络无法充分描述输入和输出变量之间的关系；若隐含层神经元个数较多，则将导致网络的学习时间过长，产生过拟合等问题。为了确定BP神经网络中隐含层神经元的个数，本文作者针对不同神经元个数对模型性能的影响进行了仿真，其评价指标为决定系数R²

   (17)

式中：y_i为第i个训练数据的预测值；Y_i为第i个训练数据的理想值；n为训练数据的数量。由于此时的BP神经网络还未经PSO算法优化，初始阈值和权值是随机选取的，每次选取的阈值和权值都将引起R²的波动，所以这里选取的评价指标为对程序运行20次后对应R²的平均值，其计算的结果如图3所示。

图3  隐含层神经元个数对决定系数的影响

Fig. 3  Influence of nodes numbers in hidden layer to R²

由图3可知，当隐含层神经元为44时，决定系数R²的平均值达到最大值，同时，单隐层BP神经网络在进行ERT反演时，随机阈值和权值的选取使得BP神经网络的输出结果不够稳定，需要进行进一步优化。

综上所述，建立电阻率层析成像二维反演的BP神经网络模型，该模型为单隐层BP神经网络模型，输入节点为54，隐含层节点为44，输出节点为10，其他参数的设置参考文献[14]。具体如下：传递函数为‘logsig’(对数S型传递函数)，训练函数为‘trainrp’(弹性BP算法)，学习函数为‘learngdm’(梯度下降动量学习函数)。

2.2  基于混沌振荡的PSO算法性能验证

由于惯性权重w值对PSO算法的性能有较大影响，通过电阻率层析成像所提供的训练数据，对式(3)、(4)、(5)、(7)所提出的不同w采用PSO-BP算法进行了训练。其主要参数如下：c₁=c₂=1.194 45，种群规模为30，进化代数为300。训练的结果如图4所示。

图4  使用不同惯性权重的适应度下降曲线

Fig. 4  Curves of fitness using different inertia weights

由图4和表1可知，w₂和w_c均能较快地收敛，但是w₂陷入了局部极值，w₃和w_c均能够达到较小的适应度，但是w_c的收敛速度更快，计算时间更少。这是由于w_c的混沌振荡特性使得混沌振荡PSO-BP算法在早期能够跳出局部极值，在晚期能够更快地收敛于全局极值，因此验证了该算法的性能。

2.3  模型反演输出

为了验证反演算法的可行性，在两个不同的异常体模型下，使用混沌振荡PSO-BP算法和最小二乘法(RES2DINV软件的反演结果)进行了反演对比。

表1  不同惯性权重的性能比较

Table 1  Comparison on performance using different inertia weights

用于验证的模型1为两个相邻地电体，两个地电体直接的距离为8个电极距。模型1的基本参数如下：采用Wenner-Schlumberger装置，每排含37个电极，17层电阻率数据，电极距为1.0m，围岩电阻率为100 Ω·m，低阻异常体大小为2 m×3 m，电阻率为10 Ω·m；高阻异常体大小为2 m×3 m，电阻率为500 Ω·m，其顶部埋深均为2 m。用该模型的正演视电阻率作为混沌振荡PSO-BP网络的输入，对网络进行反演测试，网络输出的反演结果如图5所示。

从图5的反演结果可以看出，最小二乘法和混沌振荡PSO-BP反演算法均能够较为准确地反映高低阻异常体的位置、形态和电阻值，但混沌振荡PSO-BP反演算法的结果更加精确，细节方面也更加清晰，其结果优于最小二乘法的反演结果。

用于验证的模型2为两个垂直地电体，用来检验反演的垂直分辨率。模型2的基本参数如下：采用Wenner-Schlumberger装置，每排含37个电极，17层电阻率数据，电极距为1.0 m，围岩电阻率为100 Ω·m，低阻异常体大小为1 m×4 m，电阻率为10 Ω·m，其顶部埋深为0.5 m；高阻异常体大小为1.5 m×4 m，电阻率为500 Ω·m，其顶部埋深均为3.5 m。用该模型的正演视电阻率作为混沌振荡PSO-BP网络的输入，对网络进行反演测试，网络输出的反演结果如图6所示。

从图6中的反演结果来看，最小二乘法虽然也能够反映出垂直异常差异，但是无法分辨两个高低阻异常体的形态；混沌振荡PSO-BP反演算法的高低阻异常体的位置准确、形态与间隔清晰、电阻值与实际模型更加接近，其结果优于最小二乘法的反演结果。

表2所列为混沌振荡PSO-BP算法、标准PSO-BP算法和BP算法的性能对比，其衡量的指标为均方误差MSE和决定系数R²。从表3的数据可以看出：1) 混沌振荡PSO-BP算法较标准PSO-BP算法和BP算法具有更低的均方误差值，这是由于PSO算法优化了BP神经网络的权值和阈值，提高了算法的全局寻优性 能，同时，PSO算法采用基于混沌振荡曲线的惯性权重自适应调整，类似于模拟退火的退火过程，能够在早期更快的跳出局部极值，并加快收敛；2) 混沌振荡PSO-BP算法较标准PSO-BP算法和BP算法具有更优的决定系数值，表明其反演的结果与理论数据更加接近，误差波动小，具有较好的泛化性能和较高的稳定性；3) 算法在模型2的反演中性能优于模型1的反演，这是由于本研究神经网络的训练集和测试集是纵向分组的，因此，垂直方向上的反演准确度更高。

图5  模型1的模型及两种方法的反演结果

Fig. 5  Model 1 and inversion results by two methods

图6  模型2的模型及两种方法的反演结果

Fig. 6  Model 2 and inversion results by two methods

表2  3种反演方法的结果比较

Table 2  Comparison of results from three inversion algorithms

3  结论

1) 合理的设计和规划神经网络的结构和训练样本，能够使得训练出来的神经网络基本准确的反映阻率层析成像非线性反演的输入输出特性，取得较好的反演效果。

2) 合理的惯性权重w的选择能够使得PSO-BP算法较快的逼近全局最优解，优于标准PSO-BP和BP算法。

3) 基于混沌振荡的PSO-BP算法避免了对初始模型的依靠和计算偏导数矩阵的问题，具有较强的适应性。

REFERENCES

[1] SHIMA H, SAKAYAMA T. Resistivity tomography: An approach to 2-D resistivity inverse problem[C]// Expanded Abstracts of 57th SEG Annual Meeting. New Orleans: Society of Exploration Geophysicists, 1987: 59-61.

[2] 卢元林, 王兴泰, 王 若, 孙仁国, 王劲松. 电阻率成像反演中的模拟退火方法[J]. 地球物理学报, 1999, 42(S1): 225-233.

LU Yuan-lin, WANG Xing-tai, WANG Ruo, SUN Reng-guo, WANG Jin-song. The simulated annealing method of electrical resistivity tomography[J]. Chinese J Geophys, 1999, 42(S1): 225-233.

[3] SHIMA H. 2-D and 3-D resistivity image reconstruction using crosshole data[J]. Geophysics, 1992, 57(10): 1270-1281.

[4] LOKE M H, BARKER R D. Least-squares deconvolution of apparent resistivity pseudosections[J]. Geophysics, 1995, 60(6): 1682-1689.

[5] ZOHDY A A R. A new method for the automatic interpretation of Schlumbeger and Wenner sounding curves[J]. Geophysics, 1989, 54(2): 245-253.

[6] LESUR V, CUER M, STRAUB A. 2-D and 3-D interpretation of electrical tomography measurements. Part 2: The inverse problem[J]. Geophysics, 1999, 64(2): 396-402.

[7] 徐海浪, 吴小平. 电阻率二维神经网络反演[J]. 地球物理学报, 2006, 49(2): 584-589.

XU Hai-lang, WU Xiao-ping. 2-D resistivity inversion using the neural network method[J]. Chinese J Geophys, 2006, 49(2): 584-589.

[8] MADAN K J, KUMAR S, CHOWDHURY A. Vertical electrical sounding survey and resistivity inversion using genetic algorithm optimization technique[J]. Journal of Hydrology, 2008, 359(1): 71-87.

[9] P C O. Feasibility analysis of the use of binary genetic algorithms as importance samplers application to a 1-D DC resistivity inverse problem[J]. Mathematical Geosciences, 2008, 40(4): 375-408.

[10] LIU B, LI S C, NIE L C. 3D resistivity inversion using an improved genetic algorithm based on control method of mutation direction[J]. Journal of Applied Geophysics, 2012, 87: 1-8.

[11] M J L, ESPERANZA G G. PSO: A powerful algorithm to solve geophysical inverse problems application to a 1D-DC resistivity case[J]. Journal of Applied Geophysics, 2012, 71(1): 13-25.

[12] GAD E Q, KEISUKE U. Inversion of DC resistivity data using neural network[J]. Geophysical Prospecting, 2001, 49(4): 417-430.

[13] MAITI S, ERRAM V C, GUPTA G, TIWARI R K. ANN based inversion of DC resistivity data for groundwater exploration in hard rock terrain of western Maharashtra (India)[J]. Journal of Hydrology, 2012, 464: 294-308.

[14] AHMAD N, SAMSUDIN T, WAN ABDULLAH W A T. Using artificial neural networks to invert 2D DC resistivity imaging data for high resistivity contrast regions: A MATLAB application[J]. Computers & Geosciences, 2009, 35(11): 2268-2274.

[15] KENNEDY J, EBERHART R. Particle swarm optimization[C]// IEEE International Conference on Neural Networks. Washington DC: IEEE Computer Society, 1995, 4: 1942-1948.

[16] 史 峰, 王小川, 郁 磊. Matlab神经网络案例分析[M]. 北京: 北京航空航天大学出版社, 2010.

SHI Feng, WANG Xiao-chuan, YUN Lei. Matlab neural network case study[M]. Beijing: The Beijing University of Aeronautics & Astronautics Press, 2010.

[17] 师学明, 肖 敏, 范建柯, 杨国世, 张旭辉. 大地电磁阻尼粒子群优化反演研究[J]. 地球物理学报, 2009, 52(4): 1114-1120.

SHI Xue-ming, XIAO Min, FAN Jian-ke, YANG Guo-shi, ZHANG Xu-hui. The damped PSO algorithm and its application for magnetotelluric sounding data inversion[J]. Chinese J Geophys, 2009, 52(4): 1114-1120.

[18] 闫永利, 陈本池, 赵永贵, 陈 贇, 马晓冰, 孔祥儒. 电阻率层析成像非线性反演[J]. 地球物理学报, 2009, 52(3): 758-764.

YAN Yong-li, CHEN Ben-chi, ZHAO Yong-gui, CHEN Yun, MA Xiao-bing, KONG Xiang-ru. Nonlinear inversion for electrical resistivity tomography[J]. Chinese J Geophys, 2009, 52(3): 758-764.

 (编辑  李艳红)

基金项目：国家自然科学基金资助项目(41374118)；教育部博士点基金资助项目(20120162110015)

收稿日期：2012-10-24；修订日期：2013-04-25

通信作者：戴前伟，教授，博士；电话：13574816998；E-mail: qwdai@csu.edu.cn