基于變權重免疫克隆選擇算法的新型力矩電機多目標優化設計

李 勇 趙新華 安躍軍 任明旭（沈陽工業大學電氣工程學院 沈陽 110870）

?

基于變權重免疫克隆選擇算法的新型力矩電機多目標優化設計

李 勇 趙新華 安躍軍 任明旭
（沈陽工業大學電氣工程學院 沈陽 110870）

摘要針對傳統的電機設計方法工作量大、數據處理繁瑣和設計偏差大的不足，提出了新型力矩電機的優化設計方法。該方法利用APDL語言對電機進行參數化建模并進行有限元分析，在此基礎上以Matlab為平臺同Ansys實現數據連通，將目標“角度-轉矩”曲線分解為多個基于差值二次方的目標函數，建立了基于聚合函數的新型力矩電機多目標優化模型。通過分析各個目標函數的權重，同時根據所得曲線中各部分擬合的難易程度，確定變權重機制，與免疫克隆選擇算法（ICSA）相結合，提出了變權重免疫克隆選擇算法（VW-ICSA），用VW-ICSA對上述模型進行求解，對該電機的參數進行優化。仿真結果表明，優化的結果更符合設計的要求，人工處理數據的工作量大大減少。

關鍵詞：新型力矩電機 參數化建模 多目標優化設計 變權重免疫克隆選擇算法

國家自然科學基金（61102124）和遼寧省教育廳科學技術研究一般項目（L2015387）資助。

Multi-Objective Optimal Design of New-Type Torque Motor Based on Variable Weight Immune Colonial Selection Algorithm

Li Yong Zhao Xinhua An Yuejun Ren Mingxu
（College of Electrical Engineering Shenyang University of Technology Shenyang 110870 China）

Abstract Aiming at the problems of the immense data processing and large design error about the traditional motor design method, an optimal design method for the new-type torque motor is proposed. In this method, APDL language is used to perform parameterized modeling and finite element analysis of the new-type torque motor. Matlab is used to realize data transmission with ANSYS software. The objective curve of angle-torque is decomposed into multiple objective functions according to the squared deviations. A multi-objective optimization model of the new-type torque motor is established based on aggregate function. Through analyzing the weight of each objective function and the difficult degree of the curve fitting, the mechanism of weight variation is determined. Combined with the immune clonal selection algorithm (ICSA), variable weight immune clonal selection algorithm (VW-ICSA) is proposed. Then, VW-ICSA is used to identify the proposed model and its parameters. The simulation and experimental results show that the optimization results are more suitable for the design requirements, and the workload of manual data processing is greatly reduced.

Keywords：New-type torque motor, parametric modeling, multi-objective optimization design, variable weight immune colonial selection algorithm

0 引言

操縱負荷系統作為飛行模擬器中重要的組成部分是提供飛行員訓練和工程研究的操縱力人感系統，在電動式操縱負荷系統中，駕駛員操縱力的感覺主要由力矩電機提供[1,2]。由于傳統的力矩電機，轉子電阻高、損耗大以及所產生的熱量也大，特別在低速運行和堵轉時更為嚴重，因此在電機后端蓋上還要裝有同軸風扇或離心式風扇，作強迫通風冷卻用[3,4]。

針對直流力矩電機的不足，設計出的新型力矩電機，定子及轉子主要由永磁體組成，不用電流驅動從而不需要相應復雜的控制電路，既大大降低了熱能的損耗同時又節省了能源的消耗。

新型力矩電機設計運用傳統的方法即運用Ansoft軟件進行建模和磁場的計算，不斷地修改電機各部分的參數使其達到設計要求，該傳統方法工作量比較大，數據處理繁瑣，人工耗時也隨之增多，同設計要求的偏差也較大。

相對于電機傳統的優化設計方法，以全局優化算法為代表的新型優化算法能跳出局部最優解，以較少的計算代價得到最優解。文獻[5]提出了一種基于多種群遺傳算法的環形繞組形式的無鐵心永磁直線同步電機優化設計方法，在建立磁場分析模型的基礎上，推導了關鍵參數的解析表達式。以永磁體、環形繞組的尺寸為變量，以推力體積比、電機常數和推力波動為優化目標，建立了基于權重系數的多目標優化函數。最后將優化的結果通過有限元仿真和樣機實驗測試進行了驗證。文獻[6]利用有限單元法（Finite Element Method, FEM）對電機進行優化設計的參數化網格剖分，并結合表面響應法構建優化模型，再利用粒子群優化算法尋求目標函數的最優結果。文獻[7]采用Ansoft進行有限元磁場分析計算，取代傳統的直線電機數學模型分析計算，結合梯度遺傳算法通過編程自動調用、改寫Ansoft軟件下的電機有限元模型進行優化，使得計算結果更精確。文獻[8]建立了電機的數學模型，結合電機設計的工程實際對粒子群優化算法做出了適當改進，應用改進的算法分別對電機進行了優化計算，并對優化設計前后的性能指標及有效材料用量進行對比分析。

由于新型力矩電機的定子及轉子均為永磁體組成，并且設計要求電機符合特定的“角度-轉矩”曲線，而通過Ansoft處理起來比較繁瑣，優化設計結果偏差大。本文通過Ansys自帶的APDL語言建立電機的有限元模型并進行電磁分析，其分析結果均被調用至Matlab數據庫中，同時將目標“角度-轉矩”曲線分解為多個基于差值二次方的目標函數，建立基于聚合函數的新型力矩電機多目標優化模型，根據所得曲線中各部分擬合的難易程度，確定各個目標函數的變權重機制，最后用提出的變權重免疫克隆選擇算法（Variable Weight Immune Clonal Selection Algorithm, VW-ICSA）求解該多目標模型，優化該電機所選參數。

1 新型力矩電機有限元模型及設計目標

新型力矩電機的模型區別于傳統的電機，電機中沒有繞組，定子及轉子僅由永磁體構成，建模參數見表1。

表1 新型力矩電機主要參數Tab.1 Main parameters of new-type torque motor

所建立的新型力矩電機的三維模型如圖1所示。

圖1 新型力矩電機三維模型Fig.1 3D model of new-type torque motor

該電機結構簡單，具有較好的對稱性，采用Ansys對其二維模型進行電磁分析，分析之前假設磁場在z軸方向沒有變化，并且忽略端部效應，其二維結構如圖2所示。

圖2 電機二維結構簡圖Fig.2 2D structure diagram of motor

模型中定子軛與轉子軛均采用10號鋼，其磁導率按照磁特性曲線變化，定轉子磁極材料為釹鐵硼（NdFeBN50），其中Br=1.425T，Hc=11kOe。

電磁轉矩的計算應用麥克斯韋張量法，通過對位于氣隙封閉表面范圍內的面積積分可得到電機的總轉矩為[9-11]

式中，r為該點指向旋轉軸線的半徑矢量；n為表面S的法線矢量。

根據電機磁力線分布的性質，可以將電機等效成直線部分的橫截面，這樣就能將三維模型簡化成一個二維模型來處理，相應的電磁轉矩方程也得到簡化，原積分為沿一曲面積分，簡化后為沿氣隙周圍封閉線l的線積分，這里假設曲線是以半徑為r的一個圓周，電磁轉矩的表達式變為[11]

式中，lFe為電機鐵心長度。

在有限元電磁分析中只要把積分路徑選擇在定、轉子之間的氣隙中，新型力矩電機的電磁轉矩便可得出。

應用于飛行模擬操縱負荷系統的力矩電機要求其轉矩特性即電磁轉矩與定、轉子相對位置之間符合特定的關系曲線如圖3所示。其具體的函數表達

圖3 總體目標曲線Fig.3 Overall target curve

式為

式中，θ為轉子相對定子轉過的角度，如圖4所示。T理想為轉子在θ角度時電機的轉矩。

圖4 優化變量Fig.4 Optimization variables

圖4中虛框位置為轉子的初始位置，此時θ= 0°，定、轉子上對應的永磁體的軸線重合。同時設定轉子順時針旋轉時θ 為正值，逆時針旋轉時為負值，由于轉子在兩個方向旋轉時對應的轉矩曲線關于原點中心對稱，故僅需考慮θ∈[0°,30°]。

通過APDL語言建立電機模型并且使用循環語句計算出角度范圍內的轉矩，同目標曲線進行比較。

2 新型力矩電機多目標優化模型的建立

2.1 優化變量

考慮到電機設計、永磁體安裝等因素，為使優化設計更接近合理化，給出了設計變量及相應的限制條件。區別于一般永磁電機轉子上各永磁體的圓心角保持相等的設計，本文選取轉子上6塊永磁體的圓心角x=(1α,2α,3α,4α,5α,6α)作為優化變量，如圖4所示。6個變量iα各自獨立取值，取值范圍均為[20°,50°]。另外設定定子永磁體圓心角α= 42°；定子永磁體厚度d1=0.0035m；氣隙長度δ=0.0008m。

2.2 目標函數

將圖5中的目標“角度-轉矩”曲線分解為多個基于差值二次方的目標函數，采用聚合函數法通過線性加權和的方式將多個目標聚合為一個目標，建立目標函數。其值越小說明優化得到的曲線更接近目標曲線，目標函數為其中

式中，n為采樣點個數；w1,w2,…,wn為各個單目標函數的權重；（x）為當優化變量為x，轉角為θk時，用Ansys通過場計算得出的電機轉矩；為式（3）在θ=θk時所得的理想電機轉矩；θk為第k個采樣點對應的當前轉子相對定子轉過的角度，即如圖5所示在[0°,30°]內均勻選取的15個采樣點（2°,4°,6°,8°,10°,12°,14°,16°,18°,20°, 22°,24°,26°,28°,30°）。

圖5 目標“角度-轉矩”曲線Fig.5 Aim curve of angle-torque

2.3 約束條件

由于新型力矩電機不需要電力供應，因此不包括功率因數、起動轉矩、起動電流和電機效率等通常的電機約束條件。該電機優化的約束條件記為gm，m=1～4，g1、g2、g3和g4為優化的約束函數，如式（6）所示約束條件均在要求的變化范圍內。考慮到電機設計、安裝等因素，為使優化設計更接近合理化，同時減小轉子永磁體間尺寸不同產生的影響，選取轉子上三組兩兩相對的永磁體圓心角之差1β=1α-4α、2β=2α-5α和3β=3α-6α，定子軛部磁通密度Bd1為約束條件。

綜上，建立的電機多目標優化模型為

式中，F（ x）為聚合的目標函數；f1（x）,…,fn（x）為n個單目標函數；w1,w2,…,wn為各個單目標函數的權重；x為六維設計變量，xl和xu分別為x的下限和上限；gm（x）為約束條件。

3 基于VW-ICSA的新型力矩電機多目標優化設計

3.1 VW-ICSA中的變權重機制

在非線性程度高的情況下，權重固定的目標函數在全局尋優過程中，其優化的結果與設計要求之間的偏差較大，且隨著優化算法運行代數的增加也不能有效地減小。同權重固定不變相比，變權重機制能夠利用當前種群中的有用信息來重新調整權重，獲得朝向理想點的搜索壓力，使收斂速度更快，得到的最優解更好[12,13]。通過多次仿真計算，式（7）中n個單目標函數的適應性權重及變權重機制如下式中，Ni（x）為第i個采樣點的實際轉矩與理想轉矩的差值，gen為優化算法的進化代數；C為增長系數，其取值根據實際情況而定。

變權重機制如下：當代數gen=1時，設定增長系數C=1，若采樣點上實際轉矩與理想轉矩的差值函數，則變大，從而權重w增

i大；當代數gen＞1時，C為略大于1的某個定值，若差值函數較上一代變化不大，則會呈指數倍增大，從而權重wi會更大，進一步加大搜索壓力，直至差值減小。若差值較上一代大大減小，則的增大幅度也會大大減小，從而權重wi基本保持不變。

根據曲線中各采樣點擬合的難易程度，變權重機制具體分為以下三類：

（1）針對曲線中擬合難度大的點，需要進一步加大搜索壓力，加大該點的權重，增長系數C取大于1的值。

（2）曲線中擬合難度一般的點，則不需要太大的權重，故C的取值為1。

（3）由于方案自身的不足，某些采樣點的值變化幅度很小，因此將其權重固定為1來避免對其他采樣點的擬合產生影響。

采樣點擬合的難易程度通過對固定權重的目標函數的分析獲取。首先保持各個單目標函數的權重相等，調用ICSA進行尋優，經過5次優化之后，各個采樣點上實際轉矩與理想轉矩之差的絕對值見表2。

表2 轉矩差值比較Tab.2 Comparison of torque difference(單位：N·m)

為更直觀地分析表中的數據，圖6將表2中各個采樣點對應轉矩的平均差值連成曲線。圖中的參考線反映了曲線擬合的難易程度，高于參考線的部分，說明擬合難度較大，而低于參考線的部分則說明擬合難度一般。由于采樣點的選取具有代表性，平均轉矩差值曲線能夠大致地反映出采樣點所在的整條“角度-轉矩”曲線同目標曲線之間各部分的差值變化情況。

圖6 平均轉矩差值曲線Fig.6 Curve of average torque difference

由圖6可見，采樣點1由于方案本身的原因平均轉矩差值較大，而采樣點2～6及14、15的平均轉矩差值大部分高于1.5N·m，說明該部分采樣點的擬合難處較大，采樣點7～13的平均轉矩差值則大部分遠低于1.5N·m，其擬合難度一般。由此可以根據具體的變權重機制對其進行設定，采樣點1對應的單目標函數采用機制三固定權重；采樣點2～6 及14、15對應的單目標函數采用機制一加大權重，其中設定C=1.05（增長系數不宜過大，否則會使權重增大過快從而降低了整體的多目標尋優準確性）；采樣點7～13對應的單目標函數則采用機制二來調整權重。

通過采用變權重機制，差值函數減小時相應的權重w也會減小，而差值大的單目標函數的權重主動增加。通過這種方法能使電機優化過程中得到的數據更加準確，收斂速度加快從而能更快地找到最優值。

3.2 優化設計流程

首先確定電機優化模型中的目標函數，選取優化設計的變量，設立約束條件。最后結合Matlab調用Ansys電磁計算的數據利用VW-ICSA進行優化。

具體的優化步驟敘述如下：

（1）設定電機轉子6個永磁體圓心角變量x尋優的范圍，確定種群規模Npop、進化代數gen及終止條件，記憶單元的規模0.2Npop。

（2）隨機產生Npop組優化變量x（抗體）即, i=1,2,…,Npop。

（3）調用Ansys宏，計算出以上Npop組抗體在各個采樣點下的實際轉矩，根據式（5）計算出各個單目標函數值，同時根據式（8）確定各個單目標函數的權重值，最后由式（4）得出Npop組目標函數值F( x)。根據式（7）中的約束條件計算懲罰值，及加入懲罰因子的抗體-抗原親合度函數Feval。在當前種群Npop中給定一個抗體x，構造適應性罰函數為式中，ci（x）為將式（7）中等式與不等式約束轉化為ci（x）≤bi統一形式后的“≤”的左邊部分；Δbi（x ）是當前抗體對第i個約束的違背值；是當前種群中對約束i的最大違背值；ε是一個小正數，用來避免罰函數中出現被零除的情況。加入懲罰因子的抗體-抗原親合度函數為

（4）將種群中的抗體按照其對應的Feval值從小到大的順序排序。選擇前0.15Npop個抗體加上剩余的抗體中隨機選出的0.05Npop個抗體，放入記憶單元中，作克隆操作。設對抗體xi做克隆操作得到yi，用cloneop表示克隆算子則

式中，Ii為元素為1的qi維行向量，qi=0.2Npop，i=1～0.2Npop。

（5）對每個yi進行概率為1的多項式變異得到muti。多項式變異操作用Pmop表示

（6）計算muti及xi的加入懲罰因子的抗體-抗原親合度函數（用步驟（3）的方法計算）。

（7）更新記憶單元。選出Feval（muti）中最小值的個體記為muti，*，比較與Feval（xi）的大小，如果則令xi=muti，否則xi保持不變，i=1～0.2Npop。

（8）代數加1，判斷是否滿足終止條件，若滿足則輸出最小的值Feval（muti）及對應的x，若不滿足則返回（4）。

4 優化結果分析

以Matlab為平臺循環調用Ansys對模型的有限元分析結果以便VW-ICSA優化。表3為VW-ICSA運行3次之后的優化結果。

表3 VW-ICSA最優解和對應的變量Tab.3 Optimal results and the corresponding variables of VW-ICSA

由表3可見VW-ICSA運行3次的結果基本一致，說明ICSA引入變權重機制并沒有影響算法的穩定性。第3次得到的解最小，因此選其為最終的最優解，其所選采樣點的實際轉矩與理想轉矩之間差的絕對值如圖7所示。

圖7 轉矩差值曲線Fig.7 Curve of torque difference

引入變權重機制之后，通過VW-ICSA得到最優解更小，在圖7中，絕大多數采樣點對應的轉矩差值較平均轉矩差值曲線都有不同程度的減小：擬合難度較大的采樣點2～6及14、15采用變權重機制一后其差值大部分降到了參考線以下；擬合難度一般的采樣點7～13采用變權重機制二后，得到的差值比原來更小；因方案原因差值較大且變化幅度小的采樣點1采用變權重機制三后差值也有所減小。由此可以看出采用的變權重機制的VW-ICSA在尋找最優解過程中起到了一定的促進作用。

表4給出了三種方法即傳統方法、ICSA及引入變權重機制的VW-ICSA經過多次運行實驗后得到的最優解。

表4中傳統方法不僅效率低而且得到的結果偏差大，而通過VW-ICSA得到的最優解較ICSA更小。應用三種方法所得到的“角度-轉矩”曲線圖8所示。

表4 三種方法對應的最優解Tab.4 Optimal solution corresponding to three methods

圖8 角度-轉矩曲線對比圖Fig.8 Curve comparison of angle-torque chart

新型力矩電機應用于飛行模擬器的操縱負荷系統中，為飛行員提供操縱力感覺，其性能直接影響飛行員的訓練效果及對飛行模擬器飛行品質的評價，所以要求實際的“角度-轉矩”曲線與目標曲線之間的偏差越小越好。

通過圖8可見，傳統方法所得到的“角度-轉矩”曲線同目標曲線之間的偏差很大，在引入人工進化算法優化之后，兩種曲線之間的偏差大大減小。雖然經ICSA優化后所得的特性曲線中15°～25°之間的部分與目標曲線對應的部分擬合度較高，但是其整條曲線同目標曲線的方均差（也就是F）比采用VW-ICSA得到的要大。

從以上的結果分析可知，通過Matlab與Ansys聯合，基于VW-ICSA的新型力矩電機的多目標優化設計方法與傳統方法相比計算時間少，優化過程人工參與少，減少了工作量，同目標曲線的擬合效果好。同時與采用ICSA相比其所得結果的偏差更小。

5 結論

本文提出了一種以Matlab與Ansys聯合實現數據連通，基于變權重免疫克隆算法的新型力矩電機的優化設計方法，并將優化結果與傳統方法進行了比較。該優化設計的特點如下：

1）采用Ansys參數化建模并進行有限元分析取代傳統的電機數學模型分析計算，使得計算結果更精確，所選參數可以通過編程循環調用。

2）提出的變權重機制同權重固定相比，單目標函數之間的權重會隨著算法的運行而不斷地調整，其中接近于設計要求的單目標函數的權重會變小，偏離要求較大的單目標函數的權重繼續加大。從而加快了算法的收斂速度，節省了大量的不必要的尋優時間，更快地得到最優結果。

3）與ICSA相比VW-ICSA收斂速度更快，能有效地保持種群的多樣性，能夠同時實現全局搜索和局部搜索，有助于防止進化早熟和搜索陷于局部極小值，可以解決復雜的非線性問題。

參考文獻

[1] 王行仁. 飛行實時仿真系統及技術[M]. 北京: 北京航空航天大學出版, 1998.

[2] 盧穎, 王勇亮, 孫方義. 電磁式操縱負荷系統的設計與實現[J]. 現代電子技術, 2011, 34(11): 161-163. Lu Ying, Wang Yongliang, Sun Fangyi. Design and implementation of electromagnetic control loading system[J]. Modern Electronics Technique, 2011, 34(11): 161-163.

[3] 李永祥, 林莘, 徐建源. 高壓斷路器有限轉角永磁電機操動機構兩種定子結構[J]. 電工技術學報, 2010, 25(5): 61-68.

Li Yongxiang, Lin Xin, Xu Jianyuan. Two stator structures of limited angle permanent magnet motors for operating mechanism on high voltage circuit breaker[J]. Transactions of China Electrotechnical Society, 2010, 25(5): 61-68.

[4] 彭兵, 王成元, 夏加寬, 等. 磁動勢法五相永磁力矩電機轉矩分析[J]. 中國電機工程學報, 2012, 32(21): 105-111.

Peng Bing, Wang Chengyuan, Xia Jiakuan, et al. Torque analysis of five-phase permanent magnet torque motors based on magnetic motive force[J]. Proceedings of the CSEE, 2012, 32(21): 105-111.

[5] 李立毅, 唐勇斌, 劉家曦, 等. 多種群遺傳算法在無鐵心永磁直線同步電機優化設計中的應用[J].中國電機工程學報, 2013, 33(15): 69-77.

Li Liyi, Tang Yongbin, Liu Jiaxi, et al. Application of the multiple population genetic algorithmin optimum design of air-core permanent magnet linear synchronous motors[J]. Proceedings of the CSEE, 2013, 33(15): 69-77.

[6] 劉慧娟, 傅為農. 一種參數化網格剖分方法及其在電機優化設計中的應用[J]. 中國電機工程學報, 2013, 32(21): 125-130.

Liu Huijuan, Fu Weinong. A parameterized mesh generation and refinement method for finite element magnetic field computation and its application in optimal design of electric motors[J]. Proceedings of the CSEE, 2013, 32(21): 125-130.

[7] 崔皆凡, 劉娜. 有限元結合并行梯度遺傳算法及其應用[J]. 組合機床與自動化加工技術, 2012(5): 17-20.

Cui Jiefan, Liu Na. The finite element combined with parallel gradient genetic algorithm and its application[J]. Modular Machine Tool & Automatic Manufacturing Technique, 2012(5): 17-20.

[8] 孟大偉, 張羽, 趙成. 粒子群算法在電機優化設計中的應用[J]. 防爆電機, 2011, 46(5): 1-4.

Meng Dawei, Zhang Yu, Zhao Cheng. Application of particle swarm algorithm to motor optimization design[J]. Explosion-Proof Electric Machchine, 2011, 46(5): 1-4.

[9] 崔鵬, 張錕, 李杰. 基于許-克變換的懸浮電磁鐵力與轉矩解析計算[J]. 中國電機工程學報, 2010, 30(24): 129-133.

Cui Peng, Zhang Kun, Li Jie. Calculation of electromagnetic force and torque of suspension electromagnet based on Schwarz-Christoffel transform[J]. Proceedings of the CSEE, 2010, 30(24): 129-133.

[10] 信建國, 夏長亮, 李洪鳳, 等. 基于解析模型的永磁球形電機永磁體優化設計[J].電工技術學報, 2013, 28(7): 87-95.

Xin Jianguo, Xia Changliang, Li Hongfeng, et al. Optimization design of permanent magnet array forspherical motor based on analytical model[J]. Transactions of China Electrotechnical Society, 2013, 28(7): 87-95.

[11] 王群京, 李爭, 夏鯤, 等. 新型永磁球形步進電動機結構參數及轉矩特性的計算與分析[J]. 中國電機工程學報, 2006, 26(10): 158-165.

Wang Qunjing, Li Zheng, Xia Kun, et al. Calculation and analysison configuration parameters and torque characteristies of a novel spherical steppermotor[J]. Proceedings of the CSEE, 2006, 26(10): 158-165.

[12] 李勇, 劉建昌, 王昱. 改進權重自適應GA及冷軋軋制規程多目標優化[J]. 控制理論與應用, 2009, 26(6): 687-693.

Li Yong, Liu Jianchang, Wang Yu. An improved adaptive weight approach GA for optimizing multiobjective rolling schedules in a tandem cold rolling[J]. Control Theory & Applications, 2009, 26(6): 687-693.

[13] 劉磊, 安斯光, Lu Junwei, 等. 高維多目標電磁場逆問題計算的改進多重單目標 Pareto 采樣算法[J].電工技術學報, 2013, 28(6): 9-15.

Liu Lei, An Siguang, Lu Junwei, et al. An improved multiple single objective pareto sampling algorithm applied to many-objective inverse problems[J]. Transactions of China Electrotechnical Society, 2013, 28(6): 9-15.

李 勇 男，1980年生，博士，碩士生導師，研究方向為特種電機的控制與優化設計以及復雜工業過程的控制與優化。

E-mail: liyong.sut@live.com（通信作者）

趙新華 男，1987年生，碩士研究生，研究方向為特種電機及其控制。

E-mail: zhaoxinhua_2008@163.com

作者簡介

收稿日期2014-03-28 改稿日期 2014-12-15

中圖分類號：TM302