黃 燦, 桂衛(wèi)華, 謝永芳, 陽春華
(中南大學 信息科學與工程學院,長沙 410083)
基于改進互相關(guān)函數(shù)的氧化鋁碳分過程多重時滯辨識
黃 燦, 桂衛(wèi)華, 謝永芳, 陽春華
(中南大學 信息科學與工程學院,長沙 410083)
為了解決復雜工業(yè)生產(chǎn)過程多重時滯辨識難題,提出一種基于改進互相關(guān)函數(shù)的多重時滯辨識方法。對于工業(yè)過程中多個受控制信號作用的過程變量,確定一個參考變量,分別考慮其他各變量和參考變量之間的相關(guān)性,選擇變量某個時間段內(nèi)的一組數(shù)據(jù)作為辨識對象,通過計算兩個變量的數(shù)據(jù)組在不同相對時間延遲對應的互相關(guān)矩陣,比較互相關(guān)矩陣的奇異值,其最大奇異值對應的滯后時間,即為所要辨識的時滯。將所提方法應用于某鋁廠連續(xù)碳分過程多重時滯的辨識,基于工業(yè)現(xiàn)場采集的生產(chǎn)數(shù)據(jù),分析變量間的關(guān)聯(lián)關(guān)系,辨識出多重時滯,然后將辨識所得的時滯代入碳分過程模型。結(jié)果表明:計算值和實測值的最大平均相對誤差僅為3.23%,驗證了所提方法辨識多重時滯的有效性。
互相關(guān)函數(shù);多重時滯;辨識;碳分過程
復雜工業(yè)生產(chǎn)過程通常由不同的生產(chǎn)工序組成,各工序的生產(chǎn)設備地域跨度大。為確保產(chǎn)品質(zhì)量,需要對各工序生產(chǎn)指標進行控制以保證工藝指標的穩(wěn)定。由于物料傳質(zhì)輸送和完成物理化學反應需要一定的時間,對工藝指標進行閉環(huán)控制會出現(xiàn)反饋信息延時。為提高產(chǎn)品質(zhì)量,在生產(chǎn)工序的不同地方,設置多個控制點,而不同控制點到工序出口處的距離不同,物料傳質(zhì)所需要的時間也不同,使得控制點的控制量到出口處的輸出量產(chǎn)生不同的時滯,因此,時滯具有多重性。
時滯辨識是復雜工業(yè)過程實現(xiàn)工藝指標閉環(huán)控制的前提條件,是控制研究的重要領(lǐng)域之一。RAD等[1]和孫建平等[2]提出互相關(guān)函數(shù)辨識方法,通過分析輸入/輸出之間的互相關(guān)函數(shù),以互相關(guān)函數(shù)的最大值辨識輸入/輸出之間的時滯。該方法根據(jù)兩個時間序列相關(guān)特性進行辨識,僅適用于單時滯辨識。SVANTE和LANNART[3]、RAMES[4]和 LIU 等[5]采用基于頻域特性的繼電反饋方法辨識時滯,但該方法需要檢測Nyquist曲線信號,而實際生產(chǎn)過程中的Nyquist曲線難以獲取。NI等[6]和孟昕元及薛東亮[7]通過小波變換的方法辨識時滯,由于實際工業(yè)生產(chǎn)過程信號變化緩慢,小波變換方法的母函數(shù)選擇困難,難以得到實際工業(yè)應用。神經(jīng)網(wǎng)絡、遺傳算法等智能辨識方法[8-10],因其算法復雜且運行時間較長,主要用于單時滯辨識,很少看到這些方法應用于多重時滯辨識的報道。
連續(xù)碳酸化分解過程(簡稱碳分過程)是氧化鋁生產(chǎn)的關(guān)鍵工序過程,具有典型的多重時滯特性。碳分過程一般由多個分解槽串級連接,由于物料反應和傳輸過程需要一定的時間,造成了控制量和輸出量之間存在很大的時滯。并且,各個槽中的物料分解程度不同,反應時間各不相同,使得各個分解槽的時滯不同,形成多重時滯。在實際碳分過程中,由于多重時滯的存在,難以實現(xiàn)分解梯度的閉環(huán)控制。因此,對碳分過程進行多重時滯的辨識,是實現(xiàn)該過程有效閉環(huán)控制、提高生產(chǎn)效率的重要前提。針對氧化鋁碳分過程,許多學者做了大量的研究工作[11-14],但對其多重時滯特性的研究較少。鄧燕妮等[15]在碳分過程機理分析的基礎(chǔ)上,建立了連續(xù)碳分過程非線性多重時滯動態(tài)數(shù)學模型,考慮了碳分過程中的多重時滯,但僅給出了時滯的經(jīng)驗取值范圍,沒有辨識出其具體值。
本文作者提出改進的互相關(guān)函數(shù),基于工業(yè)現(xiàn)場數(shù)據(jù),考查多個變量之間的時滯關(guān)系。以包含更多過程變化信息的一個時間段數(shù)據(jù)作為對象,定義互相關(guān)矩陣,研究多個變量在不同時滯下的關(guān)聯(lián)程度,然后通過計算互相關(guān)矩陣的奇異值來比較變量之間的關(guān)聯(lián)程度。因此,可以將多重時滯辨識問題轉(zhuǎn)化為求解一系列互相關(guān)矩陣的最大奇異值問題?;谔挤诌^程現(xiàn)場的運行數(shù)據(jù),運用所提方法辨識出該過程多個分解槽的變量之間的相對時滯,并將所求時滯代入鄧燕妮等[15]建立的碳分模型中驗證所提方法的有效性。
傳統(tǒng)的互相關(guān)函數(shù)法是描述兩個隨機信號在不同時刻的取值之間的相關(guān)程度,對于單輸入、單輸出系統(tǒng),一般通過輸入、輸出信號的相關(guān)性來辨識時滯。
在實際工業(yè)過程中,一個過程往往由幾個環(huán)節(jié)串聯(lián)而成,前一個環(huán)節(jié)的輸出作為后一個環(huán)節(jié)的輸入,這一類同時充當不同環(huán)節(jié)的輸出和輸入的變量,可以通過考查多個變量之間的關(guān)聯(lián)關(guān)系來辨識時滯。對于工業(yè)過程中由反應或者傳輸過程的緩慢變化造成的時滯,用一個連續(xù)時間段內(nèi)的數(shù)據(jù)辨識其大小,通過對反映變量變化過程的一組數(shù)據(jù)的辨識,更能準確獲取變量之間的關(guān)聯(lián)關(guān)系。
令數(shù)據(jù)矩陣X=[x1, x2,…, xN]表示多變量時滯系統(tǒng)中N個受控制信號作用的過程變量,其中
xi=[xi,1, xi,2, …, xi,K, xi,K+1, xi,K+2, …, xi,P-1, xi,P]T(i=1, 2, …, N; P>K)表示過程變量xi在連續(xù)P個采樣時間點的值。選取參考過程變量x1,其余N-1個變量相對于 x1的時滯表示為[τ2, τ3, …, τN]。將現(xiàn)場采集的原始數(shù)據(jù)標準化,以消除數(shù)值型屬性大小不一造成計算的偏差。
采用最大最小標準化方法對初始數(shù)據(jù)進行處理。對樣本數(shù)據(jù)集xi=[xi,1, xi,2, …, xi,P]T(i=1, 2, …, N),令ximax和ximin分別為集合xi中的最大值和最小值,則xi,l(l=1, 2, …, P)在區(qū)間[0, 1]的映射值為
將數(shù)據(jù)矩陣標準化得到
根據(jù)工業(yè)現(xiàn)場經(jīng)驗、物料傳輸和反應速度分析,可以獲得多重時滯區(qū)間[τjmin, τjmax](j=2, 3, …, N),其中,τjmax和τjmin分別為變量xj(j=2, 3, …, N)相對于x1的時滯估計值的最大值和最小值。
為分析參考變量x1和其他變量xj(j=2, 3, …, N)之間的關(guān)聯(lián)關(guān)系,取過程變量xj一個連續(xù)時間段的fj個數(shù)據(jù)K),滿足fjT≥τjmax,其中T為采樣周期,保證所取數(shù)據(jù)段能夠包含一個時滯范圍內(nèi)的有效信息。為使tn=K時能取到一組完整的數(shù)據(jù),K+fj-1≤P,即K<P+1-fj-1。當變量xj取fj個數(shù)據(jù)時,參考變量x1相應的取fj個數(shù)據(jù),K)。由于xj的變化比x1滯后時間τj,所取數(shù)據(jù)組應滿足 tn>tm。
定義x1和其他變量xj(j=2, 3, …, N)之間的互相關(guān)矩陣
2.1 氧化鋁連續(xù)碳分過程工藝分析
氧化鋁連續(xù)碳分過程是在鋁酸鈉溶液中通入CO2氣體,中和溶液中的苛性堿,使溶液的苛性比值降低,使鋁酸鈉溶液過飽和度增大,產(chǎn)生鋁酸鈉溶液自發(fā)分解的析出反應,從而析出氫氧化鋁沉淀。其主要化學反應為
實際生產(chǎn)中連續(xù)碳分過程由 6個連續(xù)分解槽組成,整個過程從原料入槽到出料一般需要 3~4 h。如圖1所示,經(jīng)脫硅送來的合格精液首先進入高位槽,從高位槽底部自壓進入1#分解槽,用低壓風提料進2#分解槽,采用同樣的方法依次將料提入后面槽。前 5臺槽子根據(jù)分解率要求通入一定量 CO2氣體進行分解,通過調(diào)節(jié)CO2氣體的通入量,控制各槽的分解率。6#槽作為出料槽,檢測合格后由出料泵打到沉降槽,沉降底流送往氫氧化鋁過濾機過濾,得到氫氧化鋁產(chǎn)品。
圖1 氧化鋁連續(xù)碳分過程工藝流程圖Fig.1 Flowsheet of continuous alumina carbonation decomposition process
由于每個分解槽內(nèi)的反應處于晶體析出的不同時期,溶液分解程度不同,導致各槽的反應時間不相同,即CO2氣體通入分解槽中到其作用到分解率的變化上所需的時間不相同,故每個槽的時滯不同,為多重時滯過程。
碳分過程每個分解槽有不同的分解率,各分解槽要滿足一定的分解梯度,才能保證末槽的分解率要求。分解率是自鋁酸鈉溶液中分解析出的氧化鋁數(shù)量與精液中所含氧化鋁數(shù)量的百分比,實際上,常用溶液分解前后的苛性比值來計算,如式(7)所示:
式中:ρ(Al2O3)J為精液中氧化鋁濃度,g/L;ρ(Al2O3)母為母液(精液和CO2氣體反應后的溶液,即1#~6#槽中的溶液)中的氧化鋁濃度,g/L;ρ(Na2OT)M為精液中的全堿濃度,g/L;ρ(Na2OT)M為母液中的全堿濃度,g/L。
2.2 碳分過程多重時滯辨識
生產(chǎn)過程中,采集了與分解率密切相關(guān)的全堿濃度和氧化鋁濃度數(shù)據(jù)計算分解率的大小,由于在連續(xù)碳分過程中,全堿濃度變化不明顯,而氧化鋁濃度變化較大,因此用氧化鋁濃度作為過程變量進行多重時滯辨識。令原液的氧化鋁濃度為x1,1#~5#槽反應后的氧化鋁濃度分別為x2、x3、x4、x5、x6,則初始數(shù)據(jù)矩陣為X=[ x1, x2, …, x6]。將原液的氧化鋁濃度x1作為參考變量,分別考查1#~5#槽的氧化鋁濃度x2、x3、x4、x5、x6與原液的氧化鋁濃度x1之間的關(guān)聯(lián)關(guān)系。
以某鋁廠的工業(yè)現(xiàn)場數(shù)據(jù)為例,選取穩(wěn)定工況條件下連續(xù)15 h的生產(chǎn)數(shù)據(jù),每10 min采樣1次,剔除生產(chǎn)異常狀況下的數(shù)據(jù),經(jīng)預處理后得到90組有效數(shù)據(jù),則P=90。對數(shù)據(jù)矩陣X進行標準化處理,得
根據(jù)現(xiàn)場經(jīng)驗,物料在每個槽中的反應和傳輸造成的時間滯后范圍是[30, 60] min。以1#分解槽為例,數(shù)據(jù)采樣周期T=10 min,取f2=6,則f2T=60 min,即取變量60 min時間段內(nèi)的數(shù)據(jù),滿足f2T≥τ2max的條件??疾樵?0 min時間段內(nèi)各變量之間的關(guān)聯(lián)程度,可得數(shù)據(jù)矩陣為:
為保證所有變量辨識的最后一組數(shù)據(jù)能取到完整的fj(j=2, 3, …, 6)個數(shù)據(jù),取K=60。辨識各槽之間的相對滯后時間,由式(4)可得當原液數(shù)據(jù)取第 1組值,原液和1#分解槽內(nèi)的氧化鋁濃度數(shù)據(jù)之間的互相關(guān)矩陣最大奇異值為
圖2所示為1#槽的 λ (x1, x2)1,tn變化圖,其中橫坐標表示 1#分解槽內(nèi)的氧化鋁濃度所取的第 tn(tn=1,2, …, 60)組數(shù)據(jù),縱坐標表示對應的互相關(guān)矩陣的最大奇異值 λ (x1, x2)1,tn。相對原液的氧化鋁數(shù)據(jù),1#槽在tn取5時,對應數(shù)據(jù)矩陣的互相關(guān)矩陣的奇異值為最大值,表示此時原液和分解槽內(nèi)的氧化鋁具有最大相關(guān)性,因此,可以辨識出1#槽相對于原液的時滯為50 min。
圖2 原液和1#槽內(nèi)氧化鋁濃度數(shù)據(jù)的相關(guān)性Fig.2 Correlation of alumina concentration data between stock solution and tank 1
粗略估計 2#~5#分解槽內(nèi)的氧化鋁濃度相對于原液的氧化鋁濃度的時滯范圍。因為每個槽的滯后時間范圍為[30, 60] min,分別取辨識數(shù)據(jù)組長度 f3=12、f4=18、f5=24、f6=30。滿足 fjT≥τjmax(j=3, 4, 5, 6)。用同樣的方法辨識2#~5#槽的時滯,圖3~6所示分別為2#~5#槽的 λ (x1, xj)1,tn(j=3, 4, …, 6; tn=2, 3, …, 60)變化圖。由圖3~6 可知,2#、3#、4#、5#分解槽在tn分別取9、14、20、25時, λ (x1, xj)1,tn為最大值。因此,可以辨識出1#~5#槽相對于原液的時滯分別為 50、90、140、200和 250 min。
為驗證改進互相關(guān)函數(shù)方法辨識時滯的準確性,以2#槽為例辨識單個槽的時滯。分析1#和2#分解槽氧化鋁數(shù)據(jù)之間的關(guān)聯(lián)關(guān)系,由圖7可知,2#槽的時滯為40 min,而由1#、2#分解槽相對于原液的時滯分別為50和90 min可知,2#槽單個槽的時滯為40 min。由此可見,兩者結(jié)果一致。
圖3 原液和2#分解槽內(nèi)氧化鋁濃度數(shù)據(jù)的相關(guān)性Fig.3 Correlation of alumina concentration data between stock solution and tank 2
圖4 原液和3#分解槽內(nèi)氧化鋁濃度數(shù)據(jù)的相關(guān)性Fig.4 Correlation of alumina concentration data between stock solution and tank 3
圖5 原液和4#分解槽內(nèi)氧化鋁濃度數(shù)據(jù)的相關(guān)性Fig.5 Correlation of alumina concentration data between stock solution and tank 4
圖6 原液和5#分解槽內(nèi)氧化鋁濃度數(shù)據(jù)的相關(guān)性Fig.6 Correlation of alumina concentration data between stock solution and tank 5
圖7 1#和2#分解槽內(nèi)氧化鋁濃度數(shù)據(jù)的相關(guān)性Fig.7 Correlation of alumina concentration data between tank 1 and tank 2
以上是原液氧化鋁取第一組數(shù)據(jù)時的相關(guān)性,當原液依次取1到60組數(shù)據(jù)時,考查原液和其他分解槽內(nèi)氧化鋁濃度數(shù)據(jù)最大相關(guān)性對應的相對滯后周期數(shù)。圖8所示為1#分解槽的數(shù)據(jù)相對于原液數(shù)據(jù)滯后周期數(shù)(D)隨原液所取數(shù)據(jù)組數(shù)變化的情況。由圖8可以看出,當原液的考查數(shù)據(jù)組由第 1組變化到第 60組時,它和1#槽內(nèi)的氧化鋁濃度的相關(guān)性最大值,基本上保持在5個周期后取得,進一步驗證了1#槽的時滯為50 min。
圖8 原液與1#分解槽的相對滯后周期變化曲線Fig.8 Relatively delay period change curve of stock solution and tank 1
將辨識出來的時滯代入鄧燕妮等[15]建立的模型,計算氧化鋁濃度值,對比模型估計濃度和實際測量濃度值,驗證時滯的精度,其計算值和實際值對比如表1所示。
由表1可以看出,1#~5#槽的氧化鋁濃度實際值和計算值之間的平均相對誤差最大值為 3.23%,表明所辨識的時滯能夠較好地反映實際碳分過程,本文作者提出的方法能夠較好地應用于復雜工業(yè)過程的多重時滯辨識。
表1 連續(xù)碳分過程各槽溶液氧化鋁濃度的實際測量值與模型計算值Table 1 Measured values and calculated values of alumina concentration of carbonation decomposition tanks (g/L)
1) 類似氧化鋁碳分過程的復雜工業(yè)生產(chǎn)過程,其多個變量具有多重時滯的特點。改進傳統(tǒng)的互相關(guān)函數(shù),從多變量中確定一個參考變量,分別考查其他變量與參考變量之間的關(guān)聯(lián)關(guān)系,選取變量某段時間內(nèi)的一組數(shù)據(jù)作為辨識對象,通過計算數(shù)據(jù)組之間的互相關(guān)矩陣,比較互相關(guān)矩陣的奇異值,其最大的奇異值對應的時滯,即為所要辨識值。
2) 基于某鋁廠的碳分過程現(xiàn)場數(shù)據(jù),分析各分解槽和原液之間氧化鋁濃度的關(guān)聯(lián)關(guān)系,從而辨識出各個分解槽相對于首槽入口處的多重時滯值。將該多重時滯值代入氧化鋁碳分過程橫型中,計算各分解槽的氧化鋁濃度,計算值和實測值的比較結(jié)果表明,計算值的精度較高。
REFERENCES
[1] RAD A B, LO W, TSANG K M. Simultaneous online identification of rational dynamics and time delay: A correlation-based approach[J]. IEEE Transaction on Control Systems Technology, 2003, 11(6): 957-959.
[2] 孫建平, 閆永躍, 于樹新, 李慶周. 時滯時變對象參數(shù)辨識方法[J]. 電光與控制, 2008, 15(1): 94-96.SUN Jian-ping, YAN Yong-yue, YU Shu-xin, LI Qing-zhou. A parameter identification algorithm for time-varying/time-delay system[J]. Electronic Optics and Control, 2008, 15(1): 94-96.
[3] SVANTE B, LENNART L. An improved phase method for time delay estimation[J]. Automatic, 2009, 45(10): 2467-2470.
[4] RAMES C P. Estimation of parameters of under-damped second order plus dead time processes using relay feedback[J].Computers and Chemical Engineering, 2006, 30(5): 832-837.
[5] LIU Tao, GAO Fu-rong. A generalized relay identification method for time delay and non-minimum phase processes[J].Automatic, 2009, 45(4): 1072-1079.
[6] NI B Y, XIAO D Y, SHAH S L. Time delay estimation for MIMO dynamical systems—With time-frequency domain analysis[J]. Journal of Process Control, 2010, 20(1): 83-94.
[7] 孟昕元, 薛東亮. 大滯后控制系統(tǒng)性能的小波變換分析方法研究[J]. 機電工程技術(shù), 2006, 35(3): 47-50.MENG Xin-yuan, XUE Dong-liang. A study of the performance of long time-delayed control system using wavelet transform analysis method[J]. Mechanical & Electrical Engineering Technology, 2006, 35(3): 47-50.
[8] 陸 燕, 杜繼紅, 李春文. 延遲時間未知的時延系統(tǒng)神經(jīng)網(wǎng)絡補償控制[J]. 清華大學學報, 1998, 38(9): 67-69.LU Yan, DU Ji-hong, LI Chun-wen. Neural networks compensate control for unknown systems with time delay[J].Journal of Tsinghua University: Sci & Tech, 1998, 38(9): 67-69.
[9] 趙仕俊, 李 遜, 左光遠. 灰色預估神經(jīng)網(wǎng)絡在時滯系統(tǒng)控制中的應用[J]. 微電子學與計算機, 2008, 25(2): 25-27.ZHAO Shi-jun, LI Xun, ZUO Guang-yuan. Research on neural network controller with grey-model prediction applied in controlling systems with uncertain time delay[J].Microelectronics & Computer, 2008, 25(2): 25-27.
[10] YADAIAH N, DEEKSHATULU B L, SIVAKUMAR L, RAO V S H. Neural network algorithm for parameter identification of dynamical systems involving time delays[J]. Applied Soft Computing, 2007, 7(3): 1084-1091.
[11] 王 志, 楊毅宏, 畢詩文, 謝雁麗. 鋁酸鈉溶液碳酸化分解過程的影響因素[J]. 有色金屬, 2002, 54(1): 43-45.WANG Zhi, YANG Yi-hong, BI Shi-wen, XIE Yan-li.Influencing factors of sodium aluminate solution carbonation decomposition[J]. Nonferrous Metals, 2002, 54(1): 43-45.
[12] 彭志宏, 李小斌, 茍中入, 劉桂華, 周秋生, 丁安平, 李光柱,李 明. 鋁酸鈉溶液碳酸化分解產(chǎn)品中的 Na2O[J]. 中國有色金屬學報, 2002, 12(6): 1285-1289.PENG Zhi-hong, LI Xiao-bin, GOU Zhong-ru, LIU Gui-hua,ZHOU Qiu-sheng, DING An-ping, LI Guang-zhu, LI Ming.Impurity Na2O in carbonization precipitation from sodium aluminate solution with high Al2O3concentration[J]. The Chinese Journal of Nonferrous Metals, 2002, 12(6): 1285-1289.
[13] 李小斌, 陳 濱, 周秋生, 劉桂華, 彭志宏, 劉祥民. 鋁酸鈉溶液碳酸化分解過程動力學[J]. 中國有色金屬學報, 2004, 14(5):848-853.LI Xiao-bin, CHEN Bin, ZHOU Qiu-sheng, LIU Gui-hua,PENG Zhi-hong, LIU Xiang-min. Kinetics of carbonation decomposition of sodium aluminate solution[J]. The Chinese Journal of Nonferrous Metals, 2004, 14(5): 848-853.
[14] 胡志坤, 桂衛(wèi)華, 陽春華. 基于遺傳算法的連續(xù)碳酸化分解過程模糊優(yōu)化[J]. 金屬材料與冶金工程, 2007, 35(5): 34-37.HU Zhi-kun, GUI Wei-hua, YANG Chun-hua. Fuzzy optimization for the process of continuous carbonation decomposition based on genetic algorithm[J]. Metal Materials and Metallurgy Engineering, 2007, 35(5): 34-37.
[15] 鄧燕妮, 桂衛(wèi)華, 陽春華, 謝永芳. 氧化鋁碳酸化分解動態(tài)過程建模及非線性分析[J]. 中國有色金屬學報, 2008, 18(9):1736-1741.DENG Yan-ni, GUI Wei-hua, YANG Chun-hua, XIE Yong-fang.Dynamic process modeling and nonlinear analysis for alumina carbonation decomposition[J]. The Chinese Journal of Nonferrous Metals, 2008, 18(9): 1736-1741.
Multi-delays identification for alumina carbonation decomposition process based on improved cross-correlation function
HUANG Can, GUI Wei-hua, XIE Yong-fang, YANG Chun-hua
(School of Information Science and Engineering, Central South University, Changsha 410083, China)
In order to resolve the problem of multi-delays identification for complex industrial process, an improved cross-correlation function was proposed. In all the process variables affected by control signals, the reference variable was selected by considering the correlation with the other variables respectively. For the considered variable, a set of data in a continuous time segment sampled was selected as identification object. The cross-correlation matrix for the data sets of the reference variable and the other variables was calculated. By comparing the singular values of cross-correlation matrix, the delay corresponding to the maximum singular value was the required delay. The proposed method was applied to identifying the multi-delays of alumina carbonation decomposition process using the field data. At last, the identified delays were substituted to alumina carbonation decomposition process model. The results show that the maximum average relative error between the calculated and tested results is only 3.23%, and the proposed multi-delays identification method based on improved cross-correlation function is effective.
cross-correlation function; multi-delays; identification; carbonation decomposition process
TF821
A
1004-0609(2011)05-1186-06
國家自然科學基金資助項目(60634020, 61074117)
2010-11-22;
2011-03-15
桂衛(wèi)華,教授,博士;電話:0731-88830765;E-mail: gwh@mail.csu.edu.cn
(編輯 何學鋒)