带有落角约束的间接Gauss伪谱最优制导律

陈琦，王中原，常思江

（南京理工大学能源与动力工程学院，江苏南京210094）

带有落角约束的间接Gauss伪谱最优制导律

陈琦，王中原，常思江

（南京理工大学能源与动力工程学院，江苏南京210094）

针对带有落角约束的末制导问题，提出了一种基于极小值原理和Gauss伪谱法的最优制导律。以期望落角方向为坐标轴定义了落角坐标系，并在其中建立了线性化的导引运动关系方程。将控制系统简化为1阶惯性环节，利用极小值原理得到正则方程，然后引入Gauss伪谱法进行离散，将其转化为代数方程，结合边界条件，推导出最优制导律的解析表达式，无需任何积分过程，避免了求解黎卡提微分方程。仿真结果表明，所提出的算法运算量小，计算效率高，同时也能方便地求解出复杂加权矩阵下的最优制导律，能够在满足落角约束的条件下更快地收敛到落角参考线，并且具有更小的末端需用过载。

兵器科学与技术；落角约束；最优控制；末制导律；极小值原理；间接Gauss伪谱法

0 引言

在导弹实施末端精确打击时，为了最大限度提高毁伤效果，往往希望导弹以最小的脱靶量命中目标的同时，还能以特定的落角攻击目标要害或薄弱的部位。因此，带有落角约束的末制导律目前已逐渐成为国内外学者的研究热点［1-4］。文献［5］利用圆轨迹导引思想，设计出了一种带有落角约束的三维制导律；文献［6］和文献［7］采用模型预测静态规划技术，得到了可满足终端落角要求的非线性次优制导律；文献［8］基于滑模控制技术，研究了落角约束导引律设计问题，并同时考虑了自动驾驶仪1阶动力学延迟的情况；文献［9］提出了一种非奇异终端滑模制导律，避免了传统终端滑模制导律在控制量输入饱和时出现奇异现象。

随着最优控制理论的发展，越来越多的学者将最优控制理论应用于带有落角约束的末制导律的设计中。文献［10］将剩余飞行时间的函数作为性能指标，利用最优控制理论得到了带有落角约束的最优制导律。文献［11］利用线性二次型最优控制理论研究了带落角约束的制导律，并与微分对策理论进行了对比，验证了最优控制理论在求解最优制导律方面具有较大的优势。文献［12］通过引入碰撞三角，将导引运动关系方程进行了线性化，之后依然采用最优控制理论研究了最优制导律的问题。在将末制导问题转化为最优控制问题后，上述研究均采用极小值原理进行求解，通过积分正则方程得到协态变量，然后利用极值条件获得最优解。这种方法虽然精度很高，但推导过程较为繁琐，特别是对于复杂的问题，显示求解正则方程将会变得非常困难。为此，文献［13］通过求解黎卡提微分方程得到了最优制导律，这在一定程度上降低了计算难度，但在一般情况下，要得到黎卡提微分方程的解析解依然是一个不小的挑战。考虑到这一点，文献［14］利用黎卡提微分方程的稳态解作为代替，得到了最优制导律，但这只是一种近似的处理方法，局限性较大。与上述研究不同，文献［15］采用Gauss伪谱法，以实时弹道优化的思想研究了带有落角约束的最优制导律，完全避免了积分正则方程和求解黎卡提微分方程，很大程度上简化了推导过程，并取得了不错的制导效果。但是该制导律需要不断地调用优化算法求解非线性规划问题，计算量较大，求解效率较低，实际应用时有一定的局限性。文献［16］采用解析方法求解了考虑高阶自动驾驶仪情况下的最优制导律问题，但其只考虑了简单的加权矩阵，在某些任务需求中，如果需要设计复杂的加权矩阵，该方法求解起来则会变得较为繁琐。

针对以上问题，本文结合极小值原理和Gauss伪谱法，提出了一种新型带有落角约束的最优制导律，称为间接Gauss伪谱最优制导律。利用极小值原理推导出正则方程，引入Gauss伪谱法将正则方程离散为代数方程，结合边界条件，得到了最优制导律的解析表达式。与其他文献的制导律相比，文中所提出的方法不需要任何积分或优化迭代过程，避免了求解黎卡提微分方程，同时在推导过程中不对加权矩阵的具体形式进行限制，可以很方便地处理具有复杂加权矩阵的最优制导律设计问题，具有较强的灵活性。

1 问题描述

带有落角约束的平面末制导几何关系如图1所示。其中:M和T分别表示导弹和目标；OxIyI为地面坐标系；Txfyf为落角坐标系，其原点与目标固连，Txf轴与期望落角方向重合；θd为期望落角的大小；vM和aM分别为导弹飞行速度及法向过载；θM为地面坐标系下的弹道倾角；R和q分别为弹目距离和视线角；θ为落角坐标系下弹道倾角，表示落角误差；y为导弹和Txf轴的距离偏差；为导弹和Tyf轴的距离偏差；为弹目连线和Txf轴的夹角。

图1 导弹与目标几何关系图Fig.1 Geometrical relationship between missile and target

图1中用于表示角度的箭头为逆时针方向时，相应的角度定义为正，反之为负。由此可得到如下的角度关系:

当存在落角约束时，在落角坐标系Txfyf下建立相应的导引运动关系可以简化问题的求解难度。为此，结合图1可得如下的导引运动方程:

假设导弹的飞行速度vM为常数，落角坐标系Txfyf下的弹道倾角θ较小，并定义v=vMθ，则（2）式可线性化为

进一步将导弹控制系统简化为1阶惯性环节，

式中:aM，c和Tc分别表示为制导指令和动力学系统时间常数。结合（3）式和（4）式，可得落角坐标系下的导引运动方程为

式中:x=［yvaM］T；x0=［y0v00］T；

根据图1可知，当导弹到达目标时，如果速度vM和Txf轴重合，即y（tf）=0，v（tf）=0，其中tf为末端时刻，那么即可满足终端零脱靶量和期望的落角约束。因此，本文带有落角约束的末制导问题可以描述为如下的有限时间最优控制问题P:在［t0，tf］时间内，确定控制量u，使得导弹M在满足约束（5）式的条件下，从初始状态x0转移到Txf轴上，并使如下的性能指标最小。

式中:F、Q和R均为正定的对角矩阵，分别表示对末端状态量、状态量及控制量的加权。

2 最优制导律设计

针对第1节中的最优控制问题P，构造如下形式的Hamiltonian函数:

式中:λ（t）为协态向量。根据正则方程可得

根据极值条件

可得最优控制量为

将（10）式带入（8）式，并根据横截条件，可得如下的两点边值问题:

当加权矩阵Q、R比较简单时，可以通过解析方法获得问题P的最优解［10，16］，但是对于复杂的加权矩阵，解析法需要重新推导演算，计算过程较为繁琐，实际应用有诸多不便。此外，反向积分黎卡提微分方程的方法也可以获得最优控制量，但这种方法耗时较长，实际使用时也存在一定的局限性。为此，本文提出了一种间接Gauss伪谱法，通过全局多项式逼近（11）式中的λ（t）和x（t），引入微分矩阵近似相应的导数项，进而将（11）式转换为一系列的代数方程，求解这些代数方程便可得到解析形式的最优控制量。采用Gauss伪谱法需要将时间区域［t0，tf］转换到［-1，1］上，为此引入变量子对时间t进行变换，即

因此，（11）式变为

由于（13）式中第1式已知初值，第2式已知末值，因此采用Gauss伪谱法对二者进行离散的方式略有不同。对于x（子），用N个Gauss节点子1，子2，…，子N和初始点子0=-1上的离散状态构造Lagrange插值多项式近似x（子），

对（14）式微分得

结合（16）式和（17）式，（13）式的离散格式为

式中:xk≡x（子k）；λk≡λ（子k）.根据Gauss求积公式［17］，状态量的末值x（1）和协态量的初值λ（-1）可通过下式得到:

式中:wk为Gauss求积系数。

证明 引入常数变量p（t）=c，利用（16）式对其进行离散可得

接下来引入如下形式的分步积分:

式中:tk（k=1，2，…，N）为Gauss节点。结合（16）式和（17）式，将（23）式中的导数项用微分近似矩阵表示，可得

由于证明过程中并未限制f（t）和g（t）的具体形式，因此（24）式对所有的多项式均成立。取f（t）为N次Lagrange插值基函数Ll（t），l=0，1，…，N（其所使用的N+1个节点为初始点t0=-1加上N个Gauss节点）；g（t）为N次Lagrange插值基函数，j=1，2，…，N+1（其所使用的N+1个节点为N个Gauss节点加上末端点tN+1=1）.以上的取法使f（t）和g（t）各有N+1种形式，如果只使用其中的l=1，2，…，N对应的Ll（t）和j=1，2，…，N对应的，则有

因此可得

因此有

由定理可知，通过微分插值基函数Lk（t）得到矩阵D后，由（20）式便可非常方便地计算出其他3个矩阵

式中:

将（25）式写成简洁形式为

给定ΛN+1和X0，通过（28）式可以解析地求出N个Gauss节点上的协态量λi和状态量xi，i=1，2，…，N，将其带入到（19）式中便可得到端点上的协态量λ（-1）和状态量x（1）.至此，所有节点上的协态量和状态量均已得出，根据（10）式便可得到每个节点上的最优控制量，如（29）式所示，

从以上的推导过程可以看出，最优控制量的计算中不需要任何的积分或迭代过程，只要已知初始偏差X0，即可解析地确定［t0，tf］时间历程内各节点上的控制量，相比于传统的反向积分黎卡提微分方程的方法，这可在很大程度上提高运算速度。由于（29）式中的最优控制量只依赖于初始偏差，所以（29）式得到的是开环最优解，无法抑制外界干扰的影响。为了解决这一问题，本文引入了滚动时域方法。流程如下:

步骤1 初始化Gauss节点和初始时刻t0，给定初始状态量x0.

步骤2 根据当前状态估计剩余飞行时间tgo，得到tf=t0+tgo.

步骤4 只取第1个控制量，即t0时刻对应的控制量u0，将其作用于（5）式中得到下一时刻t′以及状态量x.

步骤5 判断是否命中目标。如果没有命中目标，以当前时刻t′作为初始时刻t0，当前状态作为初始状态x0，重复步骤2～步骤4，直至命中目标。

值得注意的是，（25）式中每个微分矩阵的计算都需要tf的值，因此，步骤2通过估计剩余飞行时间tgo，对tf进行了近似，结合图1，tgo采用下式进行估计:

3 算例仿真与分析

为了验证所提出算法的准确性并考察其性能，3.1节将（28）式的计算结果和目前成熟的最优控制问题求解软件GPOPS的结果进行了对比。3.2节在不同条件下对所提出的末制导算法进行了仿真分析。3.3节将所提出的末制导算法和其他带有落角约束的最优制导律进行了对比分析。

3.1 准确性验证

结合极小值原理和Gauss伪谱法，（28）式给出了最优控制问题P的解析解（包括状态量和协态量），有必要对其准确性进行验证。GPOPS是由A. Patterson和Anil V.Rao等学者开发的最优控制问题通用求解软件，该软件利用Gauss积分伪谱方法将连续时间的最优控制问题直接离散为非线性规划问题，然后调用相关优化算法进行求解。目前该软件广泛应用于轨迹优化、交会对接以及轨道转移等领域，并且其计算精度已被诸多学者所验证［18-20］。因此，3.1节选用GPOPS求解最优控制问题P，并与（28）式的计算结果进行比较。问题P的计算参数假设如下:t0=0 s；tf=30 s；x0=［500，50，0］T；离散节点数目取为40；时间常数假设为Tc=1.0 s；性能指标中的加权矩阵取F=diag（1000，1 000，1 000），Q=diag（1/4002，1/1502，1/502），R=diag（1/502）.对比结果如图2和图3所示，其中x和λ表示（28）式的计算结果，x*和λ*表示GPOPS的计算结果。

图2展示了两种方法得到的状态量和协态量的对比效果，从中可以看出，二者的计算结果吻合得非常好。图3为状态量和协态量误差曲线，从中可以明显地看出，状态量的误差不超过2×10-5，协态量的误差则更小，低于1×10-9.以上对比结果表明，解析表达（28）式具有很高的计算精度。此外，GPOPS在求解过程中必须调用SNOPT或IPOPT进行迭代寻优，这一过程会消耗较多的计算资源并产生较长的计算耗时，而（28）式则避免了这种寻优过程。表1对比了GPOPS和（28）式在MATLAB平台下的计算耗时，可以看出（28）式在计算效率上具有非常明显的优势，这种优势在滚动时域过程中将更为突出。

图2 状态量和协态量对比结果Fig.2 Comparisons of states and costates

表1 两种方法的计算耗时对比Tab.1 Comparison of computation times of two methods

图3 状态量和协态量误差曲线Fig.3 Errors of states and costates

3.2 不同条件下末制导性能

为了考察所提出算法的性能，对不同条件下的末制导性能进行了仿真分析。具体仿真参数见表2，其他参数（离散节点数目、时间常数及权重矩阵）与3.1节一致。仿真结果如图4～图6所示。

表2 末制导仿真参数Tab.2 Simulation parameters for terminal guidance

图4给出了在θ0=0°，θd=-60°条件下，不同起始距离偏差y0对应的弹道曲线及过载曲线。图4中的落角参考线与图1中的Txf轴重合，也表示零过载弹道曲线。由图4可以看出，当导弹沿着落角参考线飞行时，便可满足落角约束。从图4（a）可以看出，在不同的起始距离偏差下，导弹均能较好地趋向于落角参考线，以期望的落角去攻击目标。从图4（b）可以看出，制导初始段过载较大，并且其幅值随y0的增加而增加，这样可以使得落角更快地向期望值收敛。同时还可以看出，不同情况下的过载最终都会趋向于0，这从另外一方面验证了导弹在弹道末端将和落角参考线重合，并且在弹道末端随着需用过载的减小，导弹的抗干扰能力也会随之增强。

图4 不同起始距离偏差情况下仿真结果Fig.4 Simulation results at different initial miss distances

图5展示了不同初始偏角θ0情况下的制导效果（y0=500 m，θd=-60°）。由图5可以看出，导弹能够很好地趋向于落角参考线，可以在满足落角约束的条件下攻击目标。图6给出了在不同期望落角θd条件下制导效果。结合图6（a）和图6（b）可知导弹能够准确地命中目标，并且可以很好地满足相应的落角约束。从图6（c）可看出，随着期望落角的绝对值增加，弹道过载的幅值相应地变大，但最终均也逐渐收敛到0.

以上的仿真算例选取的控制量权函数为常数，然而在某些情况下，根据具体的任务需求，有时会选取时变的控制量权函数，如文献［10］构造的权函数为剩余时间的幂函数，文献［21］研究了一般加权函数情况下的最优制导律求解问题。复杂的权函数给问题的求解带来一定的难度，然而在本文所提方法的推导过程中，并未对权函数的具体形式有所限制，因此，本文方法可以非常容易地处理复杂的加权函数。为了验证所提方法对复杂加权函数的求解效果，参考文献［21］，选择指数函数形式的加权函数，即

式中NG为制导系数。这种形式的加权函数给问题的求解带来了很大的难度，传统的方法将会遇到很大的困难。而本文方法只需在程序中更改W值即可，其他求解公式无需作任何改变，仿真结果如图7所示。

图5 不同起始偏角情况下仿真结果Fig.5 Simulation results in the case of different initial heading errors

图6 不同期望落角情况下仿真结果Fig.6 Simulation results at different desired impact angles

从图7（a）可以看出，随着系数NG的增加，导弹将更快地趋向于落角参考线。图7（b）则展示了导弹过载随NG的变化情况，从中可以看出，随着NG的增加，开始阶段的过载逐渐增加，但末段的过载则更快地趋向于0，这有助于提高导弹的末端过载裕量。不同的加权函数可以产生不同的制导效果，针对特定的任务需求，可以合理地设计相应的加权函数。而本文方法则很容易得到复杂加权函数下的最优控制量，相比与传统的方法，具有很强的灵活性。

3.3 与其他算法的比较

为了进一步验证所提出算法的有效性，将所提算法分别和文献［15］及文献［22］中提出的算法进行对比。

文献［15］以期望落角作为末端约束，能量最优作为性能指标，采用实时弹道优化的思想求解了带有落角约束的末制导问题，并取得了不错的效果。为便于对比描述，本文将这种算法记为弹道优化制导律。

图7 不同过载加权函数情况下仿真结果Fig.7 Simulation results in the case of different acceleration weights

文献［22］所提出的偏置比例制导律算法可描述如下:

式中:uBPN为制导指令，即为导弹过载；η为制导系数；

vM为导弹飞行速度；为视线角速度；θd为期望落角。根据文献［22］，制导系数选取为NG=4，η=1.3.

各算法的对比结果如图8所示。从图8（a）可以看出，3种算法均能命中目标，但是本文方法能以更快的速度趋向于落角参考线。图8（b）的过载曲线展示出本文方法在制导初始段需要较大的过载，这样可以使得导弹更快地收敛于落角参考线，这也和图8（a）中的现象是一致的。此外在弹道末端，偏置比例制导律产生了很大的过载，弹道优化制导律次之，而本文算法则最小。由此可见，本文方法在末端具有较大的过载裕量。由图8（c）可知，3种算法均能保证落角约束，但显然本文方法能更快地收敛于期望落角。

图8 不同制导律情况下仿真结果Fig.8 Simulation results in the case of different guidance laws

4 结论

本文提出了一种带有落角约束具有解析形式的间接Gauss伪谱最优制导律，并在不同条件下进行了仿真验证。

1）在落角坐标系中建立了线性化的导引运动方程，简化了问题的求解，运用最优控制理论将末制导问题转换为最优控制问题。

2）结合极小值原理和Gauss伪谱法推出了解析形式的最优制导律，避免了间接法中的求解黎卡提微分方程和直接法中的寻优过程，在很大程度上提高了末制导效率，同时可以处理一类具有复杂加权矩阵的最优制导律设计问题，可为此类问题的求解提供了一个新的思路。

3）通过和GPOPS结果的对比，验证了所得解析表达式的准确性。不同的仿真算例验证了所提出的算法的有效性。和已有文献［15，20］相比，所提出的算法能够在满足落角约束的条件下更快地收敛到落角参考线，并且具有更小的末端需用过载。

（

）

［1］ Ashwini R，Debasish G.Impact angle constrained guidance against nonstationary nonmaneuvering targets［J］.Journal of Guidance，Control，and Dynamics，2010，33（1）:269-275.

［2］ Sun M W，Xu Q，Du S Z，et al.Practical solution to impact angle control in vertical plane［J］.Journal of Guidance，Control，and Dynamics，2014，37（3）:1022-1027.

［3］ 张友安，马培蓓.带有攻击角度和攻击时间控制的三维制导［J］.航空学报，2008，29（4）:1020-1026. ZHANG You-an，MA Pei-bei.Three-dimensional guidance law with impact angle and impact time constraints［J］.Acta Aeronautica et Astronautica Sinica，2008，29（4）:1020-1026.（in Chinese）

［4］ Harrison G A.Hybrid guidance law for approach angle and timeof-arrival control［J］.Journal of Guidance，Control，and Dynamics，2012，35（4）:1104-1114.

［5］ 胡锡精，黄雪梅.具有碰撞角约束的三维圆轨迹制导律［J］.航空学报，2012，33（3）:508-519. HU Xi-jing，HUANG Xue-mei.Three-dimensional circular guidance law with impact angle constraints［J］.Acta Aeronautica et Astronautica Sinica，2012，33（3）:508-519.（in Chinese）

［6］ 魏鹏鑫，荆武兴，高长生.具有落角约束的弹道导弹再入末制导律设计［J］.哈尔滨工业大学学报，2013，45（9）:23-30. WEI Peng-xin，JING Wu-xing，GAO Chang-sheng.Design of the reentry terminal guidance law with constraint of impact angle for ballistic missile［J］.Journal of Harbin Institute of Technology，2013，45（9）:23-30.（in Chinese）

［7］ Oza H B，Padhi R.Impact-angle-constrained suboptimal model predictive static programming guidance of air-to-ground missiles［J］. Journal of Guidance，Control，and Dynamics，2012，35（1）: 153-164.

［8］ Lee C H，Kim T H，Tahk M J.Design of impact angle control guidance laws via high-performance sliding mode control［J］.Journal of Aerospace Engineering，2013，227（2）:235-253.

［9］ Kumar S R，Rao S，Ghose D.Nonsingular terminal sliding mode guidance with impact angle constraints［J］.Journal of Guidance， Control，and Dynamics，2014，37（4）:1114-1130.

［10］ Ryoo C K，Cho H，Tahk M J.Time-to-go weighted optimal guidance with impact angle constraints［J］.IEEE Transactions on Control Systems Technology，2006，14（3）:483-492.

［11］ Shaferman V，Shima T.Linear quadratic guidance laws for imposing a terminal intercept angle［J］.Journal of Guidance，Control，and Dynamics，2008，31（5）:1400-1412.

［12］ Cho H，Ryoo C K，Tsourdos A，et al.Optimal impact angle control guidance law based on linearization about collision triangle［J］. Journal of Guidance，Control，and Dynamics，2014，37（3）: 958-964.

［13］ 窦磊，杨新民.大着地角卫星制导炸弹最优制导律研究［J］.南京理工大学学报:自然科学版，2010，34（3）:314-318. DOU Lei，YANG Xin-min.Optimal guidance law of satellite guided bombs with large landing angle［J］.Journal of Nanjing U-niversity of Science and Technology:Natural Science，2010，34（3）:314-318.（in Chinese）

［14］ 侯明善.精确对地攻击姿态约束最优末制导设计［J］.兵工学报，2008，29（1）:63-67. HOU Ming-shan.Optimum terminal guidance for air-to-ground missile with impact angle constraint［J］.Acta Armamentarii，2008，29（1）:63-67.（in Chinese）

［15］ Zhang L M，Sun M W，Chen Z Q，et al.Receding horizon trajectory optimization with terminal impact specifications［J］.Mathematical Problems in Engineering，2014（5）:604705.

［16］ Rusnak I，Meir L.Modern guidance law for high-order autopilot［J］.Journal of Guidance，Control，and Dynamics，1991，14（5）:1056-1058.

［17］ 关冶，陆金甫.数值分析基础［M］.北京:高等教育出版社，2005. GUAN Ye，LU Jin-pu.Fundamentals of numerical analysis［M］.Beijing:High Education Press，2005.（in Chinese）

［18］ Darby C L，Rao A V.Minimum-fuel low-earth-orbit aeroassisted orbital transfer of small spacecraft［J］.Journal of Spacecraft and Rockets，2011，48（4）:618-628.

［19］ Furfaro R，Cersosimo D，Wibben D R.Asteroid precision landing via multiple sliding surfaces guidance techniques［J］.Journal of Guidance，Control，and Dynamics，2013，36（4）:1075-1092.

［20］ Boyarko G，Romano M，Yakimenko O.Time-optimal reorientation of a spacecraft using an inverse dynamics optimization method［J］. Journal of Guidance，Control，and Dynamics，2011，34（4）: 1197-1208.

［21］ 张友安，黄诘，孙阳平.带有落角约束的一般加权最优制导律［J］.航空学报，2014，35（3）:848-856. ZHANG You-an，HUANG Jie，SUN Yang-ping.Generalized weighted optimal guidance laws with impact angle constraints［J］. Acta Aeronautica et Astronautica Sinica，2014，35（3）:848-856.（in Chinese）

［22］ Kim B S，Lee J G，Han H S.Biased PNC law for impact with angular constraint［J］.Aerospace and Electronic Systems，1998，34（1）:277-287.

Optimal Guidance Law with Impact Angle Constraints Based on Indirect Gauss Pseudospectral Method

CHEN Qi，WANG Zhong-yuan，CHANG Si-jiang
（School of Energy and Power Engineering，Nanjing University of Science and Technology，Nanjing 210094，Jiangsu，China）

A novel optimal guidance law is proposed for the terminal guidance with impact angle constraints by using the combination of the minimal principle and Gauss pseudospectral method.An impact angle coordinate system is defined with an coordinate axis in the direction of the desired impact angle，and the linear engagement kinematics is established using this coordinate system.The control system of missile is simplified into a first-order inertial system.The canonical equation is obtained via the minimal principle，and then translated into a set of algebraic equations by employing the Gauss pseudospectral method.According to the boundary conditions，an analytical solution is finally derived for the optimal guidance law with impact angle constraints without any integral process or solving the Riccati differential equation.Numerical simulations show that the proposed guidance law ensures the much fast convergence of impact angle to the reference line，and has smaller required terminal acceleration compared with other guidance laws.In addition，the proposed guidance law can easily tackle with the guidance problem with complex weighting matrices.

ordnance science and technology；impact angle constraint；optimal control；terminal guidance；minimal principle；indirect Gauss pseudospectral method

TJ765.3

A

1000-1093（2015）07-1203-10

10.3969/j.issn.1000-1093.2015.07.008

2014-09-09

国家自然科学基金项目（11402117）

陈琦（1989—），男，博士研究生。E-mail:qiychan@126.com；王中原（1958—），男，研究员，博士生导师。E-mail:zywang@njust.edu.cn