郭健彬,杜紹華,王 鑫,曾聲奎
(1.北京航空航天大學(xué)可靠性與系統(tǒng)工程學(xué)院,北京100191;2.可靠性與環(huán)境工程技術(shù)重點實驗室,北京100191;3.北京環(huán)境特性研究所,北京100854)
動態(tài)系統(tǒng)故障的混雜傳播特征及建模方法
郭健彬1,2,杜紹華1,王 鑫3,曾聲奎1,2
(1.北京航空航天大學(xué)可靠性與系統(tǒng)工程學(xué)院,北京100191;2.可靠性與環(huán)境工程技術(shù)重點實驗室,北京100191;3.北京環(huán)境特性研究所,北京100854)
動態(tài)系統(tǒng)的故障傳播過程是由離散事件、連續(xù)特性及其相互作用共同驅(qū)動的,具有顯著的混雜特征,為故障規(guī)律認(rèn)知與建模帶來了較大的難度。現(xiàn)有研究將故障視為離散隨機事件,分析由單元隨機故障發(fā)生到系統(tǒng)失效的離散故障過程,卻忽略了連續(xù)特性對故障傳播的影響,本質(zhì)上是對故障混雜傳播的工程簡化處理,不能真實地描述動態(tài)系統(tǒng)的故障規(guī)律。首先在分析動態(tài)系統(tǒng)故障規(guī)律的基礎(chǔ)上,給出了離散與連續(xù)雙維度下的動態(tài)系統(tǒng)故障混雜傳播的定義,分析了其混雜影響要素以及混雜傳播特征;其次,為了完整準(zhǔn)確描述混雜特征,借鑒混雜理論在描述離散事件與連續(xù)參數(shù)相互作用方面的優(yōu)勢,提出了基于隨機混雜自動機(stochastic hybrid automata,SHA)的故障混雜傳播建模方法;最后通過對某溫度控制系統(tǒng)的故障混雜傳播過程進行建模和仿真,驗證了動態(tài)系統(tǒng)故障傳播過程中的混雜特征,以及建模方法的可行性。
故障傳播;混雜;可靠性;動態(tài)系統(tǒng);隨機混雜自動機
動態(tài)系統(tǒng)故障傳播過程是由離散事件(如故障、重構(gòu)等)、連續(xù)特性變化(如應(yīng)力改變、物理參數(shù)退化、性能波動等)以及兩者之間的相互作用共同決定的。例如,單元的離散故障會導(dǎo)致系統(tǒng)構(gòu)型變化或余度減少,進一步引起應(yīng)力的重新分布或者關(guān)鍵性能指標(biāo)的突變;而系統(tǒng)應(yīng)力(如機械應(yīng)力、熱應(yīng)力、電應(yīng)力、流體力等)的動態(tài)變化以及物理參數(shù)(如磨損、裂紋等)的持續(xù)退化又對故障產(chǎn)生加速或減速作用,當(dāng)連續(xù)特性超過閾值時會直接導(dǎo)致故障的發(fā)生。
現(xiàn)有研究將單元故障作為隨機離散事件,重點關(guān)注從單元故障發(fā)生到系統(tǒng)失效的故障離散傳播過程。文獻[1]以核電廠供電系統(tǒng)為例,利用動態(tài)故障樹研究了考慮維修備件情況下系統(tǒng)的離散故障傳播過程,文獻[2]提出的動態(tài)貝葉斯網(wǎng)絡(luò)能夠較好地描述離散故障動態(tài)相關(guān)性;文獻[3]提出了基于Petri網(wǎng)的故障傳播模型,仿真了故障傳播及其診斷的過程;文獻[4]提出基于連續(xù)時間Markov鏈的連鎖故障相關(guān)性建模方法。
對于一般動態(tài)系統(tǒng),經(jīng)常出現(xiàn)性能降級等軟故障,這些故障對環(huán)境載荷、性能波動、參數(shù)退化、制造偏差等連續(xù)特性比較敏感。上述研究側(cè)重描述故障的離散傳播過程,解決系統(tǒng)規(guī)模大、相關(guān)性復(fù)雜等問題,大部分研究將故障完全視為隨機事件,卻忽略了故障對連續(xù)特性的作用,以及連續(xù)特性對故障的反作用。隨著研究的不斷深入,電網(wǎng)領(lǐng)域的部分學(xué)者開始研究連續(xù)特性變化(電力負(fù)載)對故障傳播的影響,如文獻[5- 6]提出了基于分支過程的負(fù)載相關(guān)連鎖故障建模方法。
結(jié)合故障傳播技術(shù)的發(fā)展趨勢,論文提出在動態(tài)系統(tǒng)故障傳播研究中,應(yīng)該從離散單一維度下的傳播逐漸擴展到離散和連續(xù)雙維度下的傳播,即故障的混雜傳播,充分考慮離散故障和連續(xù)特性之間的相互作用。這樣不僅能更真實地描述動態(tài)系統(tǒng)的故障傳播規(guī)律,還能分析得到對各種故障敏感的連續(xù)特性,為開展故障原因分析、可靠性優(yōu)化以及故障診斷提供支持。
論文第1節(jié)給出了動態(tài)系統(tǒng)故障混雜傳播的定義,并分析了混雜傳播過程的影響要素與特征;第2節(jié)探討了基于隨機混雜自動機的故障混雜傳播建模方法,利用模型離散層、連續(xù)層、交互層以及層次之間的動態(tài)相關(guān)性來描述故障的混雜傳播過程;第3節(jié)通過對某溫度控制系統(tǒng)的故障混雜傳播過程進行建模和仿真,驗證了動態(tài)系統(tǒng)故障傳播過程中的混雜特征,以及建模方法的可行性。
1.1 混雜傳播定義
動態(tài)系統(tǒng)故障傳播過程中,單元故障、物理參數(shù)退化和廣義應(yīng)力變化(如機械應(yīng)力、電應(yīng)力、熱應(yīng)力、液體力等)之間表現(xiàn)為兩兩相互作用,并且受到控制行為的有利調(diào)整以及環(huán)境因素的不利影響,使故障傳播變得更為復(fù)雜,表現(xiàn)出混雜性和不確定性。典型故障混雜傳播過程如圖1所示。
單元故障的發(fā)生(狀態(tài)2)使廣義應(yīng)力x1發(fā)生突變,隨著應(yīng)力在系統(tǒng)內(nèi)的非線性傳播,可能導(dǎo)致廣義應(yīng)力xn突破廣義強度,造成相關(guān)單元的過應(yīng)力故障或者性能降級故障(狀態(tài)3)。另一方面,廣義應(yīng)力x1的突變可能導(dǎo)致與其相關(guān)的物理參數(shù)p1的退化加劇,使p1超過閾值的時間變短,即p1相關(guān)的耗損性故障提前發(fā)生(狀態(tài)4)。在應(yīng)力、參數(shù)和故障三者的耦合影響下,驅(qū)動故障過程不斷發(fā)展直至系統(tǒng)故障(狀態(tài)5)。
將這種由離散和連續(xù)因素及其相互作用驅(qū)動的動態(tài)系統(tǒng)故障傳播,稱為動態(tài)系統(tǒng)故障的混雜傳播。傳統(tǒng)離散傳播將故障視為隨機離散時間,只能給出從狀態(tài)2到狀態(tài)4(或狀態(tài)3)的轉(zhuǎn)移概率;而本文綜合考慮連續(xù)傳播和混雜傳播后,即可給出狀態(tài)2到狀態(tài)4(或狀態(tài)3)的故障演化過程、故障相關(guān)應(yīng)力、敏感參數(shù)以及狀態(tài)4發(fā)生的確切時間等,為后續(xù)可靠性設(shè)計分析及故障診斷工作提供豐富的參考數(shù)據(jù)。
圖1 動態(tài)系統(tǒng)故障混雜傳播示例
1.2 混雜要素及特征分析
論文總結(jié)了影響動態(tài)系統(tǒng)故障傳播的5類主要離散和連續(xù)因素,即混雜要素:
(1)組成單元的故障事件(離散);
(2)廣義應(yīng)力的變化(短期連續(xù)波動、長期連續(xù)降級);
(3)物理參數(shù)由于磨損、腐蝕、老化等機理導(dǎo)致持續(xù)退化(連續(xù));
(4)環(huán)境影響(離散突發(fā)影響、連續(xù)漸變影響);
(5)自動及人工控制行為(離散重構(gòu)/維修、連續(xù)伺服調(diào)整)。
基于上述混雜因素及其耦合作用,可以描述動態(tài)系統(tǒng)故障傳播的混雜特征。這些特征在動態(tài)系統(tǒng)中普遍存在,但在傳統(tǒng)離散傳播中很難描述。典型特征如下:
(1)單元故障導(dǎo)致物理參數(shù)退化過程的加速,如余度系統(tǒng)某單元故障會導(dǎo)致其他單元分擔(dān)更多的應(yīng)力,進一步導(dǎo)致工作單元參數(shù)的加速磨損或者老化;另一方面物理參數(shù)的退化超過容差極限后,又會反過來導(dǎo)致單元故障,進一步導(dǎo)致系統(tǒng)故障。
(2)單元故障會導(dǎo)致某個廣義應(yīng)力的降級或者突變,系統(tǒng)性能降級超過閾值后,會反過來導(dǎo)致系統(tǒng)故障。
(3)物理參數(shù)的退化會導(dǎo)致廣義應(yīng)力的退化,而廣義應(yīng)力的變化又對物理參數(shù)退化起到加速或減速過程。
(4)環(huán)境變化直接導(dǎo)致廣義應(yīng)力的波動,并且會對物理參數(shù)的退化速率和單元故障率產(chǎn)生影響。
(5)控制手段會使系統(tǒng)向著遠(yuǎn)離故障的狀態(tài)演變,伺服控制可以容忍小擾動,容錯控制可以容忍故障發(fā)生。
混雜理論[7]是近年來在控制領(lǐng)域興起的一門新技術(shù),已被廣泛應(yīng)用于交通系統(tǒng)、機器人系統(tǒng)、電力系統(tǒng)等領(lǐng)域[8-9],通過混雜模型可以將“切換”等離散控制信號與連續(xù)動態(tài)信號統(tǒng)一到同一建??蚣芟拢瑢崿F(xiàn)仿真過程中離散模型和連續(xù)模型的動態(tài)交互。
但現(xiàn)有混雜模型注重描述系統(tǒng)正常工作時的確定性演化過程,對故障狀態(tài)及隨機特征描述不足。如文獻[10]中雖給出了隨機混雜自動機(stochastic hybrid automata,SHA)的定義,但該模型只考慮了離散狀態(tài)的切換隨機性,在離散狀態(tài)內(nèi)部的演化由確定性微分方程驅(qū)動。事實上,在離散狀態(tài)內(nèi)部,系統(tǒng)的連續(xù)演化規(guī)律也是受隨機性影響的,如環(huán)境擾動、輸入噪聲以及制造誤差。
為了能將控制領(lǐng)域的SHA運用到故障混雜傳播建模中,論文提出了一種新的SHA模型,更全面地考慮了各種隨機因素及其對故障傳播的影響。
定義SHA為一個9元組:
式中,Q為離散狀態(tài)集Q={q1,q2,…,qm};X為連續(xù)狀態(tài)空間X=[x1,x2,…,xn];E為有限的事件集,包括確定性事件和隨機性事件;W=(u,w)為外部輸入集合,u為控制信號,w為擾動信號;Act是指定在每個離散狀態(tài)q∈Q下的映射,反映系統(tǒng)在該離散狀態(tài)下的運動規(guī)律(如應(yīng)力變化規(guī)律、參數(shù)退化規(guī)律等),可由隨機微分方程和數(shù)值仿真模型表示。由于各離散狀態(tài)表征的故障或工作模式不同,因此各微分方程或仿真模型的結(jié)構(gòu)或參數(shù)也不同。Λ為遷移強度函數(shù),Λ→λij(x,t)表示在t時刻,x點處發(fā)生從狀態(tài)qi到qj遷移的瞬時強度;C為遷移關(guān)系集合C=(qi,eij,Γij(x,t),qj),qi,qj分別為起始狀態(tài)和目標(biāo)狀態(tài),eij為對應(yīng)的事件,Γij(x,t)為跳變條件。對于確定性跳變
式中,G(x)表示邊界判據(jù),對應(yīng)廣義應(yīng)力x或者物理參數(shù)A的閾值。對于隨機性跳變
2.1 隨機混雜自動機定義
混雜故障傳播雖然能更真實地描述系統(tǒng)故障過程,但同時會為建模帶來難度,需要將離散的故障狀態(tài)轉(zhuǎn)移模型與連續(xù)的應(yīng)力模型和參數(shù)退化模型有機整合起來,不僅需要實現(xiàn)各類模型之間的動態(tài)作用,還要考慮各類混雜要素隨機性的影響,這給混雜傳播建模帶來了很大的難度。論文基于混雜理論初步探討了動態(tài)系統(tǒng)故障混雜傳播的建模方法。其中,函數(shù)g表示執(zhí)行狀態(tài)切換所需滿足的瞬時遷移強度λij(x(t))的積累限度。P=[]為概率分布矩陣,p表示從狀態(tài)qi到ql的概率;Inv=(q0,x0,P0)對應(yīng)初始條件。
2.2 基于SHA的故障混雜傳播模型
基于上述SHA模型,系統(tǒng)故障混雜傳播模型可表述為離散層、連續(xù)層和交互層,結(jié)構(gòu)如圖2所示。
圖2 系統(tǒng)故障混雜傳播模型結(jié)構(gòu)
離散層由自動機模型描述,主要表征長周期尺度下系統(tǒng)正常、降級、維修、故障運行狀態(tài)q之間的邏輯關(guān)系,與傳統(tǒng)離散傳播模型類似。連續(xù)層由隨機微分方程描述,主要表征系統(tǒng)在各個離散狀態(tài)下系統(tǒng)的連續(xù)演化規(guī)律,包括廣義應(yīng)力模型、退化模型以及故障率模型。廣義應(yīng)力模型描述在離散狀態(tài)qi下,應(yīng)力x在物理參數(shù)A和B、控制輸入u以及環(huán)境擾動(Wiener過程d w(t))共同驅(qū)動下的動態(tài)規(guī)律;退化模型描述物理參數(shù)A和B在應(yīng)力x和擾動d w(t)下的變化過程;故障率模型給出由離散狀態(tài)qi跳變到qj的概率λij,考慮了應(yīng)力x對故障率的影響。3類模型相互作用,驅(qū)動當(dāng)前離散狀態(tài)qi下的系統(tǒng)演化。
交互層借助SHA的跳變機制實現(xiàn)離散模型與連續(xù)模型間的動態(tài)交互。一方面根據(jù)離散層的不同故障狀態(tài)qi執(zhí)行連續(xù)層的模型重構(gòu);另一方面實時檢測連續(xù)層的應(yīng)力、物理參數(shù)及故障率,動態(tài)決定離散層跳變到故障狀態(tài)qj的準(zhǔn)確時間,并能判斷出故障原因是廣義應(yīng)力(如屈服應(yīng)力超過強度、性能降級)還是物理參數(shù)(如磨損或裂紋超過閾值),或者是偶然故障。
某溫度控制系統(tǒng)包括兩個互為冗余的控制回路(PI控制和開關(guān)控制)和一個烤箱箱體。PI控制模式是烤箱系統(tǒng)的主要運行模式,保證烤箱溫度快速穩(wěn)定到指定溫度Tref=190°C。當(dāng)PI控制器故障時,切換為開關(guān)控制模式,通過在加熱和自然降溫之間實現(xiàn)切換,使溫度維持在[170℃,210℃]之間,并啟動PI控制器的維修過程。
案例做出以下假設(shè):
(1)僅考慮烤箱、PI控制器和開關(guān)控制器的故障。當(dāng)溫度超過270℃時,烤箱故障;將PI控制器、開關(guān)控制器視為隨機事件,故障率隨溫度變化。λ(x,t)=h(x)λ0(t),溫度系數(shù)[11]為
(2)考慮烤箱物理參數(shù)的退化。加熱系數(shù)τ服從均值退化的正態(tài)分布τ~N(0.5×e-0.0001t,0.1)。
(3)考慮控制信號噪聲影響。設(shè)信號噪聲為白噪聲,信噪比為65 dB。
(4)考慮控制器維修過程修復(fù)如舊,衰減系數(shù)α=0.9,累計衰減超過0.5時更換控制器,烤箱不可修。
(5)考慮伺服控制對溫度的調(diào)節(jié)作用。
(6)考慮故障檢測,當(dāng)溫度超出[140℃,240℃],傳感器報故。
建模過程如下:
步驟1 根據(jù)系統(tǒng)的運行原理,定義系統(tǒng)的主要離散狀態(tài),如表1所示。
表1 系統(tǒng)離散狀態(tài)
步驟2 根據(jù)系統(tǒng)工作原理、故障狀態(tài)、退化特征以及噪聲分布等,得到系統(tǒng)各個離散狀態(tài)下的連續(xù)狀態(tài)方程。
步驟3 確定各個離散狀態(tài)之間的跳變條件,包括隨機跳變和確定跳變,建立起連續(xù)狀態(tài)與離散狀態(tài)間的交互關(guān)系。
步驟4 根據(jù)上述分析結(jié)果,利用Simulink/Stateflow建立溫控系統(tǒng)的故障傳播模型。在Stateflow中建立離散層和交互層模型,刻畫溫控系統(tǒng)的狀態(tài)轉(zhuǎn)移關(guān)系,在Simulink中建立連續(xù)層模型,刻畫系統(tǒng)在不同運行狀態(tài)下的動態(tài)性能。
步驟5 設(shè)置仿真時間、仿真步長等參數(shù),進行故障傳播過程仿真,輸出仿真結(jié)果。
圖3為某一條故障混雜傳播路徑。其中,左側(cè)坐標(biāo)是系統(tǒng)溫度,上方曲線為溫度曲線,兩條虛線之間為溫度控制系統(tǒng)正常工作的溫度范圍[170℃,210℃],右側(cè)坐標(biāo)是離散狀態(tài),下方曲線為狀態(tài)變化曲線。
圖3 典型混雜故障傳播過程
從圖3可以看出,系統(tǒng)故障傳播是在離散狀態(tài)和連續(xù)溫度相互作用下動態(tài)演化的。如果只將故障視為隨機事件,而不考慮連續(xù)溫度的實時影響,將難以準(zhǔn)確給出故障狀態(tài)的變化過程。
為了準(zhǔn)確觀察故障混雜傳播過程,截?。?20 h,228 h]系統(tǒng)溫度及相應(yīng)運行狀態(tài)的變化,如圖4中上、下兩條曲線所示。
圖4 [220 h,228 h]的離散狀態(tài)遷移和溫度變化過程
從圖4可以看出,離散故障、維修事件與連續(xù)溫度之間的相互作用:
(1)220 h時系統(tǒng)處于正常的PI運行狀態(tài)(狀態(tài)1),溫度穩(wěn)定在190℃。
(2)之后PI控制器故障(狀態(tài)2),溫度急劇升高至控制器故障的檢測溫度240℃,檢出故障(A點)。
(3)切換至開關(guān)控制器(狀態(tài)3),控制溫度開始下降。
(4)在開關(guān)控制器作用下,溫度控制在[170℃,210℃]范圍內(nèi)(狀態(tài)3、狀態(tài)4)。
(5)在B點開關(guān)控制器也發(fā)生了故障(狀態(tài)5),使溫度迅速上升至控制器故障的檢測溫度240℃。
(6)C點開關(guān)控制器故障被檢測單元檢出并且PI控制器已修好,系統(tǒng)切換至PI控制狀態(tài)(狀態(tài)1)。
論文提出了動態(tài)系統(tǒng)的故障傳播具有混雜特征,并初步討論了混雜傳播的建模方法?;赟HA的動態(tài)系統(tǒng)故障傳播模型能夠有效地描述不確定性因素作用下,離散單元故障和連續(xù)過程參數(shù)之間的交互作用,能更真實地描述系統(tǒng)故障傳播過程,為動態(tài)系統(tǒng)故障規(guī)律認(rèn)知、故障原因分析和故障診斷提供參考數(shù)據(jù)。
但在混雜傳播計算過程中,離散模型采樣跨度大,適合采用事件驅(qū)動機制,連續(xù)動態(tài)模型采樣跨度小,適合采用時間驅(qū)動機制。離散和連續(xù)2種采樣機制步長相差多個數(shù)量級,還需研究針對性的步長控制算法,通過變步長的方式來提高計算效率。
[1]Rao K D,Gopika V,Rao V S,et al.Dynamic fault tree analysis using Monte Carlo simulation in probabilistic safety assessment[J].Reliability Engineering&System Safety,2009,94(4):872- 883.
[2]Spackova O,Straub D.Dynamic Bayesian network for probabilistic modeling of tunnel excavation processes[J].Computer Aided Civil and Infrastructure Engineering,2013,28(1):1- 21.
[3]Georgilakis P S.Petri net based transformer fault diagnosis[C]∥Proc.of the IEEE International Symposium on Circuits and Systems,2004:980- 983.
[4]Srinivasan M I,Marvin K N,Alexandros V G.A Markovian dependability model with cascading failures[J].IEEE Trans.on Computers,2009,58(9):1238- 1249.
[5]Dobson I,Carreras B A,Newman D E.A loading-dependent model of probabilistic cascading failure[J].Probability in the Engineering and Informational Sciences,2005,19(1):15 -32.
[6]Kim J,Dobson I.Approximating a loading-dependent cascading failure model with a branching process[J].IEEE Trans.on Reliability,2010,59(4):691- 699.
[7]Van Der Schaft A T,Schumacher J M.An introduction to hybrid dynamical systems[M].London:Springer-Verlag,2000.
[8]Seah C E,Hwang I.Stochastic linear hybrid systems:modeling,estimation,and application in air traffic control[J].IEEE Trans.on Control Systems Technology,2009,17(3):563- 575.
[9]Yin A,Zhao H,Zhang B.Control of hybrid electric bus based on hybrid system theory[C]∥Proc.of the International Conference on Electric Information and Control Engineering,2011:1253- 1256.
[10]Castaneda G P,Aubry J F,Brinzei N.Stochastic hybrid automata model for dynamic reliability assessment[J].Proceedings of the Institution of Mechanical Engineers,Part O:Journal of Risk and Reliability,2011,225(1):28- 41.
[11]Marseguerra M,Zio E.The cell-to-boundary method in Monte Carlo-based dynamic PSA[J].Reliability Engineering&System Safety,1995,48(3):199- 204.
Fault hybrid propagation and modeling method for dynamic system
GUO Jian-bin1,2,DU Shao-hua1,WANG Xin3,ZENG Sheng-kui1,2
(1.School of Reliability and Systems Engineering,Beihang University,Beijing 100191,China;2.Science and Technology on Reliability and Environment Engineering Laboratory,Beijing 100191,China;3.Beijing Institute of Environment Features,Beijing 100854,China)
The failure propagation in dynamic systems is driven by discrete component fault events,continuous processes as well as their interactions.This hybrid feature of the fault propagation brings about the difficulty in fault cognition and modeling.Existing researches regard faults just as discrete events.Thus they only focus on analyzing how the discrete system failures are caused by random component failures.However,the continuous processes in fault propagation have always been ignored for the sake of simplifications in engineering,which leads to the inaccuracy of the description of dynamic system failures.This paper defines the hybrid failure propagation within two dimensions for the dynamic system,and analyzes its hybrid factors and propagation features.For accurate description of these hybrid features,a hybrid failure propagation modeling method is proposed based on the hybrid theory which is used to model the interaction of discrete events and continuous processes.The proposed method is applied to a temperature control system.The simulation results show the hybrid feature of the failure propagation,as well as the feasibility of the presented modeling method.
fault propagation;hybrid;reliability;dynamic system;stochastic hybrid automata(SHA)
TB 114.3
A
10.3969/j.issn.1001-506X.2015.01.36
郭健彬(1979-),男,講師,博士,主要研究方向為系統(tǒng)可靠性、故障機理、可靠性優(yōu)化。
E-mail:guojianbin@buaa.edu.cn
杜紹華(1987-),男,碩士研究生,主要研究方向為系統(tǒng)可靠性、靈敏度分析。
E-mail:h2oh2o2@126.com
王 鑫(1985-),男,工程師,碩士,主要研究方向為系統(tǒng)可靠性、混雜系統(tǒng)。
E-mail:w_xin_ok@126.com
曾聲奎(1968-),男,研究員,博士,主要研究方向為系統(tǒng)可靠性、故障機理、可靠性與性能一體化設(shè)計。
E-mail:zengshengkui@buaa.edu.cn
1001-506X(2015)01-0224-05
網(wǎng)址:www.sys-ele.com
2013- 12- 27;
2014- 06- 20;網(wǎng)絡(luò)優(yōu)先出版日期:2014- 08- 20。
網(wǎng)絡(luò)優(yōu)先出版地址:http://w ww.cnki.net/kcms/detail/11.2422.TN.20140820.1745.007.html
國家自然科學(xué)基金(61304218)資助課題