運載火箭的抗干擾分數階控制器設計

程昊宇，董朝陽，王 青

（1.北京航空航天大學航空科學與工程學院，北京100191；2.北京航空航天大學自動化科學與電氣工程學院，北京100191）

運載火箭的抗干擾分數階控制器設計

程昊宇1，董朝陽1，王 青2

（1.北京航空航天大學航空科學與工程學院，北京100191；2.北京航空航天大學自動化科學與電氣工程學院，北京100191）

針對火箭的姿態控制問題，基于自抗擾理論和分數階控制，提出了一種控制器設計方法。通過設計擴張狀態觀測器，充分利用慣組平臺和速率陀螺的復合量測信息對系統的內部和外部干擾進行估計。設計分數階控制器對估計的干擾進行實時補償，實現火箭的姿態控制，通過引入可調參數使得控制器參數調整范圍變大，保證控制器具有更好的控制精度。考慮到控制器參數多且難以被確定，采用一種基于遺傳算法和粒子群優化的混合算法對參數進行調整；通過仿真驗證了所提算法可以獲得良好的動態性能、抗擾性和魯棒性。

姿態控制；擴張狀態觀測器；分數階控制器；遺傳算法；粒子群優化

0 引 言

隨著空間發射任務的多樣化和復雜化，運載火箭的結構越來越復雜，性能要求也越來越高。傳統運載火箭采用比例微分+校正網絡的控制方式［1］，對于不同的飛行狀態，通過改變動靜態增益和校正網絡對火箭進行控制，但是運載火箭在飛行過程中會受到陣風等擾動的影響，導致姿態控制系統性能下降。傳統比例微分（proportional-differential，PD）控制不易滿足系統高性能的要求。為了滿足運載火箭的機動性、靈活性以及衛星高精度入軌的要求，需要設計控制精度高、抗干擾能力強的姿態控制器。

近年來，許多學者對火箭的姿態控制進行了研究［2-4］，文獻［3］采用自適應方法研究了火箭的姿態控制問題；文獻［4］采用二階滑模實現了運載火箭的姿態控制，在保證系統魯棒性的前提下，有效地削弱了抖振現象。但是以上方法都依賴于火箭的精確模型，由于火箭模型復雜，在實際工程中很難獲得精確的數學模型。分數階PDμ控制［59］與滑模等現代控制方法相比，不僅具有傳統PD控制不依賴于系統精確數學模型的優點，而且通過增加可調節參數擴大控制器的參數整定范圍，與傳統單純調節動態增益和靜態增益的控制方法相比，具有更好的控制效果和控制精度，可以很好地滿足系統對控制精度和魯棒性的要求［8］。

運載火箭在飛行過程中會受到各種內外不確定性的影響，這些不確定性會影響運載火箭的控制器性能甚至導致失穩，因此有必要研究對內外不確定性具有強魯棒性的控制方法。自抗擾控制［1014］可以實時估計和補償系統內部不確定動態，具有不依賴于被控對象的精確模型、魯棒性強、精度高和抗干擾能力強的優點，在航空航天領域獲得了廣泛的應用［15-19］，文獻［15］研究了航天器的高度控制問題，將滑模控制與自抗擾控制結合起來設計控制器；文獻［17］針對航天器動力學姿態控制問題，提出了一種二階線性自抗擾控制方法，說明了自抗擾算法具有較強的魯棒性和控制精度。

基于以上分析，本文針對存在內外不確定性的運載火箭模型，將分數階控制和自抗擾控制結合起來，設計了自抗擾分數階控制器。首先，通過擴張狀態觀測器（extended state observer，ESO）對系統存在的內外不確定項進行實時估計。接著結合分數階PDμ控制，設計控制器對火箭進行控制并實時補償干擾。考慮到分數階控制參數多且難以被確定的問題，將粒子群算法和遺傳算法結合起來，對參數進行優化［2023］，并通過改變慣性權值來提高算法收斂速度。通過仿真驗證了所設計的控制方法具有較好的控制精度和較強的魯棒性。

1 模型描述

火箭的姿態控制是通過俯仰、偏航和滾轉3個控制通道分別實現對俯仰、偏航和滾轉方向的姿態進行控制。由于正常飛行時火箭3個通道的耦合程度較小，工程上一般將運載火箭的3個通道視為解耦的，分別進行控制器的設計。以俯仰通道為例，姿態控制系統是由箭體、慣組平臺、速率陀螺、控制器、校正網絡和伺服機構組成閉環控制回路，為了提高控制系統的魯棒性，本文在此基礎上引入了擴張狀態觀測器，設計的火箭姿態控制回路如圖1所示。

圖1 火箭姿態控制系統結構圖

在半速度坐標系中箭體的質心運動方程分量形式為

式中，m為火箭質量；V為火箭速度；θ為彈道傾角；σ為航跡偏航角；FXh、FYh、FZh為存在彈性形變情況下，作用在火箭箭體上所有外力之和在半速度坐標系X、Y、Z 3個坐標軸上的分量。

令箭體坐標系3個坐標軸為箭體的慣量主軸，由動量矩定理可得運載火箭在箭體系下的繞質心運動方程為

式中，JX1、JY1、JZ1分別為火箭箭體繞箭體系X軸、Y軸、Z軸的轉動慣量；MX1、MY1、MZ1為存在彈性形變的情況下，作用在火箭箭體上的外力之和對質心的力矩；ωX1、ωY1、ωZ1為箭體繞質心轉動角速度在箭體坐標系中的投影分量。

運載火箭在飛行過程中會產生彈性形變，引起彈性振動，其中，俯仰方向的彈性振動方程為

式中，ωiγ、ωiφ、ωi為滾轉、俯仰、偏航方向第i階振動圓頻率；為對應方向第i階阻尼比為相應的廣義質量為相應的廣義力為相應的廣義坐標；nγ、nφ、n為各個方向彈性振動考慮的模態個數。

由于控制系統主要考慮的是姿態角偏差的消除，并使箭體按照選定的軌道和給定的程序角飛行，可以近似認為箭體姿態角偏差暫態過程中方程系數是固定不變的，可以采用“固化系數法”對火箭進行建模和分析。

本文以某型火箭俯仰通道為例進行姿態控制系統的設計，其剛體姿態運動方程［1］為

彈性振動方程為

慣組平臺和速率陀螺測量方程為

式中，Δθ為彈道傾角的偏差；Δα為箭體攻角的偏差；δφ為俯仰通道的發動機擺角輸入為俯仰通道的發動機擺角加速度輸入；αw為風攻角；為火箭受到的干擾力為火箭受到的干擾力矩； 為俯仰角的偏差； 和分Δφ Δφgz別為慣組平臺和速率陀螺測量到的輸出信號。

偏航通道控制器設計和滾動通道控制器設計與俯仰通道控制器設計類似，這里不再贅述。

2 擴張狀態觀測器設計

運載火箭在飛行過程中，會受到陣風、參數漂移等多種內外干擾的影響，柔性火箭的控制器設計還需要考慮箭體彈性振動的影響，而彈性振動頻率也具有不確定性，這些不確定因素會使得姿態控制系統的穩定裕度降低，甚至失穩。本文綜合利用慣組平臺和速率陀螺的測量信息，設計擴張狀態觀測器，采用復合誤差對姿態角、姿態角速度和內外干擾進行估計和補償，簡化了設計。定義復合觀測誤差為

式中，κ和λ為測量權重，可以調節來自慣組平臺和速率陀螺兩種測量信息的比例，若κ＝1，λ＝0，則退化為僅采用慣組平臺測量信息的形式，若0＜κ＜1，0＜λ＜1，則表示綜合采用慣組平臺和速率陀螺的測量信息，可以通過改變權重，使得系統收斂速度加快，超調量減小。

設計擴張狀態觀測器為

式中，z1、z2分別為俯仰角和俯仰角速率的估計；z3為內外不確定項的估計；β1、β2、β3為ESO的參數，可以根據觀測效果進行調整；b0是決定補償強弱的“補償因子”，作為可調參數來使用；fal（·）為非線性函數；δ為設定值；μ1、μ2為其參數，可以根據誤差大小來調整控制增益，其表達式為

利用擴張狀態觀測器可以對系統中的擾動進行很好估計，然后反饋至控制器進行補償。

3 分數階控制器設計

分數階微積分是將整數階微積分拓展到實數域，通過引入可調算子，提高控制系統的控制效果和控制精度。其算子定義［89］為

式中，a和t分別為操作算子的上限和下限；α∈R為其階次。

為了提高運載火箭姿態控制系統性能，采用分數階PDμ控制器對運載火箭進行姿態控制，同時對擴張狀態觀測器的擾動估計值進行補償，設計分數階控制器為

由于分數階控制器從理論上講是無限維的，一般情況下不能直接對其進行分析和設計，需要將其轉化為整數階系統，用整數階系統去近似分數階系統，利用整數階系統的分析方法對其進行分析和設計。本文采用文獻［8］中提出的方法對其進行有理化近似，通過引入兩個調整系數來提高擬合精度，在頻率范圍（ωb，ωh）內，將分數階控制器轉化為階次為N的整數階控制器，即

式中，b、d是為了提高頻段兩端擬合精度而引入的調整系數。將分數階傳遞函數表示成整數階傳遞函數的零極點形式為

第k個零點、極點為

構造出分數階微積分算子的連續有理傳遞函數模型為

式中，K＝（ωbωh）α。

4 控制器參數優化設計

為了充分利用遺傳算法較強的全局搜索能力和粒子群算法較快的收斂速度的優點，本文采用分層結構設計優化算法。底層采用遺傳算法，對控制器參數進行全局搜索；然后選取遺傳算法搜索得到的最優解作為頂層算法的初始種群，采用粒子群算法對控制器參數進行精確局部搜索，而且為了提高粒子群算法的收斂速度，采用動態慣性權值對算法進行設計。

本文對分數階控制器和擴張狀態觀測器的參數進行優化，待優化參數為：①分數階控制器的比例項增益Kp、微分項增益Kd以及微分階次μ；②擴張狀態觀測器的κ、λ、β1、β2、β3、δ、u1，u2。數目n＝11，設待優化參數個體為xi（i＝0，1，…，n－1），中間代個體為x－i（i＝0，1，2，…，n－1），參數的適應度函數記為fi（i＝0，1，…，n－1），歷史最優個體及其適應度函數記為Ib和Fb，變異概率記為Pm，最大迭代代數記為G1。采用遺傳算法優化控制器參數的步驟如下：

步驟1 確定染色體編碼方法。染色體（個體）為

結合11組待優化參數的可能取值范圍，采用二進制對參數進行編碼，將χ編碼為相應位數的二進制串（ρk－1ρk－2… ρ0）滿足

式中，k為χ對應的二進制編碼的位數；

步驟2 如果達到最大迭代代數，則算法終止，反之則執行步驟3；

步驟3 通過式（12）選擇出中間代的n個個體。

式中，i＝0，1，…，n－1；m為［0，n－1］之間的隨機整數，且滿足m≠i；

步驟4 對產生的中間個體進行交叉運算，本文中選取交叉概率為100%；

步驟5 對所有個體進行變異操作，變異概率為Pm；

步驟6 計算個體的適應度值，控制器設計的目標為減小回路誤差，抑制系統超調量。因此選取優化目標函數為

式中，e（t）為俯仰角信號與俯仰角指令信號的偏差；h為俯仰角信號的最大值與穩態值的差；ε1和ε2分別為誤差和超調量的加權系數；

步驟7 選出當代最優個體并與歷史最優個體比較，若優于歷史最優個體，則用當代最優個體替換歷史最優個體，對歷史最優個體進行更新，否則，將當前代的一個個體用歷史最優個體替換后繼續尋優；

步驟8 判斷得到的最優解是否滿足停止準則，若滿足，則輸出最優解，否則跳到步驟2繼續進行尋優。

當迭代代數達到G1后，將采用遺傳算法優化后控制器參數的最優解取出，作為上層粒子群算法的初始種群，對遺傳算法得到的最優解進一步進行優化，運行一定代數G2之后，若滿足停止準則，則將粒子群優化的結果作為最優結果輸出，若不滿足停止準則，則用粒子群種群中的q個粒子隨機替換掉遺傳算法子群中的q個粒子，對遺傳算法子群再一次進行優化求解，不斷循環，直到控制器參數滿足停止準則，輸出最優解。采用粒子群算法對控制器參數進行優化時，慣性權值和參數收斂速度密切相關，較大的慣性權值有利于全局搜索，較小的慣性權值有利于算法的局部開發，加速算法的收斂。本文通過設計動態慣性權值來提高粒子收斂速度，更新粒子的位置和速度的表達式為

式中，r1和r2為［－1，1］內相互獨立且服從均勻分布的隨機數；g1和g2為給定常數，分別為其靠近自身局部最優解的權重和靠近全局最優解的權重；vi（j）為粒子i在第j代的速度；li（j）為粒子i在第j代的位置；Pi（j）為粒子i在第j代的最優位置；Pg（j）為第j代全局最優位置。設動態慣性權值的表達式為

從式（16）可以看出，隨著迭代次數的增多，慣性權值減小，而較小的慣性權值可以加快算法的收斂速度。

上層粒子群算法的流程如下：

步驟1 對種群進行初始化，選取底層遺傳算法的最優解作為粒子群算法種群的位置，并隨機產生每個粒子的速度；

步驟2 判斷是否到達最大迭代數，若到達，則算法終止，反之則執行步驟3；

步驟3 計算粒子群算法種群每個粒子的適應度函數值；

步驟4 更新粒子群算法種群的個體極值和全局極值；

步驟5 按照式（14）和式（15）更新粒子群算法中每個粒子的速度和位置；

步驟6 判斷優化得到的最優解是否滿足停止準則，若滿足，則輸出最優解，反之，跳轉到步驟2繼續進行參數優化。

5 仿真驗證

本文以某型號火箭發射過程中相關秒點的系統參數為例，考慮一階彈性振動，對系統進行仿真，標稱模型的參數取值為

采用改進的Oustaloup濾波方法對分數階控制器進行近似，擬合頻率ωb＝0.000 1，ωh＝1 000；階次N＝5。采用基于遺傳算法和粒子群優化的分層算法對分數階PDμ和擴張狀態觀測器的參數進行整定。設采用遺傳算法迭代500代；采用粒子群算法迭代100代；優化得到的分數階控制器參數為Kp＝1.52，Kd＝2.83，μ＝0.89；擴張狀態觀測器的參數κ＝0.86；λ＝0.15；β1＝9.21；β2＝90.32；β3＝0.91；δ＝0.06；u1＝0.26，u2＝0.17。

與運載火箭控制器的傳統離線設計方法類似，本文提出的方法也主要用于離線設計，避免了由于優化算法耗時帶來工程應用上的不便。由于火箭姿態控制系統主要考慮的是姿態角偏差的消除，并使箭體按照選定的軌道和給定的程序角飛行，分別在不存在干擾和不確定性的標稱情況和存在內外干擾的不確定情況下，給定系統一個初始俯仰角偏差，對系統進行仿真，并將本文所提控制算法與傳統動態增益+靜態增益的控制方式以及單純的分數階控制器進行對比：

情況1 標稱情況下，俯仰角存在5°的初始偏差；

情況2 存在風干擾，俯仰通道干擾力，參數向下攝動30%，其中＝0.005 8s－1，αw＝57.25°＝0.007 9 s－2，俯仰角存在5°初始偏差。

在標稱情況下得到的仿真結果如圖2～圖4所示，其中，PD表示采用傳統的動態增益+靜態增益的控制方式；FOC表示采用普通的分數階控制器；DRFOC表示采用本文提出的自抗擾分數階控制器。可以看出，ESO可以很好地跟蹤俯仰角信號，自抗擾分數階控制具有較快的調節時間、較小的超調量和較小的穩態誤差；在穩態條件下角速率信號的振幅也明顯小于傳統的PD控制和單純的分數階控制，具有更好的動態性能和控制精度，可以滿足運載火箭對控制系統精度和動態性能的要求。

圖2 標稱情況下俯仰角響應曲線

圖3 標稱情況下俯仰角速率響應曲線

圖4 標稱情況下ESO跟蹤曲線

在不確定情況下得到的仿真結果如圖5～圖7所示，可以看出，在不確定情況下，系統性能均有所下降，但是由局部放大圖可知，自抗擾分數階控制系統受到擾動的影響更小，性能下降較小，可以保證系統良好的動態性能和控制精度，采用的ESO可以很好地跟蹤俯仰角信號，而采用傳統的PD控制和分數階控制器，在受到外界干擾時，角度存在穩態誤差，角速率信號在穩態也存在振蕩現象，表明了自抗擾分數階控制器具有較強的魯棒性和抗擾性。

圖5 不確定情況下俯仰角響應曲線

圖6 不確定情況下俯仰角速率響應曲線

圖7 不確定情況下ESO跟蹤曲線

6 結 論

本文針對運載火箭的姿態控制問題，設計了分數階自抗擾控制器。設計ESO對系統中存在的內外不確定項進行估計。采用分數階控制器對火箭姿態進行控制并實時補償ESO估計的干擾，考慮到分數階控制器和ESO參數多且難以被確定的問題，將粒子群算法和遺傳算法結合起來，設計分層優化算法對參數進行優化。仿真驗證了所設計的控制方法具有較好的控制精度和較強的魯棒性。

［1］Li J W.Modeling，design and analysis of large strap-on launch vehicle's attitude control system［D］.Changsha：National University of Defense Technology，2011.（李家文.大型捆綁火箭姿態控制系統的建模、設計與分析［D］.長沙：國防科學技術大學，2011.）

［2］Roshanian J，Saleh A R，Jahed-Motlagh M R.On the design of adaptive autopilots for a launch vehicle［J］.Proc.of the Institution of Mechanical Engineers，Part I：Journal of Systems and Control Engineering，2007，221（1）：27-38.

［3］Oh C S，Bang H，Park C S.Attitude control of a flexible launch vehicle using an adaptive notch filter：ground experiment［J］.Control Engineering Practice，2008，16（1）：30-42.

［4］Wang Q，Wang Z，Dong C Y.Attitude control of flexible launch vehicles based on second-order sliding mode［J］.Journal of System Simulation，2009，21（7）：2006-2009.（王青，王昭，董朝陽.一種基于二階滑模的柔性運載火箭姿態控制［J］.系統仿真學報，2009，21（7）：2006-2009.）

［5］Zamani M，Karimi-Ghartemani M，Sadati N，et al.Design of a fractional order PID controller for an AVR using particle swarm optimization［J］.Control Engineering Practice，2009，17（12）：1380-1387.

［6］Huang L L，Zhou X L.The adaptive design of the fractional order PID controller parameters［J］.Systems Engineering and Electronics，2013，35（5）：1064-1069.（黃麗蓮，周曉亮.分數階PID控制器參數的自適應設計［J］.系統工程與電子技術，2013，35（5）：1064-1069.）

［7］Zhang B T，Pi Y G.Integration of fuzzy and sliding mode control based on fractional calculus theory for permanent magnet synchronous motor［J］.Przeglad Elektrotchniczny，2011，87（11）：251-255.

［8］Li G Y，Guo C，Li Y X.Fractional-order control of USV course based on improved PSO algorithm［J］.Systems Engineering and Electronics，2014，36（6）：1146-1151.（李光宇，郭晨，李延新.基于改進粒子群算法的USV航向分數階控制［J］.系統工程與電子技術，2014，36（6）：1146-1151.）

［9］Huang L L，Zhou X L，Xiang J H.Self-adjusting design on parameters of the fractional order PID controller［J］.Systems Engineering and Electronics，2013，35（5）：1064-1069.（黃麗蓮，周曉亮，項建弘.分數階PID控制器參數的自適應設計［J］.系統工程與電子技術，2013，35（5）：1064-1069.）

［10］Han J Q.From PID to active disturbance rejection control［J］.IEEE Trans.on Industrial Electronics，2009，56（3）：900-906.

［11］Li S L，Yang X，Yang D.Active disturbance rejection control for high pointing accuracy and rotation speed［J］.Automatica，2009，45（8）：1854-1860.

［12］Xia Y Q，Dai L，Fu M Y，et al.Application of active disturbance rejection control in tank gun control system［J］.Journal of the Franklin Institute，2014，351（4）：2299-2314.

［13］Guo B Z，Zhao Z.On the convergence of an extended state observer for nonlinear systems with uncertainty［J］.Systems& Control Letters，2011，60（6）：420-430.

［14］Guo B Z，Zhao Z L.On convergence of non-linear extended state observer for multi-input multi-output systems with uncertainty［J］.IET Control Theory&Applications，2012，6（15）：2375-2386.

［15］Xia Y Q，Zhu Z，Fu M Y，et al.Attitude tracking of rigid spacecraft with bounded disturbances［J］.IEEE Trans.on Industrial Electronics，2011，58（2）：647-659.

［16］Fang Y C，Shen H，Sun X Y，et al.Active disturbance rejection control for heading of unmanned helicopter［J］.Control Theory&Applications，2013，30（11）：238-243.（方勇純，申輝，孫秀云，等.無人機航向自抗擾控制［J］.控制理論與應用，2013，30（11）：238-243.）

［17］Kang Y，Li D H，Lao D Z.Performance comparison of active disturbance rejection control and sliding mode control in spacecraft attitude control［J］.Control Theory&Applications，2013，30（12）：1623-1629.（康瑩，李東海，老大中.航天器姿態的自抗擾控制與滑模控制的性能比較［J］.控制理論與應用，2013，30（12）：1623-1629.）

［18］Zhu Z，Xu D，Liu J M，et al.Missile guidance law based on extended state observer［J］.IEEE Trans.on Industrial Electronics，2013，60（12）：5882-5891.

［19］Ma Y Y，Tang S J，Guo J，et al.High angle of attack control system design based on ADRC and fuzzy logic［J］.Systems Engineering and Electronics，2013，35（8）：1711-1716.（馬悅悅，唐勝景，郭杰，等.基于自抗擾與模糊邏輯的大攻角控制系統設計［J］.系統工程與電子技術，2013，35（8）：1711-1716.）

［20］Jiao B，Lian Z，Gu X.A dynamic inertia weight particle swarm optimization algorithm［J］.Chaos，Solitons&Fractals，2008，37（3）：698-705.

［21］Jin M，Lu H X.A multi-subgroup hierarchical hybrid of genetic algorithm and particle swarm optimization［J］.Control Theory&Applications，2013，30（10）：1231-1238.（金敏，魯華祥.一種遺傳算法與粒子群優化的多子群分層混合算法［J］.控制理論與應用，2013，30（10）：1231-1238.）

［22］Zhang Y，Su G S，Yan L B.Cooperative optimization algorithm based on particle swarm optimization and Gaussian process［J］.Systems Engineering and Electronics，2013，35（6）：1342-1347.（張研，蘇國韶，燕柳斌.基于粒子群優化與高斯過程的協同優化算法［J］.系統工程與電子技術，2013，35（6）：1342-1347.）

［23］Ran M P，Wang Q.Spacecraft rendezvous trajectory optimization method based on EPSO［J］.Journal of Astronaustic，2013，34（9）：1195-1201.（冉茂鵬，王青.一種基于EPSO的航天器交會軌跡優化方法［J］.宇航學報，2013，34（9）：1195-1201.）

Disturbance rejection fractional-order controller design for launch vehicle

CHENG Hao-yu1，DONG Chao-yang1，WANG Qing2
（1.School of Aeronautic Science and Technology，Beihang University，Beijing 100191，China；2.School of Automation Science and Electrical Engineering，Beihang University，Beijing 100191，China）

A novel controller based on active disturbance rejection control and fractional-order control is proposed for the attitude control of launch vehicle.An extended state observer is designed to estimate the impact from parametric uncertainties and disturbances by taking full advantage of the information from the gyros and the attitude sensors.A fractional-order controller is proposed to realize the attitude control by compensating the disturbances.A variable parameter is brought to expand the tuning range and obtain higher control accuracy.A hybrid algorithm based on genetic algorithm and particle swarm optimization is proposed to tune the parameters of the controller.Simulation results show that the controller exhibits better dynamic performance，higher ability in disturbance rejection and stronger robustness against the disturbances.

attitude control；extended state observer（ESO）；fractional-order controller；genetic algorithm；particle swarm optimization

V 448.1 文獻標志碼：A DOI：10.3969／j.issn.1001-506X.2015.09.23

程昊宇（1990-），男，博士研究生，主要研究方為飛行器制導與控制、飛行器總體設計、魯棒控制。

E-mail：chenghaoyu＠buaa.edu.cn

董朝陽（1966-），男，教授，博士研究生導師，主要研究方向為飛行器建模、飛行器制導與控制、魯棒控制。

E-mail：dongchaoyang＠buaa.edu.cn

王 青（1968-），女，教授，博士研究生導師，主要研究方向為導航制導與控制、魯棒控制、自抗擾控制。

E-mail：wangqing＠buaa.edu.cn

2014-08-27；

2015-03-01；網絡優先出版日期：2015-04-28。

網絡優先出版地址：http：／／www.cnki.net／kcms／detail／11.2422.TN.20150428.1651.003.html

國家自然科學基金（61273083）資助課題