基于非高斯仿真的風(fēng)壓系數(shù)極值計(jì)算方法

李壽科，李壽英，陳政清，孫洪鑫

（1.湖南科技大學(xué)土木工程學(xué)院，湖南 湘潭 411201；2.北京交通大學(xué)結(jié)構(gòu)風(fēng)工程與城市風(fēng)環(huán)境北京市重點(diǎn)實(shí)驗(yàn)室，北京：100044；3.湖南大學(xué)風(fēng)工程試驗(yàn)研究中心，長沙 410082）

基于非高斯仿真的風(fēng)壓系數(shù)極值計(jì)算方法

李壽科1，2，李壽英3，陳政清3，孫洪鑫1

（1.湖南科技大學(xué)土木工程學(xué)院，湖南 湘潭 411201；2.北京交通大學(xué)結(jié)構(gòu)風(fēng)工程與城市風(fēng)環(huán)境北京市重點(diǎn)實(shí)驗(yàn)室，北京：100044；3.湖南大學(xué)風(fēng)工程試驗(yàn)研究中心，長沙 410082）

以多變量相關(guān)非高斯過程仿真方法為基礎(chǔ)，發(fā)展了一種基于單次采樣的多變量非高斯仿真極值計(jì)算方法。首先介紹開孔屋蓋的風(fēng)洞試驗(yàn)概況和多變量相關(guān)非高斯過程仿真的基本理論，對屋蓋上一組測點(diǎn)風(fēng)壓進(jìn)行了非高斯仿真，結(jié)果表明基于譜修正的多變量相關(guān)非高斯過程仿真方法得到的時程在功率譜密度，相干函數(shù)，高階矩三方面與目標(biāo)值接近，仿真效果較好，然后采用經(jīng)典極值理論對多次仿真的非高斯時程進(jìn)行極值計(jì)算，將該方法得到的峰值因子與以往常用方法的結(jié)果進(jìn)行比較，結(jié)果表明：Davenport峰值因子法高估氣流分離區(qū)左偏風(fēng)壓的正峰值因子60%，低估負(fù)峰值因子43%；Sadek－Simiu峰值因子法低估了高峰度風(fēng)壓的峰值因子50%；而基于單次樣本進(jìn)行仿真的非高斯仿真峰值因子法，其估計(jì)的開孔屋蓋的峰值因子最為準(zhǔn)確，與觀察峰值因子總體上最為接近。

非高斯仿真；風(fēng)壓極值；峰值因子；風(fēng)洞試驗(yàn)

根據(jù)風(fēng)災(zāi)調(diào)查，建筑圍護(hù)結(jié)構(gòu)的破壞在風(fēng)災(zāi)破壞損失中占很大的比例。建筑圍護(hù)結(jié)構(gòu)的設(shè)計(jì)由其表面的極值風(fēng)荷載決定，通常得到極值風(fēng)荷載的方法主要有兩類。第一類為基于多次獨(dú)立采樣的經(jīng)典極值理論方法，該方法對每次采樣的極值進(jìn)行分布擬合，從而確定出具有一定保證率的極值風(fēng)壓。該方法對非高斯過程和高斯過程具有同樣的適應(yīng)性，被認(rèn)為是較準(zhǔn)確的一種極值統(tǒng)計(jì)方法，但此種方法需對實(shí)際信號進(jìn)行多次重復(fù)獨(dú)立的采樣，將耗費(fèi)較大的人力和物力，在實(shí)際過程中卻往往較少采用。第二類為基于單次采樣的零值穿越理論的方法。對于服從高斯分布的隨機(jī)信號，Davenport［1］在零值穿越理論的基礎(chǔ)上，獲得了服從窄帶高斯分布隨機(jī)信號的峰值因子，此種方法使用簡單方便，在早期工程實(shí)際中被廣泛采用。但很多時候風(fēng)壓并不服從通常的高斯分布，風(fēng)壓的非高斯特征對風(fēng)壓極值的估計(jì)具有較大的影響，Holmes［2］及Gioffre［3－4］對此做了相關(guān)研究，結(jié)果表明來流分離區(qū)風(fēng)壓的非高斯特性對風(fēng)壓極值的影響很大，按照高斯假定估計(jì)的峰值因子明顯偏小。非高斯風(fēng)壓具有比高斯風(fēng)壓更大的破壞性，Holmes［2］發(fā)現(xiàn)非高斯分布風(fēng)壓導(dǎo)致的結(jié)構(gòu)破壞要比高斯分布的風(fēng)壓大15%～30%，所以對于屋蓋表面風(fēng)壓極值的估計(jì)需考慮其概率分布特性，兼顧其非高斯特性。對于非高斯過程的極值風(fēng)壓，Sadek［5］以Gamma分布和高斯分布作為母體分布，提出了一種非高斯轉(zhuǎn)換過程的極值風(fēng)壓計(jì)算方法，而Kareem等［6］在Davenport方法的基礎(chǔ)上，通過將非高斯過程展開為標(biāo)準(zhǔn)高斯變量的Hermite級數(shù)，將僅適用于窄帶高斯過程的基于零值穿越理論的峰值因子法擴(kuò)展到了非高斯過程。全涌等［7］提出了一種對單次標(biāo)準(zhǔn)長度的非高斯風(fēng)壓時程數(shù)據(jù)進(jìn)行分段，然后通過子段的極值分布規(guī)律估算出母段的期望極值的改進(jìn)經(jīng)典極值法，在此方法中選擇合適的采樣長度較為重要。

Gurley［8］在進(jìn)行非高斯風(fēng)壓場仿真研究時采用了Hermite級數(shù)來描述非高斯風(fēng)壓的概率分布，結(jié)果表明Hermite級數(shù)可以較好的體現(xiàn)隨機(jī)過程的高階矩特性從而達(dá)到可描述非高斯過程概率密度函數(shù)的目的。本文針對多次獨(dú)立采樣的經(jīng)典極值方法需耗費(fèi)大量人力物力的缺點(diǎn)，研究多變量的相關(guān)非高斯隨機(jī)過程仿真方法，提出一種基于單次采樣結(jié)合多變量相關(guān)非高斯隨機(jī)過程仿真的峰值因子估計(jì)方法。最后，基于一個屋蓋開孔建筑的剛性模型測壓試驗(yàn)數(shù)據(jù)，通過與以往幾種方法的比較，顯示本方法的進(jìn)步之處。

1 剛性模型測壓試驗(yàn)概況

試驗(yàn)在湖南大學(xué)HD－2風(fēng)洞的高速試驗(yàn)段進(jìn)行。試驗(yàn)?zāi)Ｐ筒捎糜袡C(jī)玻璃制作，實(shí)際結(jié)構(gòu)在主要尺寸上與TTU建筑（13.72×9.14×3.96 m）保持幾何相似，在屋蓋的中心進(jìn)行15%的開孔，模型的幾何縮尺比為1∶50，試驗(yàn)照片如圖1所示。模型立墻表面布置56個測點(diǎn)，屋蓋上下表面布置120個測點(diǎn)，上下表面測點(diǎn)位置對應(yīng)，詳細(xì)測點(diǎn)布置見圖2。試驗(yàn)的風(fēng)向角定義見圖2，風(fēng)向角間隔5°，共72個測試風(fēng)向角。采樣時長33 s，采樣頻率330 Hz，共采集10 000個數(shù)據(jù)點(diǎn)，所有采樣后數(shù)據(jù)采用管道頻響函數(shù)進(jìn)行畸變修正，且對典型風(fēng)向角0°、45°、90°、270°、315°進(jìn)行10次獨(dú)立重復(fù)采樣。試驗(yàn)風(fēng)速11.0 m/s，參考高度為8 cm，相當(dāng)于實(shí)際高度4 m。試驗(yàn)?zāi)M了B類地貌風(fēng)場［9］，風(fēng)場比例為1∶50，平均風(fēng)剖面指數(shù)為0.15，湍流度剖面也與實(shí)際大氣中的情況基本一致，10 m高度處為0.20左右。

圖1 試驗(yàn)?zāi)Ｐ驼掌現(xiàn)ig.1 Photo of testmodel

圖2 試驗(yàn)?zāi)Ｐ蜏y點(diǎn)布置圖Fig.2 Tap location of testmodel

2 數(shù)據(jù)處理方法

風(fēng)壓系數(shù)是結(jié)構(gòu)風(fēng)壓的無量綱表現(xiàn)形式，測點(diǎn)i的風(fēng)壓系數(shù)CPi（t）定義如下：

其中：Pi（t）為風(fēng)洞試驗(yàn)中壓力掃描閥測得的風(fēng)壓時程；P0為風(fēng)洞試驗(yàn)段處的靜壓，采用皮托管測得；ρ為空氣密度，取ρ＝1.225 kg/m3；uh為屋蓋最高點(diǎn)處的平均風(fēng)速。CPi（t）的平均值CPi＿mean為平均風(fēng)壓系數(shù)，CPi（t）的脈動值CPi＿rms為脈動風(fēng)壓系數(shù)，極大值風(fēng)壓系數(shù)為)CPi和極小值風(fēng)壓系數(shù)為(CPi。風(fēng)壓系數(shù)的正負(fù)峰值因子定義為

上表面減下表面風(fēng)壓系數(shù)為凈風(fēng)壓系數(shù)，根據(jù)式（2）可得相應(yīng)的峰值因子。

3 基于多變量相關(guān)非高斯過程仿真的極值計(jì)算理論

3.1 基于譜修正的多變量相關(guān)非高斯過程的仿真

仿真多變量非高斯相關(guān)過程，需要保證仿真結(jié)果的前四階矩、功率譜密度和相干函數(shù)與測量目標(biāo)值一致，首先需要得到多變量的高斯相關(guān)過程，繼而進(jìn)行非高斯過程轉(zhuǎn)換，其具體步驟如下。

（1）基于譜分解方法的多變量高斯仿真

已知n個相關(guān)過程的自譜和它們之間的互譜，其互譜密度矩陣可以用切比雪夫分解表示為：

式（3）中Gii（ω）為第i個變量頻率點(diǎn)ω處的自譜，Gij（ω）為第i個變量和第j個變量頻率點(diǎn)ω處的互譜，其中：

式（4）中Rij（ω）為互相干函數(shù)，參考Gurley的方法［8］，產(chǎn)生兩個獨(dú)立的0均值和Δω根方差白噪聲序列η和ζ組成復(fù)數(shù)序列ξ＝η＋kζ，則這n個變量的傅里葉序列可以表示為：

對式（5）進(jìn)行傅里葉逆變換則可得到n個相關(guān)的高斯時間序列Yi（t）。

（2）單變量基于譜修正方法的非高斯仿真

在得到n個相關(guān)的高斯過程后，需將n個高斯過程利用Hermite級數(shù)轉(zhuǎn)換方法轉(zhuǎn)換到非高斯過程，轉(zhuǎn)換后的變量的前四階矩和功率譜密度需要保證與目標(biāo)變量值一致，如不滿足則需要進(jìn)一步迭代修正，具體步驟如下。

（?。┫蚯胺歉咚罐D(zhuǎn)換

應(yīng)用向前Hermite轉(zhuǎn)換方法將高斯過程Y（t）產(chǎn)生非高斯過程Xng（t）：

轉(zhuǎn)換后的非高斯變量滿足目標(biāo)值的前四階矩，但其功率譜密度將被扭曲。

（ⅱ）功率譜修正

保持非高斯過程的相位不變，利用目標(biāo)譜值進(jìn)行修正形成新的非高斯序列：

式（7）中GT為目標(biāo)功率譜值。修正后序列的偏度和峰度將可能不再滿足目標(biāo)矩，設(shè)定偏度和峰度的誤差判別標(biāo)準(zhǔn)為：

設(shè)定一誤差限，假如誤差超過其誤差限，則進(jìn)行下一步驟（ⅲ），如果在誤差限內(nèi)則退出，則得到了多條相關(guān)性可能不滿足的非高斯過程，繼續(xù)步驟（ⅲ）。

（ⅲ）向后高斯轉(zhuǎn)換

由于沒有達(dá)到目標(biāo)非高斯過程，需再次進(jìn)行仿真迭代，則需將步驟（ⅱ）得到的非高斯過程采用式（6）的反函數(shù)進(jìn)行Hermite高斯過程轉(zhuǎn)換，然后再返回到步驟（?。?，進(jìn)行迭代求解，直至達(dá)到誤差限。將式（6）求反函數(shù)，可得u（x）為：

（3）相干函數(shù)修正

采用單變量非高斯仿真方法可以得到多條非高斯過程，但各個變量之間的相關(guān)性將得不到滿足，此時可對其相干函數(shù)也進(jìn)行迭代修正，其設(shè)定誤差判別標(biāo)準(zhǔn)為：

式（11）中：RXng（f）為迭代前非高斯過程的相干函數(shù)；RT

ijij（f）為目標(biāo)相干函數(shù)；tol為相干函數(shù)誤差限。如果結(jié)果在誤差限內(nèi)則退出，如果結(jié)果不在誤差范圍內(nèi)，更新相干函數(shù)，返回到開始步驟（i）繼續(xù)仿真，其更新相干函數(shù)表達(dá)如下：

3.2 非高斯仿真峰值因子計(jì)算

為利用多變量相關(guān)過程非高斯仿真方法獲得峰值因子，首先可利用多變量非高斯仿真方法對標(biāo)準(zhǔn)化測點(diǎn)時程進(jìn)行多次仿真，可得到多條時程，對每條時程求得觀察極值，對所有觀察極值進(jìn)行經(jīng)典極值分布（Gumbel分布）擬合，則可得到一定保證率下的峰值因子，在后文中稱為非高斯仿真峰值因子，進(jìn)行平均即可得到與Davenport峰值因子法［1］、Sadek-Simiu方法［5］對應(yīng)的峰值因子。如果仿真時程不進(jìn)行標(biāo)準(zhǔn)化，則可直接得到相應(yīng)的風(fēng)壓系數(shù)極大和極小值。

4 非高斯仿真極值算例

4.1 多變量相關(guān)非高斯過程仿真實(shí)例

為驗(yàn)證基于譜修正多變量相關(guān)過程的非高斯仿真的有效性，基于本文試驗(yàn)結(jié)果，對開孔屋蓋0°風(fēng)向角時的76、77、78測點(diǎn)的風(fēng)壓時程進(jìn)行仿真，仿真測點(diǎn)布置如圖2所示，設(shè)定誤差上限為10%，當(dāng)然，測點(diǎn)數(shù)越多，誤差上限設(shè)定越小，需要迭代計(jì)算的時間越長。圖3給出了76、77、78測點(diǎn)的風(fēng)壓系數(shù)功率譜和相干函數(shù)的仿真和試驗(yàn)結(jié)果比較。由圖3可以看出，測點(diǎn)風(fēng)壓之間的相干函數(shù)仿真結(jié)果和試驗(yàn)結(jié)果在低頻處非常接近，在高頻處稍有偏差，測點(diǎn)風(fēng)壓的功率譜密度函數(shù)的仿真結(jié)果和試驗(yàn)結(jié)果比較一致。表1給出了仿真時程的高階矩與目標(biāo)矩之間的對比，從表1可以看出，仿真風(fēng)壓的偏度和峰值均與目標(biāo)試驗(yàn)值較為接近，其誤差在設(shè)定的誤差限10%以內(nèi)。由此可以看出，多變量的相關(guān)非高斯過程仿真方法可以較準(zhǔn)確的仿真出多條相關(guān)的非高斯風(fēng)壓時程，其仿真結(jié)果的功率譜密度、相干函數(shù)以及高階矩均在目標(biāo)偏差范圍內(nèi)。

圖3 三變量相關(guān)過程非高斯仿真結(jié)果和測量結(jié)果頻域統(tǒng)計(jì)比較Fig.3 Comparison between simulation and target in frequency domain

表1 目標(biāo)時程和仿真時程的2、3、4階矩比較Tab.1 Com parsion of deviation，skewness，ku rtosis between simulation and target

4.2 非高斯仿真峰值因子法與以往峰值因子法的比較

圖4給出了開孔屋蓋在0°、45°、90°風(fēng)向角時的橫向中軸線測點(diǎn)凈風(fēng)壓系數(shù)的峰度和偏度分布規(guī)律。由圖4可以看出，對于0°風(fēng)向角，迎風(fēng)區(qū)測點(diǎn)風(fēng)壓左偏，尾流區(qū)測點(diǎn)風(fēng)壓右偏，邊緣分離處測點(diǎn)風(fēng)壓出現(xiàn)大負(fù)偏度（－1.23）和高峰度（6.21），表現(xiàn)出明顯非高斯特性，背風(fēng)尾流區(qū)測點(diǎn)表現(xiàn)出中等非高斯（偏度約為0.5，峰度值大于4.2）；45°風(fēng)向角時，測點(diǎn)風(fēng)壓偏度較0°風(fēng)向角小，其相應(yīng)的峰度值也減小；90°風(fēng)向角時，屋蓋橫向中軸線測點(diǎn)風(fēng)壓偏度均為負(fù)，在三個典型風(fēng)向角中達(dá)到最小，此時測點(diǎn)處風(fēng)壓分離程度最小，但峰度值明顯偏離3，在三個典型風(fēng)向角中（0°、45°、90°）為最大。由此可以看出，屋蓋中軸線測點(diǎn)明顯偏離標(biāo)準(zhǔn)正態(tài)高斯分布（偏度為0峰度為3），風(fēng)壓的偏度主要由氣流的分離程度決定，氣流分離越明顯，其風(fēng)壓偏度越大。

圖4 屋蓋橫向中軸線測點(diǎn)風(fēng)壓系數(shù)偏度和峰度分布規(guī)律Fig.4 Skewness and kurtosis ofwind pressures on opening roof

圖5給出了開孔屋蓋在0°、45°、90°風(fēng)向角時的橫向中軸線測點(diǎn)凈風(fēng)壓在基于多變量相關(guān)過程非高斯仿真峰值因子法、Sadek-Simiu的轉(zhuǎn)換過程法、Davenport峰值因子法和觀察極值法下的正負(fù)峰值因子分布規(guī)律，各種方法對應(yīng)的結(jié)果為非高斯仿真峰值因子、Sadek-Simiu峰值因子、Davenport峰值因子和觀察峰值因子，其中非高斯仿真峰值因子為通過多變量相關(guān)非高斯過程仿真方法進(jìn)行非高斯仿真16次取極值的平均值獲得。在0°風(fēng)向角時，迎風(fēng)區(qū)分離測點(diǎn)風(fēng)壓左偏，即風(fēng)壓概率密度函數(shù)左邊表現(xiàn)為長尾部，由圖5（a）可以看出，Davenport峰值因子法低估了負(fù)峰值因子達(dá)43%，高估了正峰值因子達(dá)60%，非高斯仿真峰值因子稍大于觀察峰值因子，而Sadek-Simu峰值因子則與觀察峰值因子較為接近；估計(jì)負(fù)峰值因子時，Davenport峰值因子法明顯低估了迎風(fēng)區(qū)屋蓋的負(fù)峰值因子，Sadek-Simu峰值因子小于觀察峰值因子，非高斯仿真峰值因子與觀察峰值因子較為接近。在45°風(fēng)向角時，視風(fēng)壓為高斯分布的Davenport峰值因子與觀察峰值因子的誤差較大，Sadek-Simu峰值因子則小于觀察峰值因子，非高斯仿真峰值因子與觀察峰值因子較為接近；在90°風(fēng)向角時，Davenport峰值因子明顯低估負(fù)偏度高峰度的屋蓋風(fēng)壓正負(fù)峰值因子，其偏差高達(dá)50%，而Sadek-Simu峰值因子則明顯小于觀察峰值因子，非高斯仿真峰值因子與觀察峰值因子較為接近。

綜合不同方法對開孔屋蓋的峰值因子計(jì)算結(jié)果比較可以得出以下結(jié)論：①Sadek-Simiu峰值因子法對0°風(fēng)向角的正峰值因子具有較好的估計(jì)，但在一定程度上低估了開孔屋蓋多數(shù)風(fēng)向角下的峰值因子，且對于高峰度的風(fēng)壓峰值有明顯的低估；②Davenport峰值因子法由于忽略風(fēng)壓的非高斯特性，對開孔屋蓋的左偏風(fēng)壓，高估正峰值因子，低估負(fù)峰值因子；③基于觀察樣本進(jìn)行仿真的非高斯仿真峰值因子法，其估計(jì)的開孔屋蓋的峰值因子最為準(zhǔn)確，與觀察峰值因子總體上最為接近。

圖5 不同計(jì)算方法的峰值因子比較Fig.5 Comparison of peak factors by differentmethods

5 結(jié) 論

（1）基于譜修正的多變量相關(guān)非高斯過程仿真方法仿真出的多條非高斯時程的功率譜密度，相干函數(shù)，根方差，偏度和峰度與目標(biāo)值吻合得很好，誤差在目標(biāo)范圍內(nèi)；

（2）開孔屋蓋測點(diǎn)風(fēng)壓明顯偏離標(biāo)準(zhǔn)正態(tài)高斯分布，風(fēng)壓的偏度主要由氣流的分離程度決定，氣流分離越劇烈，其風(fēng)壓偏度越大；

（3）Davenport峰值因子法由于忽略風(fēng)壓的非高斯特性，對開孔屋蓋氣流分離區(qū)的左偏風(fēng)壓，高估正峰值因子60%，低估負(fù)峰值因子43%；

（4）Sadek-Simiu峰值因子法對于高峰度的風(fēng)壓峰值有明顯的低估，偏差高達(dá)50%；基于觀察樣本進(jìn)行仿真的非高斯仿真峰值因子法，其估計(jì)的開孔屋蓋的峰值因子最為準(zhǔn)確，與觀察峰值因子總體上最為接近。

［1］Davenport A G.Note on the distribution of the largest value of a random function with application to gust loading［J］.In：Proc.ICE，1964（28）：187－195.

［2］Holmes JD.Wind action on glass and Brown's integral［J］.Engineering Structures，1985，7（4）：226－230.

［3］GioffrèM，Gusella V，Grigoriu M.Non-Gaussian wind pressure on prismatic buildings.I：Stochastic field［J］.Journal of Structural Engineering，2001，127（9）：981－989.

［4］GioffrèM，Gusella V，Grigoriu M.Non-Gaussian wind pressure on prismatic buildings.II：Numerical simulation［J］.Journal of Structural Engineering，2001，127（9）：990－995.

［5］Sadek F E S.Peak non-gaussian wind effects for databaseassisted low-rise building design［J］.Journal of Engineering Mechanics-Asce，2002，5（128）：530－539.

［6］Kareem A，Zhao J，TognarelliMA.Surge response statistics of tension leg platforms under wind and wave loads：a statistical quadratization approach［J］.Probabilistic Engineering Mechanics，1995，10（4）：225－240.

［7］全涌，顧明，陳斌.非高斯風(fēng)壓的極值計(jì)算方法［J］.力學(xué)學(xué)報，2010（3）：560－566.

QUAN Yong，GUMing，CHEN Bin.Study on the extreme value estimationmethod of non-gaussian wind pressure［J］.Chinese Journal of Theoretical and App lied Mechanics，2010，42（3）：560－566.

［8］Gurley K R.Modelling and simulation of non-Gaussian processes［M］.University of Notre Dame，Notre Dame.1997：166－203.

［9］中華人民共和國建設(shè)部.建筑結(jié)構(gòu)荷載規(guī)范50009－2012［S］.北京：建筑結(jié)構(gòu)出版社，2012：36－50.

Prediction of wind pressure peak factor with non-gaussian simulation

LIShou-ke1，2，LIShou-ying3，CHEN Zheng-qing3，SUN Hong-xing1
（1.School of Civil Engineering，Hunan University of Science and Technology，Xiangtan 411201，China；
2.Beijing's Key Laboratory of StructuralWind Engineering and Urban Wind Environment，Beijing Jiaotong University，Beijing100044，China；
3.Wind Engineering Research Center，Hunan University，Changsha 410082，China）

Based on the Non-Gaussian Simulation of multivariate stochastic processes method，one prediction method for wind pressure extreme value was proposed.The wind pressure time histories of several opening roofs were simulated with multivariate Non-Gaussian simulation method based on wind tunnel test data.It was shown that power spectral density，coherence，deviation，skewness and kurtosis of simulated Non-Gaussian time histories are very close to the destination values.Then the peak factors ofwind pressures on the opening roofswere predicted from themultivariate non-gaussian simulation time histories for several times with the typical extreme value theory，and the results were compared with those several generalmethods.Itwas shown that Davenportmethod overestimates the positive peak factor by 60%，its skewness is negative，and it underestimates the negative peak factor by 43%；Sadek-Simiu method underestimates the peak factor by 50%，it has a higher kurtosis；the proposed method can predict the peak factor effectively，and the overall error is smallest.

non-Gaussian simulation；extreme value；peak factor；wind tunnel test

TU119＋.21

A

10.13465/j.cnki.jvs.2014.24.020

國家自然科學(xué)基金（51248001）資助；湖南省教育廳科學(xué)研究一般項(xiàng)目（14C0431）資助；湖南省高校創(chuàng)新平臺開放基金（湘教通（2012）595號）資助

2013－09－23 修改稿收到日期：2014－01－02

李壽科男，博士，講師，1981年生