多自由度非线性系统动态参数化模型建模方法研究*

刘昊鹏朱云鹏罗忠†韩清凯

（1.东北大学机械工程与自动化学院，沈阳 110819）（2.谢菲尔德大学自动控制与系统工程系，英国 S13JD）（3.大连理工大学机械工程学院，大连 116023）

多自由度非线性系统动态参数化模型建模方法研究*

刘昊鹏1朱云鹏2罗忠1†韩清凯3

（1.东北大学机械工程与自动化学院，沈阳 110819）（2.谢菲尔德大学自动控制与系统工程系，英国 S13JD）（3.大连理工大学机械工程学院，大连 116023）

针对多自由度非线性系统的动态模型辨识问题，基于NARX（Non-linear Autoregressive with Exogenous inputs）模型的建模方法，考虑系统的物理设计参数，建立非线性系统动态参数化模型.首先，根据系统输入、输出数据建立系统不同参数下的NARX模型，并通过EFOR（Extended Forward Orthogonal Regression）算法对不同参数下NARX模型进行修正，以统一辨识得到的系统模型结构.随后，建立NARX模型系数与物理设计参数间的函数关系，得到多自由度非线性系统的动态参数化模型.以单输入、单输出两自由度非线性系统为例，根据数值仿真结果，对系统的动态参数化模型建模过程进行说明.最后，以带非线性涂层阻尼的悬臂梁作为试验对象，建立其动态参数化模型以反映其动力学特性.试验结果表明，非线性系统动态参数化模型能准确预测多自由度非线性系统的输出响应，为非线性系统的分析与优化设计提供了理论基础.

多自由度非线性系统，动态参数化模型，NARX模型，涂层阻尼，悬臂梁

引言

在工程实际中，随着对系统分析与设计的要求不断提高，线性模型往往无法满足实际要求.因此在系统建模过程中，常需要考虑系统的非线性特征［1－3］.例如自动导航仪，桥梁结构及生物医药等复杂非线性系统在分析与设计过程中，都需要考虑其模型非线性特征的影响［4］.

20世纪30年代，Volterra提出了利用Volterra级数描述非线性系统，并基于Volterra级数对非线性系统模型进行辨识［5］.此后，各种形式的模块化非线性模型，如Wiener模型［6］、Hammerstein模型［7］等，均被广泛应用于非线性系统的分析与设计中.作为一种更为普遍的非线性系统模型，Billings等于1981年提出了NARX（Non-linear Autoregressive with Exogenous inputs）模型的基本结构及其辨识方法［8］.此后，Chen与Billings等引入多项式型的NARX模型［9］，分别利用正交最小二乘法（OLS，Orthogonal Least Squares）［10］及向前正交最小二乘法（FROLS，Forward Regression Orthogonal Least Squares）［11］，建立了多项式结构下NARX模型辨识方法.正交最小二乘法的优势在于能够筛选并排序出模型中重要的项，从而更好地反映系统动力学特性；而向前正交最小二乘法则是在其基础上，每一步运算都对模型剩余项进行全面筛选，避免了用最小二乘法估计模型系数可能造成的误差.

NARX模型可用于多数非线性系统建模且模型结构简单，它被广泛应用于工程实际当中.如Chiras等将NARX模型用于飞机燃气轮机的建模上，通过试验验证了由系统输入（燃料消耗），输出（轴转速）建立的NARX模型能很好地反映燃汽轮机的动力学特性，此建模方法也为许多其他非线性发动机的设计分析提供了新思路［12］；Espinoza等则将NARX模型用于电力系统的短期负载预测上，由试验验证了不同负载下建立的NARX模型可对相应系统是否过载进行很好预测，因而在安全性和节约生产成本方面拥有很大优势，但是由于当前模型的预测周期有限，进一步的优化设计还在研究中［13］.作为动态模型的一种，由于NARX模型建立了当前响应与过去输入、输出之间的关系且基于的数据足够多，因而不仅可以更准确地描述系统，还可用来对系统未来某时刻的响应进行预测［14］.但是，因为NARX体现的仅是输入、输出间的关系，其模型系数不具有实际物理意义，因而很难对系统特性作进一步的分析和设计.

为解决这一问题，通常希望在建模过程中，将需要设计的物理参数显示表达于模型中，即建立系统参数化模型.例如：De Hoff等引入一种参数化模型，该模型由多个随系统内部参数变化的子模型组成，利用该建模思想对简单非线性系统（多自由度弹簧阻尼系统）进行建模，通过改变系统的非线性阻尼观测其输出响应，再由进一步的分析设计达到使系统减振的目的，最终验证了该建模方法的准确性［15］；Cyrot-Normand等引入状态仿射模型概念，它由多个不同工况下对应的线性模型组成，随后他们利用该建模思想对复杂非线性系统（发电厂中发电机）进行建模，通过改变系统的工况（发电机转子转速）观测其负载状况，进而对发电机进行优化设计以提升其负载能力［16］.Wei基于NARX模型，针对单自由度颗粒阻尼系统及某种热塑性泡沫单元建立了相应的动态参数化模型［17］，然而由于其模型允许参数变化范围较小，对于具有多频率、多模态的非线性多自由度及连续体系统仍需进一步研究.

考虑到传统模型辨识方法的局限性，本文针对多自由度系统，提出一种非线性动态参数化模型的建模方法，并通过对非线性连续系统的试验验证了该建模方法的正确性，为多自由度非线性系统的分析与设计提供了理论基础，具有实际意义.

1 非线性系统的NARX模型

对于在时域上的离散非线性系统，根据输入、输出数据，其NARX模型可写为：

其中，u（t）、y（t）分别为t时刻系统的输入、输出；nu、ny为输入、输出的时间滞后；F［·］为非线性函数.

考虑非线性函数F［·］为幂型多项式结构，则NARX模型（1）可写为［8］：

其中，n＝nu＋ny；l为非线性多项式的最高阶数；θi1i2…il为模型系数.

一般情况下，系统的NARX模型可通过向前正交最小二乘法辨识得到，即在利用正交最小二乘法的基础上，对NARX模型项进行筛选，以最终得到形如式（2）的模型结构.其具体算法可参考文献［8］，本文仅以图1所示的两自由度非线性系统为例对NARX模型结构进行说明.图中，k1＝k2＝104N／m为线性弹簧，c1＝c2＝20Ns／m为线性阻尼，cnon＝500Ns2／m2为非线性阻尼；系统质量分别为m1＝m2＝1kg.

图1 两自由度非线性系统Fig.1 A 2DOF nonlinear system

系统的运动微分方程为：

若利用NARX模型对图1所示系统建模，输入u（t）采用服从［0，3］上均匀分布的随机激励信号，通过Runge-Kutta仿真得到系统的输入如图2（a）所示，系统中M1的输出如图2（b）实线所示.其中，采样频率为fs＝256Hz，采样点为N＝2000.

根据系统输入、系统中M1输出数据可得系统的NARX模型为：

利用图2（a）所示的输入信号，根据NARX模型（5）可得M1的输出如图2（b）虚线所示.与系统仿真输出结果作比较，可知NARX模型能准确描述系统的动态特性.

图2 系统输入、输出时域图Fig.2 Time domain diagram of system input and output signal

2 非线性动态参数化模型的建模方法

由以上分析可知，尽管NARX模型（5）能够描述系统的输入、输出关系，但是与系统的微分方程（4）相比，NARX模型不含有能够反映系统特性的物理参数，因而很难通过该模型直接对系统进行分析、设计，故需要进一步建立参数化的非线性动态模型.

2.1 动态参数化模型结构

由式（2）可知，NARX模型可简写为：

其中，pm（t）为模型的第m项，它由向量x（t）＝［1，u（t－1），…，u（t－nu），y（t－1），…，y（t－ny）］中各元素或元素间乘积组合构成；θm为模型的第m项系数；M为模型包含的所有可能项数目：

需要注意的是，针对同一设计参数，其不同取值下对应的NARX模型结构不一定完全相同.如对于图1所示的两自由度系统，当非线性阻尼cnon＝400Ns2／m2时，其NARX模型为：

由式（8）可知，当改变设计参数取值时，其模型结构与模型（5）的结构并不完全一致.

为方便建立非线性动态参数化模型，可先通过EFOR（Extended Forward Orthogonal Regression）算法［15］对设计参数的K个不同取值下的NARX模型结构进行统一，即对于第k（k＝1，2，…，K）个NARX模型可表示为：

其中，θk，m为模型系数，表示第k个NARX模型的第m项系数；M0为统一结构下模型的项数，M0≤M.

EFOR算法基于正交最小二乘法的基本思想，通过对设计参数下不同取值所对应的NARX模型进行修正，最终得统一的辨识结构.其具体步骤如下：

步骤1：模型正交化：

将第k个NARX模型以矩阵形式写出，有：

其中：

其中，p＝1＋max（nu，ny），pk，m（n），（n＝p，p＋1，…，N）表示模型项pk，m（t）在采样时刻的采样值.

若Pk为满秩矩阵，由线性代数相关理论知Pk可被分解为［6］：

其中，Ak为M×M的单位上三角矩阵；Wk为（N－p＋1）×M维包含正交列向量［ωk，1，ωk，2，…，ωk，M］的矩阵，

因此，将式（14）带入式（10），得：

其中，Gk＝［gk，1，gk，2，…，gk，M］T.

由施密特正交化可得［16］：

将各参数下的NARX模型正交化后，即可通过EFOR算法得到统一的动态参数化模型.

步骤2：EFOR方法

其中，上标（1）表示第一步运算.

由式（19）可计算m＝1，2，…M时每一项的误差减小率ERR（Error Reduction Ratio）［16］：

其中，误差减小率ERR表征模型项的重要程度，ERR值越大，意味着该项对于描述系统结构越重要.

其中，argmax｛·｝表示使函数最大时自变量的取值，即取最大值时模型（9）的项为第l1项，并令该项为正交型（16）的第一个正交向量ωk，1，即；设有向量αk，i，（1≤i≤M），计αk，1＝ pk，l1.计算第一步误差减小率的平均值为：

（2）以此类推，当进行到第S步时，令m≠l1∩m≠l2∩…∩m≠lS－1，根据已筛选出的S－1个正交向量ωk，1，ωk，2，…ωk，S－1，由式（17）计算得：

由式（18），得：

由式（23）、（24），计算余下每一项的误差减小率，得：

其中，m≠l1∩m≠l2∩…∩m≠lS－1

令：

则第S步的平均误差减小率为：

（3）当算法进行到第M0步时，若信号误差ESR满足：

则算法停止.其中，ρ为设定的误差大小.

将筛选出的项进行排序，最终得正交化型为：

由于经过EFOR得到的NARX模型具有统一的结构，因而式（29）可简写为：

步骤3：统一非线性动态模型结构

根据步骤1、2，在设计参数取值不同的情况下，经EFOR算法可得到统一结构后的模型正交化形式（30）.然而，式（30）并不是最终所要的模型表达式（9）.为此需对式（30）中系数gk，m进行逆变换，即施密特正交化的逆变换［8］：

其中

gk，m转化为θk，m后，其对应的模型项即为pm＝αm.至此，K个NARX模型的结构均相同，表达式均为式（9），故仅需考虑模型系数与设计参数间的函数关系，即可建立动态参数化模型.

本文采用较为简单的多项式关系描述模型系数与设计参数R间的关系，即：

其中，J为多项式变量最高阶；R1，…，Rq表示q个设计参数；βj1，…，jq为多项式系数，可通过最小二乘法计算得出.

设式（34）的多项式向量为：

设：

则式（34）的系数向量可由最小二乘法公式得到［18］：

其中：

由式（38）最终估计出系数β，得到式（34）中模型系数与设计参数间的具体关系式.并最终由式（9）和式（34）建立起系统的动态参数化模型.

2.2 算例

图1所示系统，若设计参数为线性阻尼c2，取c2＝6，12，18，24Ns／m，u（t）仍为服从［0，3］上均匀分布的随机激励信号，用Runge-Kutta仿真得系统在不同阻尼c下的四组输入、输出.

系统最大输入时滞为nu＝1、最大输出时滞为ny＝3，取系统最高阶数l＝3，由式（6）得候选项数.由式（6）可知单个数据集下的NARX模型表达式结构为：

其中，x1（t）＝u（t－1），xk（t）＝y（t－k＋1），k＝2，3，4.

基于这4组输入、输出数据集建立的系统动态参数化模型经EFOR算法辨识，其结果如表1所示.

表1 两自由度系统的辨识结果Table 1 Identification results of 2DOF nonlinear system

如前所述，表1中的AERR为筛选所得各项的平均误差减小率，表征该模型项的重要程度，因而可根据表中AERR大小判定模型项的重要性，最终给出的AERR之和表示经辨识所得模型能够反映系统动力学特性的程度.

基于动态参数化模型且通过EFOR算法统一结构得到的系统表达式为：

其中，θm（c），m＝1，2，…4利用3阶多项式写为：

由最小二乘法可得：

令c＝15Ns／m，对比系统中M1的仿真输出与模型辨识所得输出，结果如图3所示.

图3 随机输入下M1的仿真输出与辨识输出对比图Fig.3 Comparison of M1between simulation output and identification output under the random input

由图3可见，在随机输入下，经动态参数化模型辨识所得系统中M1的输出与仿真输出能很好吻合，验证了所建模型能准确反映系统动力学特性.

在一般工程环境中，系统所受激励通常为简谐激励或多谐波激励［19］，因此，考虑简谐输入：

将u（t）代入已建好的动态参数化模型中，当c＝15Ns／m时，对比系统中M1的仿真输出与模型辨识所得输出，结果如图4所示.

由图4可见，简谐输入下，所建立的动态参数化模型仍能很好地反映系统的动力学特性.

综上分析，多自由度非线性系统的动态参数化建模步骤可概括为：

（1）估计系统的最大输入、输出时滞和最高阶数，确定待选项，建立系统的NARX结构，如式（6）.

（2）将系统各参数下的输入、输出代入对应模型中进行辨识，并通过EFOR算法统一辨识所得模型结构式（9）.

（3）建立NARX模型系数与设计参数间的函数关系，最终建立多自由度非线性系统的动态参数化模型.

（4）模型验证，由动态参数化模型辨识所得输出与仿真输出的吻合程度判定所建模型能否准确反应系统动力学特性.

3 动态参数化模型的试验验证

3.1 带非线性涂层阻尼的悬臂梁振动试验

3.1.1 试验对象及仪器

图4 简谐输入下M1的仿真输出与辨识输出对比图Fig.4 ComparisonM1between simulation output and identification output under the harmonic input

试验对象为带非线性涂层阻尼的悬臂梁，完整测试系统包括LMS测试系统、PCB公司的208C03单向动态力传感器、CA－YD－1160压电式加速度传感器及ESD－805模态激振器［20］.

3.1.2 试验步骤

（1）将加速度传感器粘在悬臂梁顶端；力传感器一端与激振器顶杆通过螺栓连接，另一端粘在悬臂梁固定端；阻尼涂层用3M双面胶粘在悬臂梁表面.传感器与阻尼涂层的布置如图5所示.

图5 带阻尼涂层的悬臂梁试验台Fig.5 Test bench ofCantilever beam with damper coating

图6 测试系统装置图Fig.6 Device diagram of test system

（2）将传感器输出端连接至LMS测试系统，打开测试软件并设置采集参数，设定采样频率为2048Hz.

（3）阻尼层长度设置有50、100、150mm，且每个长度下对应有3个位置，即其上端距悬臂梁顶端分别为20、100、180mm，共9个工况.

（4）接通电源，启动功率放大器，将功率放大器输出端与激振器连接，激振器通过力传感器给悬臂梁固定端以激励，利用LMS测试系统测得九个工况下的输入、输出数据，测试系统装置如图6所示.

3.2 带非线性涂层阻尼的悬臂梁参数化模型

根据上述试验步骤，通过LMS测试系统测得上述9种工况下的输入、输出数据.利用工况1、工况2、工况3、工况4、工况5下的输入、输出数据进行建模，其中，每个工况均取前400个采样点.而工况6下的2000组输入、输出数据作为对所建动态参数化模型的验证.工况1到6的输入、输出时域图如下图示：

图7 工况1下输入、输出时域图Fig.7 Time domain diagram of input and output under Condition 1

图8 工况2下输入、输出时域图Fig.8 Input and output time domain diagram under condition 2

图9 工况3下输入、输出时域图Fig.9 Time domain diagram of input and output under Condition 3

图10 工况4下输入、输出时域图Fig.10 Time domain diagram of input and output under Condition 4

图11 工况5下输入、输出时域图Fig.11 Time domain diagram of input and output under Condition 5

图12 工况6下输入、输出时域图Fig.12 Time domain diagram of input and output under Condition 6

设阻尼片长度为A，粘贴在悬臂梁的位置为B.则各工况下参数如表2所示.

表2 各工况下参数Table 2 Parameters under different conditions

令悬臂梁系统中输入、输出最大时滞nu＝ny＝30，最高阶数l＝2.由5个工况下数据所建的动态参数化模型经EFOR算法对系统进行辨识后，得系统结构表达式为：

其中，θm（m＝1，2，…，14）与设计参数A、B间的函数关系可根据式（34）写为：

其中，参数βm，n的值可根据式（38）由最小二乘法给出.

动态参数化模型（45）具体结构如表3所示，式（46）中参数βm，n值如表4所示.

3.3 试验验证

将工况六下2000组输入、输出数据代入悬臂梁的动态参数化模型中，经辨识所得输出与实际输出的对比结果如图13所示.

如图13所示，经系统的动态参数化模型辨识所得输出与实际输出能较好吻合，从试验的角度验证了所建动态参数化模型能准确反映系统动力学特性.

表3 系统辨识结果Table 3 Identification results of system

表4 参数βm，n值Table 4 Estimates for the parametersβm，n

图13 实际输出与辨识输出对比图Fig.13 Comparisonbetween simulation output and identification output

4 结论

（1）动态参数化模型要求每个参数下NARX模型的结构均一致，因此需通过EFOR算法对不同参数下的NARX模型进行修正，以统一辨识得到的系统结构.

（2）动态参数化模型的建立包含如下三部分：首先估计系统的最大输入、输出时滞及最高阶数，确定其NARX模型结构；其次，将系统各参数下的输入、输出数据代入对应模型中进行辨识，并通过EFOR算法统一辨识所得模型结构；最后，利用最小二乘法拟合得到动态参数化模型系数与设计参数间的函数关系，并可通过所得函数表达式计算出相关设计参数下NARX模型的系数.

（3）本文所建的动态参数化模型系数与设计参数间函数关系式为多项式形式.需要注意的是，当模型设计参数较为敏感且变化很快时，多项式形式便不再适用，因而建立合适的函数表达式将在今后研究中讨论.

（4）由前文所述的建模方法，最终建立起带非线性涂层阻尼的悬臂梁动态参数化模型，不同于以往研究的单自由度简单结构，本文研究对象为多自由度连续体系统.随后通过试验验证了所建模型能准确预测系统的动态响应，为多自由度非线性系统的建模提供了新的方法.

1 李秀英，韩志刚.非线性系统辨识方法的新进展.自动化技术与应用，2004，23（10）：5～7（Li X Y，Han Z G.New development on identification methods of nonlinear systems.Automation Technology and Application，2004，23（10）：5～7（in Chinese））

2 Yao M H，ZhangW.Using the extended Melnikovmethod to study multi-pulse chaotic motions of a rectangular thin plate.International Journal of Dynamics and Control，2014，2（3）：365～385

3 Yang X D，ZhangW.Nonlinear dynamics ofaxiallymoving beam with coupled longitudinal-transversal vibrations.Nonlinear Dynamics，2014，78（4）：2547～2556

4 Åström K J，Eykhoff P.System identification a survey.Automatica，1971，7（2）：123～162

5 Li LM，Billings SA.Frequency domain truncation of Volterra series forweakly non-linear systems subject to harmonic excitation.International Journal of Control，2008，81（8）：1291～1297

6 Wigren T.Recursive prediction error identification using the nonlinear Wiener model.Automatica，1993，29（4）：1011～1025

7 Narendra K S，Gallman PG.An iterativemethod for the identification of nonlinear systems using a Hammerstein model.Automatic Control，IEEE Transactions on，1966，11（3）：546～550

8 Billings S A.Nonlinear system identification：NARMAX methods in the time，frequency，and spatio-temporal domains.John Wiley＆Sons，2013

9 Chen S，Billings S A.Representations of non-linear systems：the NARMAX model.International Journal of Control，1989，49（3）：1013～1032

10 Chen S，Billings S A，Luo W.Orthogonal least squares methods and their application to non-linear system identification.International Journal of Control，1989，50（5）：1873～1896

11 Billings S A，Chen S，Korenberg M J.Identification of MIMO non-linear systems using a forward-regression orthogonal estimator.International Journal of Control，1989，49（6）：2157～2189

12 Chiras N，Evans C，Rees D.Nonlinear gas turbinemodeling using NARMAX structures.IEEE Transactions on Instrumentation and Measurement，2001，50（4）：893～898.

13 Espinoza M，Suykens JA K，Belmans R，et al.Electric load forecasting.Control Systems，IEEE，2007，27（5）：43～57

14 Coca D，Billings SA.Identification of coupled map lattice models of complex spatio-temporal patterns.Physics Letters A，2001，287（1）：65～73

15 De Hoff R L，Rock SM.Development of simplified nonlinear models from multiple linearizations.Decision and Control including the 17th Symposium on Adaptive Processes，1978 IEEE Conference on.IEEE，1979：316～318

16 Van Mien H D，Normand-Cyrot D.Nonlinear state affine identification methods：applications to electrical power plants.Automatica，1984，20（2）：175～188

17 Wei H L，Lang ZQ，Billings SA.Constructing an overall dynamical model for a system with changing design parameter properties.International Journal of Modelling，Identification and Control，2008，5（2）：93～104

18 Worden K，Tomlinson G R.Nonlinearity in structural dynamics：detection，identification and modelling.CRC Press，2000

19 Worden K，Manson G，Tomlinson G R.A harmonic probing algorithm for the multi-input Volterra series.Journal of Sound and Vibration，1997，201（1）：67～84

20郭抗抗，曹树谦.压电发电悬臂梁的非线性动力学建模及响应分析.动力学与控制学报，2014，12（1）：18～23（Guo K K，Cao SQ.Nonlinear modeling and analysis of piezoelectric cantilever energy harvester.Journal of Dynamics and Control，2014，12（1）：18～23（in Chinese））

Received 11 May 2016，revised 30 June 2016.

*The project supported by the National Science Foundation of China（11572082），the Fundamental Research Funds for the Central Universities of China（N150304004，N140301001）and the Key Laboratory for Precision＆Non-traditional Machining of Ministry of Education，Dalian University of Technology（JMTZ201602）

†Corresponding author E-mail：zhluo＠mail.neu.edu.cn

MODELING METHOD ON DYNAMIC PARAMETRICAL MODEL OF NONLINEAR MULTI-DEGREE OF FREEDOM SYSTEMS*

Liu Haopeng1，Zhu Yunpeng2，Luo Zhong1†，Han Qingkai2
（1.School of Mechanical Engineering＆Automation，Northeastern University，Shenyang 110819 China）（2.Department of Automatic Control and Systems Engineering，University of Sheffield，UK，S1 3JD）（3.School of Mechanical Engineering，Dalian University of Technology，Dalian 116023 China）

For the identification problem of Multi-Degree of Freedom（MDOF）nonlinear systems，the dynamic parametrical model of a nonlinear system is presented based on the Non-linear Autoregressive with Exogenous inputs（NARX）model and considering the physical design parameters of interest of the systems.Firstly，the NARX model is established under the condition that different design parameters are considered，where each of the established NARX model is only related to the input and output data.These models are then unified to a common structure by using the proposed Extended Forward Orthogonal Regression（EFOR）algorithm.Secondly，the relationship between the coefficients of each unified NARX model and the corresponding physical design parameters are built in order to obtain the dynamic parametrical model of the MDOF nonlinear system.Furthermore，a Single Input Single Output（SISO）2DOF nonlinear system is taken for a case study to clarify the modeling process of dynamic parametrical model in detail by using the results of numerical simulation.Finally，a cantilever beam with nonlinear coating damper is employed for an experimental validation of the proposed modeling method.The dynamic parametrical model of the beam is established to reflect its dynamic characteristics.The results indicate that the dynamic parametrical model of the non-linear systems can accurately predict the response of the nonlinear systems，which provides the theoretical basis for the analysis and design of nonlinear systems.

MDOF nonlinear system，dynamic parametrical model，NARX model，nonlinear coating damper，cantilever beam

10.6052／1672-6553-2016-031

2016-05-11收到第1稿，2016-06-30收到修改稿.

*国家自然科学基金资助项目（11572082）、教育部基本科研业务费专项资金资助项目（N150304004，N140301001）及大连理工大学精密与特种加工教育部重点实验室研究基金资助项目（JMTZ201602）

†通讯作者E-mail：zhluo＠mail.neu.edu.cn