張 飛，徐光黎，魏志云，朱可俊

(1.中國地質大學a.巖土鉆掘與防護教育部工程研究中心;b.工程學院，武漢 430074;2.中國水電顧問集團成都勘測設計研究院，成都 610072)

1 研究背景

在地下工程設計與施工中，初始地應力狀態(tài)(自重應力和構造應力)及巖體力學性態(tài)參數(shù)的不確定性和復雜性，使得通過數(shù)值方法計算得到的圍巖應力、變形及塑性圈深度等與現(xiàn)場監(jiān)控量測值大相徑庭。巖土體的不連續(xù)性、各向異性，尺寸效應等因素導致通過室內外試驗確定的巖土體參數(shù)具有一定的局限性。單一的測試手段獲取的巖體參數(shù)可靠性不足，一般不能作為計算的有效原始數(shù)據(jù)。通常需要從不同的試驗方法中來綜合確定工程巖體的力學參數(shù)。20世紀70年代以位移量測為基礎的位移反分析方法的提出［1］，為獲取有效的巖(土)體力學參數(shù)提供了一種新的途徑。

根據(jù)地下洞室開挖引起的位移進行反分析，無論開挖過程中采取什么措施，因爆破損傷、應力調整、擾動等原因產(chǎn)生的松動圈是不可避免的。松動圈內的巖體性質已不同于開挖前巖體及松動圈以內巖體性質，其內部裂隙因開挖卸荷、擾動具有明顯的非連續(xù)性。而松動圈的存在導致位移量測值的增大，對反分析結果必然有較大的影響［2］。近年來考慮松動圈的位移反分析，也有不少學者研究，如江權［3］等考慮了松動圈厚度利用增量位移實現(xiàn)了高地應力條件下巖體力學參數(shù)的彈脆塑性反分析;李寧［4］等考慮松動圈影響探討了彈塑性位移反分析;陳秋紅［5］等研究了松動圈模型及其在數(shù)值分析中對松動圈的處理方法，并應用于錦屏一級水電站反饋分析中;李鴻博［6］研究了在隧道開挖中應用位移反分析時考慮圍巖松動圈的必要性，論述了雙介質模型在位移反分析中的實現(xiàn)。

精細反演分析是當今反分析研究的努力目標，本文結合前人考慮松動圈模型反饋分析研究成果，針對大崗山水電站地下廠房洞室群施工期圍巖開挖卸荷變形特征，考慮圍巖變形時間與空間分布特征，利用開挖過程現(xiàn)場監(jiān)測增量位移和松動圈深度信息，基于進化神經(jīng)網(wǎng)絡算法［7－10］和 FLAC3D差分程序，實現(xiàn)大崗山水電站地下廠房洞室群施工期巖體參數(shù)的精細反演，為后續(xù)施工開挖、反饋設計與施工提供合理的參數(shù)。

2 松動圈反演基本模型

在位移反分析中考慮松動圈的影響，就需要建立松動圈反饋分析模型［5，11－12］，主要包括 4 個內容:松圈巖體材料模型、松動圈深度確定、松動圈分區(qū)分級和松動圈巖體力學參數(shù)。

2.1 松動圈材料模型

松動圈深度內巖體變形主要由圍巖在爆破震動，應力重分布過程中形成塑性圈的變形。這部分巖體彈性模量降低、黏聚力減小、連續(xù)性降低，強度顯出一定程度的弱化效應。數(shù)值計算采用的是連續(xù)介質模型，為了簡單地考慮松動圈的影響，采用參數(shù)弱化的等效連續(xù)介質模擬圍巖松動圈巖體。

2.2 松動圈深度確定

松動圈深度的確定方法［12］主要依靠現(xiàn)場測試的手段，其中主要包括聲波測試法、多點位移計法、地質雷達法、地震波法、電阻率法和滲透法等。

聲波在巖體的傳播特性，在速度和振幅上都有所響應。振幅的衰減，取決于巖體對聲波的吸收作用。通常彈性波在這類巖體中傳播較快，如堅硬巖體;節(jié)理不發(fā)育和風化程度低的巖體;孔隙率小、密度大、彈性模量大的巖體;抗壓強度大的巖體;斷層和破碎帶少或其規(guī)模小的巖體。反之，彈性波的傳播速度較慢。因此，可根據(jù)巖體這一特性，利用聲波測試確定松動圈深度。其步驟如下:

(1)確定圍巖爆破開挖前的基準波速vp0;

(2)根據(jù)規(guī)范［13］用波速方法判斷爆破巖體破壞的標準，確定巖體爆破開挖后波速vp和基準波速vp0的變化率，即η=1－vp/vp0;

(3)根據(jù)經(jīng)驗，取η＞10%的圍巖范圍定為松動圈深度，對于致密脆性巖石，η=10%為破壞臨界值;對節(jié)理發(fā)育、變形具有塑性特征的巖石，取η=20%為破壞臨界值。

2.3 松動圈分區(qū)分級模型

根據(jù)聲波測試曲線和上述判定標準，將聲波曲線劃分為2段，第1段從洞壁表面到某一深度H0，該深度范圍稱為松弛深度，即縱波波速開始趨于穩(wěn)定時對應的測點深度(見圖1)，這一區(qū)間聲波逐漸上升;H0以內為第2段，聲波保持未擾動的原巖基準值，這樣把圍巖劃分為2個區(qū)域，即松動區(qū)域和未擾動區(qū)域。

然后，再根據(jù)松動圈深度－波速曲線的趨勢和形狀，應用回歸分析建立松動圈深度－波速函數(shù)關系。

設E，E0分別為松動圈和未擾動原巖區(qū)巖體的彈性模量。根據(jù)平面波在無限大彈性體中的傳播波速可以推到出

圖1 松動圈示意圖Fig.1 Sketch of the loose zone of surrounding rock

式中，H∈(0，H0］，α定義為松動圈內巖體的松弛系數(shù);當H=H0時，α等于1表示圍巖未松動。再根據(jù)松弛系數(shù)α值不同對松動圈進行分級，確定各級松動圈厚度，將同級松動圈深度用等值線連成封閉區(qū)域，即形成松動圈的分區(qū)分級模型。

2.4 松動圈巖體力學參數(shù)

松動圈巖體作為一種弱化參數(shù)的等效連續(xù)介質體，對松動圈內巖體變形(變形參數(shù))與破壞(強度參數(shù))起主要作用的力學參數(shù)一般通過反饋分析來確定。采用隨α值連續(xù)線性變化的弱化參數(shù)模型，巖體參數(shù)基于如下假設，即:

式中:c，c0為松動圈和原巖區(qū)巖體的黏聚力;φ，φ0為松動圈和原巖區(qū)巖體的內摩擦角。當H=H0時，松動圈參數(shù)等同原巖參數(shù)。根據(jù)松動圈巖體參數(shù)是沿圍巖徑向坐標的單值連續(xù)函數(shù)，可確定松動圈任一點巖體弱化參數(shù)。結合對巖體參數(shù)的敏感性研究，因此確定待反演的參數(shù)有波速梯度t，彈性模量E0，峰值黏聚力 c0，峰值內摩擦角 φ0。

3 施工期分層開挖模擬

3.1 工程與監(jiān)測簡介

建設中的大崗山水電站位于四川省雅安市石棉縣大渡河中游的中高山峽谷地區(qū)，是大渡河干流近期開發(fā)的大型水電工程之一，共有4臺發(fā)電機組，全部采用地下廠房方案。沿大渡河左岸分別布置主廠房、主變室和尾水調壓室三大廠房，軸線整體方向N55°E，垂直埋深 390～520 m，水平埋深 310～530 m。主副廠房開挖尺寸為226.6 m×30.0 m×73.8 m(長×寬 ×高)，屬于大型的地下廠房洞室群。

為工程安全與實現(xiàn)監(jiān)控量測、反饋設計與施工，三大洞室布置了大量的多點位移計。在主副廠房1#，2#，3#和 4#機組中心線斷面(如圖2)，主變室及尾水調壓室機組斷面也均布置了大量多點位移計。另外為了精確了解洞室群高邊墻松動圈的厚度，在主副廠房、主變室和尾水調壓室還布置了聲波檢測斷面，如圖3，典型聲波曲線如圖4。

本文對大崗山水電站硐室開挖第6層測試的30條聲波數(shù)據(jù)進行分析，典型聲波曲線如圖4。本文采用分段線性函數(shù)通過最小二乘擬合松動圈深度－波速關系，表示如下:

圖2 典型監(jiān)測剖面及多點位移計布置Fig.2 Typical monitored profile and layout of multi-point displacement meters

圖3 聲波測試剖面及檢測孔布置Fig.3 Profile of sound testing and layout of testing holes

圖4 典型聲波速度曲線Fig.4 Typical acoustic velocity curve

式中:t定義為波速梯度，反映了波速沿洞壁向圍巖深部變化快慢的程度，對廠房30條聲波曲線進行最小二乘擬合，可得各常數(shù)的取值范圍如表1。

表1 地下廠房松動圈線性最小二乘擬合常數(shù)Table 1 Fitted least square constants of the loose zone of underground powerhouse

3.2 施工過程與反演分析實測數(shù)據(jù)確定

主廠房共分9層開挖(包含尾水連接管部分)，模擬分層開挖施工過程可劃分為:

(1)初始步(計算初始地應力場)。

(2)第Ⅰ層開挖(第1施工步第1增量步)，假設中導洞開挖和擴挖間隔很短，幾乎同時開挖完成。

(3)第Ⅰ層開挖完噴錨支護(第1施工步第2增量步)，實際開挖過程中，開挖和噴錨支護時間相隔較短，假設每層開挖完立即設施噴錨支護，在每個增量步內，時間Δt都認為是從零開始的。

(4)隨后各層開挖，均采用上述相同的施工步與增量步，通過各層的開挖可以一次測得已開挖各層設置的多點位移計對應的每一施工步或增量步的位移。對于大多數(shù)地下洞室采取分層開挖，用增量位移作為實測反分析數(shù)據(jù)具有兩點好處:①增量位移參與的反演計算可以不用考慮位移計安裝前的測前位移損失;②大多數(shù)巖體都有不同程度的時效變形特征(尤其對于軟巖)，位移計讀數(shù)據(jù)并不是圍巖變形的最終位移，采用增量位移反演無需估計圍巖的最終變形量。

反分析的巖體參數(shù)一般來說是代表整個地下洞室群變形與強度特征的“綜合參數(shù)”，即要求實測數(shù)據(jù)也能充分體現(xiàn)這一點，主要從以下2個方面選取實測數(shù)據(jù):①考慮數(shù)據(jù)點的時間分布特征，選取不同開挖層的典型位移;②考慮數(shù)據(jù)點的空間分布，選取主廠房和主變室不同開挖部位的數(shù)據(jù)，包括頂拱、拱肩、邊墻及巖錨梁等部位。實測數(shù)據(jù)如表2。

表2 用于反演的實測位移及時空分布特征Table 2 Measured displacement for back analysis and its space-time distribution characteristics

3.3 基于FLAC3D的計算模型

模型計算域范圍:上邊界取至地表，下邊界取至高程718.0 m位置，左邊界距主廠房頂拱中心線200 m，右邊界距主變室頂拱中心線200 m。模型以計算域左、下邊界交點作坐標原點，x軸指向S35°E方向，與主廠房軸線垂直，z軸豎直向上。上部為自由邊界，左、右邊界限制水平方向位移，底邊界限制水平和豎直位移約束。計算模型如圖5。

FLAC3D模型中，圍巖與輝綠巖脈用四節(jié)點等參單元模擬，選用摩爾－庫倫(Mohr-Coulomb)彈塑性本構關系和強度判據(jù)?？紤]在開挖過程中的及時支護，初期支護形式:掛鋼筋網(wǎng)Φ8@15 cm×15 cm，噴15 cm厚C25混凝土及Φ28長6 m和Φ32長9 m相間的系統(tǒng)錨桿。掛網(wǎng)噴混凝土采用殼(shell)單元并以彈性模型進行模擬，系統(tǒng)錨桿采用錨索(cable)單元模擬。計算模型材料參數(shù)如表3。

圖5 FLAC3D數(shù)值計算模型Fig.5 Numerical calculation model of FLAC3D

表3 計算模型圍巖及支護結構力學參數(shù)Table 3 Mechanical parameters for the calculation model of surrounding rock and supporting structures

4 巖體參數(shù)的智能反演

4.1 目標函數(shù)的建立

本文選取現(xiàn)場實測位移增量值(如表1)與神經(jīng)網(wǎng)絡映射值的殘差平方和作為目標函數(shù)，也是遺傳算法的適應度函數(shù)，即

式中:X=(E0，c0，φ0，t)為需要確定的參數(shù)向量;為參數(shù)約束條件，即待反演參數(shù)應滿足本構關系理論規(guī)則和參數(shù)允許的取值范圍;Δuj和分別為第j測點某一施工步或某幾個施工步前后位移增量的網(wǎng)絡映射值和實測值;優(yōu)化反分析的目標就是求得使f(X)取得最小值時的參數(shù)值。

4.2 基于GA－BP算法的位移反分析

將神經(jīng)網(wǎng)絡和遺傳算法結合［7］，即利用了神經(jīng)網(wǎng)絡高度非線性映射、網(wǎng)絡推理和預測的功能，又利用了遺傳算法全局尋優(yōu)特性，克服了神經(jīng)網(wǎng)絡學習陷入局部最小問題，在處理巖體參數(shù)與目標變量之間無顯示表達式的復雜問題中，具有較高的應用價值。本文GA－BP主要步驟如下:

(1)基于正交試驗生成神經(jīng)網(wǎng)絡學習樣本和測試樣本，通過FLAC3D正算程序獲取樣本輸出，并進行歸一化處理。

(2)利用搜索得到的網(wǎng)絡進行樣本學習，建立巖土體力學參數(shù)與輸出位移之間的映射關系。

(3)對訓練成熟的網(wǎng)絡進行初始化設置，并確定待反演參數(shù)的取值范圍。

(4)在待反演的參數(shù)取值范圍內隨機生成可能參數(shù)群體，代入訓練好的神經(jīng)網(wǎng)絡進行位移預測，計算目標函數(shù)值(適應度)。

(5)若適應度滿足要求，則當前參數(shù)即為最優(yōu)參數(shù);若適應度不滿足要求，進行遺傳進化操作(復制、雜交、變異)，產(chǎn)生下一子代的參數(shù)群體。

(6)重復進行(4)、(5)步操作，直到獲得滿足目標函數(shù)精度的最優(yōu)參數(shù)組合。

4.3 BP神經(jīng)網(wǎng)絡樣本構造

采用正交試驗設計BP神經(jīng)網(wǎng)絡樣本，正交試驗具有均勻分散、整齊可比的特點，使試驗樣本具有充分的代表性［14］。

本試驗因素為4個待求參數(shù)，因素分為5個均勻變化水平，為使網(wǎng)絡具有較強的泛化能力。采用正交表L25(56)分別構造網(wǎng)絡的學習樣本和測試樣本，參數(shù)輸入部分如表4左部分，對于樣本輸出則采用FLAC3D正算獲取監(jiān)測點位移增量，如表4右部分。其中，前20個樣本作為網(wǎng)絡的學習樣本，后5個為測試樣本。

4.4 參數(shù)的反演計算

在Matlab7中編寫算法，通過GA搜索ANN結構，獲得最優(yōu)網(wǎng)絡結構為4－34－28－6;采用學習率 η=0.20、動量因子 α=0.45、利用學習樣本訓練網(wǎng)絡結構，選取測試樣本的系統(tǒng)誤差極小值對應的網(wǎng)絡連接權值，得到具有最佳泛化能力的網(wǎng)絡結構。

然后設置遺傳代數(shù)Igen=100、種群規(guī)模Np=20、選擇變異概率0.09，在Matlab7中調用gaot工具箱，經(jīng)GA搜索得到最優(yōu)參數(shù)，見表5。

5 反演分析結果檢驗

反演分析得到巖體參數(shù)不是主要目的，最終目的是利用獲取的參數(shù)進行后續(xù)開挖的預測(位移、應力、應變、松動圈深度等)和進行圍巖施工期穩(wěn)定性評價，最后反饋設計與施工。

反演得到的巖體參數(shù)是否可用還有待于驗證，本文利用獲取參數(shù)通過正算，獲取圍巖測點位移，借用灰色系統(tǒng)理論中的后驗差法［15］進行位移檢驗來判斷反演所得參數(shù)的可信性。

表4 BP神經(jīng)網(wǎng)絡學習樣本和測試樣本Table 4 Learning samples and testing samples for BP neural network

表5 反分析所得巖體參數(shù)Table 5 Rock mass parameters obtained from back analysis

對廠房第Ⅶ層開挖引起的增量位移進行預測，并與實測值比較，兩者數(shù)值較接近，計算后驗差比值和小誤差概率表明反演的巖體參數(shù)能夠較好代表整個地下洞室群的圍巖特征，如表6。說明了反分析結果是可用的。

表6 第Ⅶ層開挖位移增量計算值與實測值對比Table 6 Comparison between calculated and measured value of displacement increment of the 7th excavation layer

6 結論與建議

建立了模擬實際開挖的三維計算模型，進行了巖體參數(shù)的增量位移智能反演，得出以下結論:

(1)建立考慮松動圈參數(shù)弱化效應的分層開挖FLAC3D計算模型來模擬實際施工情況，基于正交試驗方法獲得神經(jīng)網(wǎng)絡學習和測試樣本，通過遺傳算法搜索BP神經(jīng)網(wǎng)絡結構，形成進化神經(jīng)網(wǎng)絡(GA-BP)算法，可以迅速地建立圍巖參數(shù)與位移之間的高度非線性映射關系，避免了參數(shù)調整時需進行的“繁瑣”數(shù)值計算，也減少了計算時間。

(2)合理選取反演參數(shù)和考慮空間效應現(xiàn)場位移增量數(shù)據(jù)，運用遺傳算法全局尋優(yōu)特性，利用訓練好的參數(shù)－位移之間非線性映射關系網(wǎng)絡，成功反演出與彈塑性模型相關3個巖體參數(shù)和波速梯度。

(3)運用反演參數(shù)，通過正算對主廠房第Ⅶ層開挖引起位移增量進行了預測，并進行了位移的后驗差法檢驗，結果表明反演參數(shù)的適用性，可作為后續(xù)開挖施工、反饋設計及圍巖穩(wěn)定性評判依據(jù)。

(4)本文采用松動圈巖體參數(shù)線性弱化效應的模型來考慮松動圈對參數(shù)反演的影響，參數(shù)弱化效應根據(jù)現(xiàn)場聲波測試曲線的擬合來確定，由于聲波曲線形狀很大程度上取決于圍巖質量，采用線性變化的分段函數(shù)擬合有一定的局限性，更具一般規(guī)律的松動圈參數(shù)弱化效應模型有待于進一步的研究。

［1］劉志寶，寧海龍，張海濤.巖土工程中位移反分析方法及其進展［J］.路基工程，2006，(5):7－10.(LIU Zhibao，NING Hai-long，ZHANG Hai-tao.Displacement Back-Analysis Methods and Progress in Geotechnical Engineering［J］.Subgrade Engineering，2006，(5):7 －10.(in Chinese))

［2］楊志法，熊順成，王存玉.關于位移反分析的某些考慮［J］.巖石力學與工程學報，1995，14(1):11－16.(YANG Zhi-fa，XIONG Shun-cheng，WANG Cun-yu.Some Considerations of the Back-Analysis on Displacement［J］.Chinese Journal of Rock Mechanics and Engineering，1995，14(1):11 －16.(in Chinese))

［3］江 權，馮夏庭，蘇國韶，等.基于松動圈-位移增量監(jiān)測信息的高地應力下洞室群巖體力學參數(shù)的智能反分析［J］.巖石力學與工程學報，2007，26(增1):2654－2662.(JIANG Quan，F(xiàn)ENG Xia-ting，SU Guo-shao，et al.Intelligent Back Analysis of Rock Mass Parameters for Large Underground Caverns under High Earth Stress Based on EDZ and Increment Displacement［J］.Chinese Journal of Rock Mechanics and Engineering，2007，26(Sup.1):2654 －2662.(in Chinese))

［4］李 寧，段小強，陳方方，等.圍巖松動圈的彈塑性位移反分析方法探索［J］.巖石力學與工程學報，2006，25(7):1304 － 1308.(LI Ning，DUAN Xiao-qiang，CHEN Fang-fang，et al.A Back Analysis Method for Elasto-plastic Displacement of Broken Rock Zone Around Tunnel［J］.Chinese Journal of Rock Mechanics and Engineering，2006，25(7):1304 －1308.(in Chinese))

［5］陳秋紅，李仲奎，張志增.松動圈分區(qū)模型及其在地下工程反饋分析中的應用［J］.巖石力學與工程學報，2010，29(增 1):3216 － 3220.(CHEN Qiu-hong，LI Zhong-kui，ZHANG Zhi-zeng.Loosening Zonal Model and Its Application to Back Analysis of Underground Engineering［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29(Sup.1):3216 －3220.(in Chinese))

［6］李鴻博.考慮松動圈的位移反分析在隧道工程中的應用［D］.上海:同濟大學，2006.(LI Hong-bo.Application of Displacement Back Analysis to Tunnel Engineering in Consideration of Loose Zone［D］.Shanghai:Tongji University，2006.(in Chinese))

［7］黃 戡，劉寶琛，彭建國，等.基于遺傳算法和神經(jīng)網(wǎng)絡的隧道圍巖位移智能反分析［J］.中南大學學報(自然科學版)，2011，42(1):213 －219.(HUANG Kan，LIU Bao-chen，PENG Jian-guo，et al.Intelligent Back-Analysis of Tunnel Surrounding Rock Displacement Based on Genetic Algorithm and Neural Network［J］.Journal of Central South University(Science and Technology)，2011，42(1):213 －219.(in Chinese))

［8］王樹棟.基于進化神經(jīng)元算法的堡鎮(zhèn)隧道軟弱圍巖施工彈塑性智能位移反分析［J］.北京交通大學學報，2010，34(3):112 － 116.(WANG Shu-dong.Elasto-Plastic Intelligence Displacement Back Analysis of Soft Rock Construction in Baozhen Tunnel Based on Evolutionary-Neuron Algorithm［J］.Journal of Beijing Jiaotong U-niversity.2010，34(3):112 －116.(in Chinese))

［9］馮夏庭，張治強，楊成祥，等.位移反分析的進化神經(jīng)網(wǎng)絡方法研究［J］.巖石力學與工程學報，1999，18(5):529 － 533.(FENG Xia-ting，ZHANG Zhi-qiang，YANG Cheng-xiang，et al.Study on Genetic Neural Network Method of Displacement Back Analysis［J］.Chinese Journal of Rock Mechanics and Engineering，1999，18(5):529 －533.(in Chinese))

［10］關永平，宋 建，王述紅，等.基于GA-BP算法的隧道圍巖力學參數(shù)反分析［J］.東北大學學報(自然科學版)，2012，33(2):276 － 278.(GUAN Yong-ping，SONG Jian，WANG Shu-hong，et al.Back Analysis of Mechanical Parameters of Surrounding Rocks Based on GA-BP Algorithm［J］.Journal of Northeastern University(Natural Science)，2012，33(2):276 － 278.(in Chinese))

［11］張建海，胡著秀，楊永濤，等.地下廠房圍巖松動圈聲波擬合及監(jiān)測反饋分析［J］.巖石力學與工程學報，2011，30(6):1191 －1197.(ZHANG Jian-hai，HU Zhuxiu，YANG Yong-tao，et al.Acoustic Velocity Fitting and Monitoring Feedback Analysis of Surrounding Rock Loosing Zone in Underground Powerhouse［J］.Chinese Journal of Rock Mechanics and Engineering，2011，30(6):1191 －1197.(in Chinese))

［12］鄒紅英，肖 明.地下洞室開挖松動圈評估方法研究［J］.巖石力學與工程學報，2010，29(3):513－519.(ZOU Hong-ying，XIAO Ming.Study of Methodology for Assessment of Excavation Disturbed Zone of Underground Caverns［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29(3):513 －519.(in Chinese))

［13］DL/T5099—2011.水工建筑物地下開挖工程施工技術規(guī)范［S］.北京:水利電力出版社，1999.(DL/T5099—2011.The Professional Standards Compilation Group of People’s Republic of China.Technical Specifications of Underground Excavation of Hydraulic Engineering［S］.Beijing:China Water Power Press，1999.(in Chinese))

［14］方開泰.均勻設計與均勻設計表［M］.北京:科學出版社，1994.(FANG Kai-tai.Uniformity Design and Design Tables［M］.Beijing:Science Press，1994.(in Chinese))

［15］鄧聚龍.灰色系統(tǒng)的GM模型［M］.武漢:華中理工大學出版社，1986.(DENG Ju-long.GM Model of Grey System［M］.Wuhan:Huazhong University of Science and Technology Press，1986.(in Chinese))