梅 寒 陳炳乾,2 王正帥 高 建 余 昊
(1.江蘇師范大學地理測繪與城鄉(xiāng)規(guī)劃學院,江蘇徐州221116;2.長安大學地質(zhì)工程與測繪學院,陜西西安710054)
概率積分法模型以隨機介質(zhì)理論為基礎發(fā)展而來,是我國礦區(qū)沉陷預計方面應用最為廣泛的方法之一,其參數(shù)取值的準確性直接關系到預計的精度[1]。隨著非線性模型方程最優(yōu)解求解方法的不斷發(fā)展,在開采沉陷領域出現(xiàn)了大量概率積分模型參數(shù)的求解算法,目前主要有線性近似法、試驗設計法和智能優(yōu)化算法[2]。
線性近似法是將非線性函數(shù)進行線性化處理,典型的方法為最小二乘法。該方法求參精度較高,但僅適合矩形工作面,且易受到初值影響導致結(jié)果難以收斂。郭廣禮等[3]將最小二乘法與穩(wěn)健估計模型相結(jié)合,通過迭代定權進行參數(shù)反演,克服了反演結(jié)果不準確或不收斂的問題。王曉輝等[4]比較了13種常用穩(wěn)健估計方法在參數(shù)反演中的穩(wěn)健特性,結(jié)果表明:對于概率積分模型參數(shù)的穩(wěn)健估計,L1法、German-Mclure方法和IGGIII方法具有較強的可靠性和穩(wěn)健性。王友等[5]針對最小二乘法擬合求參時易受到病態(tài)矩陣和異常點干擾而導致發(fā)散的問題,提出了基于抗差嶺估計的求參模型。通過L曲線法確定嶺參數(shù),然后用最小二乘法擬合求參,在一定程度上減小了病態(tài)矩陣和粗差對求參精度的影響。高超等[6]利用Matlab語言編寫了基于實測數(shù)據(jù)曲線擬合的概率積分模型參數(shù)求解程序,通過將實際工作面劃分為多個不同煤層厚度工作面疊加來進行沉陷預計參數(shù)求取,并對模型進行了優(yōu)化。試驗設計法是一種研究多因子、多水平的設計方法,通過挑選具有代表性的組合進行試驗并根據(jù)試驗結(jié)果確定最優(yōu)組合。由于概率積分模型中影響因素較多,導致試驗設計量很大、效率較低、計算機編程實現(xiàn)困難,因此適用性不強。查劍鋒等[7]針對正交設計在因素選取的水平數(shù)較多時試驗次數(shù)過多的問題,提出了均勻設計法,相比于正交試驗不僅減少了試驗次數(shù),還減小了初值選取不當?shù)挠绊?。智能?yōu)化算法是近些年興起的一系列模擬生物在自然界中的進化規(guī)律或自然法則的算法。吳侃等[8]、葛家新等[9]利用模矢法求參,能夠較為準確地反演模型參數(shù),但在初值不同時會出現(xiàn)即使擬合情況理想,但參數(shù)反演結(jié)果仍然相差很大的情況。查劍鋒等[10]將遺傳算法應用于概率積分模型參數(shù)反演,證明了該算法相較于最小二乘法和模矢法在反演準確性和可靠性方面優(yōu)勢明顯。徐孟強等[11]利用粒子群算法求取了河北省某礦區(qū)的概率積分模型參數(shù),通過將第一次反演結(jié)果回代能夠獲取更高的求參精度。賈新果[12]利用蟻群算法求取了平頂山某礦區(qū)工作面的概率積分模型參數(shù),走向和傾向擬合中誤差分別為51 mm和23 mm。陳濤等[13]利用果蠅算法對安徽省某礦區(qū)工作面的概率積分法參數(shù)進行了求取,下沉預計值與實測值的均方差為33.7 mm。王正帥等[14]將兩種算法進行了融合提出了文化-粒子群算法,將算法收斂的成功率提高到1,并對兗州礦區(qū)某工作面的概率積分法參數(shù)進行了求取,擬合相對誤差為最大下沉值的1.1%。
模矢法、遺傳算法、文化-粒子群算法、粒子群算法、果蠅算法和蟻群算法都已被證明了能夠應用于概率積分法模型參數(shù)反演,但由于各算法的理論基礎、模型各不相同,導致在參數(shù)反演過程中表現(xiàn)出不同的性能[15]。為了進一步分析上述6種算法在概率積分模型參數(shù)反演過程中的性能,本研究利用Matlab語言編寫了各算法的概率積分模型參數(shù)反演程序。通過模擬試驗,分別從反演準確性、穩(wěn)定性、抗誤差干擾能力、全局尋優(yōu)能力和運行效率5個方面進行對比分析,并通過工程實例進行驗證。
隨機介質(zhì)理論的主要思想為:隨著工作面開采,采空區(qū)上方的地表各單元之間產(chǎn)生相對移動,彼此不存在聯(lián)系即巖層中的顆粒體介質(zhì)之間的移動是隨機的[15]。我國學者劉寶琛、廖國華等在此基礎上發(fā)展了概率積分法[16],該方法理論體系成熟,計算機編程易實現(xiàn),因此被廣泛應用于礦區(qū)地表沉陷預計。
概率積分法共有5個預計參數(shù),分別為下沉系數(shù)q、主要影響角正切tanβ、水平移動系數(shù)b、開采影響傳播角θ以及拐點偏移距s。根據(jù)概率積分法原理,假設在傾斜煤層中開采單元B(u ,v)引起的地表點A(x ,y)的下沉值為
式中,W0=mqcosα為礦區(qū)充分采動時的最大下沉值,m;m為煤層采厚,m;α為煤層傾角,(°);D2為工作面傾向長度,m;D1為工作面走向長度,m。
6種智能優(yōu)化算法反演流程如圖2所示。
(1)模矢法。模矢法又稱步長加速法,是一種解決非線性函數(shù)最優(yōu)解的算法。該算法的主要思想是首先給定一個初值作為最優(yōu)解并計算函數(shù)值;然后選擇初始解的一個基點B1在其左右相鄰位置以等步長計算目標函數(shù)值并與上述最優(yōu)函數(shù)值比較,以目標函數(shù)值較小的點為迭代的臨時矢點;最后分別沿坐標軸方向搜索最優(yōu)點并記為新的基點B2。當?shù)鷿M足設置的終止步長時,迭代結(jié)束,此時的最優(yōu)解即為非線性函數(shù)尋優(yōu)的最優(yōu)解。該算法運行效率較高,但需要給定參數(shù)初值,適用于參數(shù)大致范圍已知的礦區(qū)。算法流程如圖2(a)所示。本研究參考文獻[17]設定終止步長 e=0.01,0.01,0.01,1,1,表示 q、tanβ、θ、b與s終止精度分別為0.01、0.01、0.01、1和1。
(2)遺傳算法。遺傳算法是一種模擬自然進化過程的隨機搜索算法[18],其思想來源于“遺傳選擇”和“物競天擇”。該算法主要思路是將優(yōu)化問題比作一個種群,將待求參數(shù)編碼為染色體,然后通過迭代方式進行選擇、交叉以及變異等步驟,最終生成符合目標要求的染色體。遺傳算法通用性高、魯棒性強,適合于多參數(shù)模型求解,但局部尋優(yōu)能力不足,算法流程如圖2(b)所示。本研究參考文獻[18]并通過多次試驗對比,設定迭代次數(shù)上限為100,種群數(shù)量為50,交叉率為0.9,變異率為0.001。
(3)文化-粒子群算法。文化算法提供了一種多進化過程的計算模型框架,以粒子群算法作為種群空間,以隨機粒子群算法為信念空間的進化算法,分別從種群空間和信念空間模擬文化的雙重演化過程,最大化地利用粒子群算法的開發(fā)能力和隨機粒子群算法的開拓能力,形成獨立進化且相互促進的雙重演化機制。該算法集成了隨機粒子群算法的大范圍、高效率搜索和粒子群算法的局部精細化搜索的優(yōu)點,收斂成功率為1,適用于參數(shù)范圍未知的礦區(qū)參數(shù)反演,但計算機編程實現(xiàn)較為復雜,算法流程如圖2(c)所示。本研究參考文獻[14]并通過多次試驗對比,設定慣性權重ω=0.4,學習因子c1=c2=2,迭代步數(shù)、允許接受代數(shù)、允許影響代數(shù)、最大迭代步數(shù)和目標尋優(yōu)值分別設置為100、2、6、200和0.001。
(4)粒子群算法。粒子群算法是一種模擬鳥群覓食行為的進化計算方法,其主要思想是將優(yōu)化問題看作覓食的鳥群,則“食物”為優(yōu)化問題的最優(yōu)解,鳥群中的每只覓食的“鳥”就是解空間中的一個粒子,每個粒子的位置都是一個潛在的最優(yōu)解。通過分享和學習粒子的位置信息不斷更新自身的位置和速度,直到滿足終止條件獲得最優(yōu)解。粒子群算法的優(yōu)勢在于個體粒子能夠充分利用個體經(jīng)驗和群體經(jīng)驗調(diào)整自身位置,尋優(yōu)速度較快,適用于下沉規(guī)律符合正態(tài)分布的礦區(qū)參數(shù)反演,易收斂;缺點是易早熟收斂,局部尋優(yōu)能力較差。算法流程如圖2(d)所示。本研究參考文獻[11]并通過多次試驗對比,設定迭代次數(shù)上限為200,種群數(shù)量為20,慣性權重ω=0.95,學習因子c1=c2=2,粒子的飛行速度-1≤V≤1。
(5)果蠅優(yōu)化算法。果蠅算法是一種模擬果蠅覓食行為的尋優(yōu)算法。果蠅的嗅覺和視覺非常敏銳,能夠發(fā)現(xiàn)數(shù)十千米以外的氣味源,在靠近氣味源后敏銳的視覺能夠發(fā)現(xiàn)食物和同伴的位置,然后快速地朝目標飛去。果蠅算法的優(yōu)點是尋優(yōu)精度較高、計算機編程易實現(xiàn),對于下沉規(guī)律不完全符合正態(tài)分布的礦區(qū)較為適用;缺點是當步長設定過大時,局部尋優(yōu)能力會下降,設定過小時,則全局尋優(yōu)能力較差。算法流程如圖2(e)所示。本研究參考文獻[19]并通過多次試驗對比,設定迭代次數(shù)上限為100,種群數(shù)量為20,步進值為0.01。
(6)蟻群算法。蟻群算法是一種模擬螞蟻覓食行為的隨機搜索優(yōu)化算法。其思想來源于螞蟻可以自發(fā)地找到巢穴到目標的最短路徑,并且能適應環(huán)境的改變發(fā)現(xiàn)新的最短路徑,同時在通過的路徑上留下信息素以供后面的螞蟻選擇。蟻群算法相較于其他啟發(fā)式算法,具有較強的魯棒性、計算機編程易實現(xiàn),適用于參數(shù)范圍已知的礦區(qū)參數(shù)反演;缺點是若參數(shù)設置不當,會導致求解速度慢且結(jié)果質(zhì)量差,并且收斂速度慢,易陷入局部最優(yōu)解。算法流程如圖2(f)所示。本研究參考文獻[12]并通過多次試驗對比,最終設定螞蟻移動次數(shù)上限為100,蟻群數(shù)量為30,信息素揮發(fā)系數(shù)為0.9,轉(zhuǎn)移概率常數(shù)為0.2。
以某礦區(qū)地質(zhì)采礦條件為背景,模擬工作面走向長度為D1=400 m,傾向長度為D2=200 m,煤層采厚3 m,平均采深H=300 m,煤層傾角5°,頂板管理方法為全部垮落法。地表沉陷預計參數(shù)為:下沉系數(shù)q=0.8,主要影響角正切tan β=2.5,開采影響傳播角θ=85°,水平移動系數(shù)b=0.25,拐點偏移距s1=s2=s3=s4=20 m,各參數(shù)設計值見表1。在工作面上方沿走向布設了30個監(jiān)測點,沿傾向布設了20個監(jiān)測點,模擬工作面與監(jiān)測點的位置關系如圖3所示。
以預計值作為參數(shù)反演的起始數(shù)據(jù),結(jié)合3.1節(jié)工作面信息,分別利用模矢法、遺傳算法、文化-粒子群算法、粒子群算法、果蠅算法和蟻群算法進行概率積分法參數(shù)反演。為防止試驗結(jié)果出現(xiàn)偶然性,對以上算法分別運行了10次并取平均值作為最終反演結(jié)果,如圖4(a)和圖4(b)所示,各算法擬合中誤差如圖4(c)所示。
由圖4(a)和圖4(b)可知:各算法反演概率積分法參數(shù)精度相當,下沉系數(shù)q、主要影響角正切tanβ、開采影響傳播角θ、水平移動系數(shù)b、拐地偏移距s的相對誤差分別不超過0.3%、0.5%、1.3%、1%和10%。由圖3(c)可知:模矢法的預計值與實測值的擬合中誤差最小,為1.4 mm,反演準確性最高,果蠅算法的擬合中誤差最大,為4.8 mm,反演準確性最低。
由于模矢法受參數(shù)初值的影響,運算方法固定,不具有隨機性,因此不將模矢法納入比較,利用3.1節(jié)的數(shù)據(jù),對各種算法分別運行了50次,得到各參數(shù)的反演結(jié)果波動情況如圖5所示。
分析圖5可知:各算法在50次試驗中的參數(shù)反演值均存在波動,由于概率積分法參數(shù)無法單獨求取,因此通過分析各算法在50次試驗中預測值與實測值的擬合中誤差來判斷算法的穩(wěn)定性;在50次試驗中,遺傳算法的預計值與實測值擬合中誤差平均值最小,為7.8 mm,穩(wěn)定性最強,果蠅算法的擬合中誤差平均值最大,為16.6 mm,穩(wěn)定性最差。
為了比較各算法對隨機誤差的抗干擾能力,在模擬數(shù)據(jù)中分別加入中誤差為5、10、20 mm的隨機誤差,再利用各算法進行概率積分法參數(shù)反演,參數(shù)反演結(jié)果如圖6所示。
分析圖6可知:在加入隨機誤差后,反演值與真值的誤差逐漸增大,在誤差不超過20 mm時,下沉系數(shù)q、主要影響角正切tanβ、開采影響傳播角θ、水平移動系數(shù)b、拐地偏移距s的相對誤差分別不超過0.8%、1.3%、1.1%、3.82%和12.61%;在分別加入3種不同量級的隨機誤差下,模矢法的預計值和實測值的擬合中誤差平均值為7.1 mm,抗隨機誤差干擾能力最強,蟻群算法的擬合中誤差平均值為16.3 mm,抗誤差干擾能力最弱。
由文獻[20]可知:粗差數(shù)量通常為觀測值數(shù)量的5%~10%,因此本研究選取監(jiān)測點數(shù)量的8%的粗差加入原始數(shù)據(jù)中,即粗差的數(shù)量不超過4個。根據(jù)文獻[8]可知,實測數(shù)據(jù)在拐點處和最大下沉點處的誤差會對反演結(jié)果造成一定影響,因此分別在最大下沉點和拐點處加入200 mm的粗差,并利用各算法對改造后的數(shù)據(jù)進行反演,參數(shù)反演結(jié)果如圖7所示。
分析圖7可知:各算法對于粗差出現(xiàn)位置的敏感程度不同,其中模矢法、粒子群算法、果蠅算法和蟻群算法對于拐點位置的粗差更敏感,文化-粒子群算法和遺傳算法對于最大下沉點的粗差更敏感。當粗差數(shù)量不超過3個時,粒子群算法預計值與實測值的擬合中誤差最小,為21.7 mm,抗粗差干擾能力最強,模矢法擬合中誤差最大,為29.6 mm,抗粗差干擾能力最弱。
為了防止試驗結(jié)果具有偶然性,在擴大尋優(yōu)范圍(表2)后對各程序分別運行了10次并取平均值作為最終結(jié)果,試驗得到參數(shù)反演結(jié)果如圖8所示。
分析圖8可知:模矢法的反演準確性相對于其他算法大大降低,原因是模矢法受迭代初值的影響較大,當?shù)踔到咏硞€局部最優(yōu)解的時候,模矢法難以跳出局部最優(yōu)解。遺傳算法、文化-粒子群算法、粒子群算法、果蠅算法和蟻群算法都有較好的跳出局部解的能力;遺傳算法預計值與實測值的擬合中誤差最小,為5.6 mm,全局搜索能力最強,模矢法最弱。
本研究試驗計算機處理器為Intel(R)Core(TM)i5-3470,處理器CPU主頻為3.20 GHz,運行內(nèi)存為8 GB,各算法的運行時間如表3所示。
由于模矢法需要給定初值,迭代初期參數(shù)擬合值就在真值附近,因此運算效率最高;文化-粒子群的雙重演化機制避免了尋優(yōu)過程發(fā)散的問題,同時還避免了反復試算迭代初值的盲目性,因此運算效率較高;遺傳算法、粒子群算法、果蠅算法和蟻群算法是利用初始群體迭代尋找最優(yōu)解,通過優(yōu)化代碼結(jié)構和調(diào)整系統(tǒng)參數(shù)能夠提高運行效率,但是程序基礎運行時間是固定的,優(yōu)化代碼結(jié)構和調(diào)整系統(tǒng)參數(shù)能提高的運算效率有限,因此這4種算法的運行效率較低。
將6種算法分別用于陜西省大柳塔礦區(qū)52304工作面進行概率積分法參數(shù)反演。該礦區(qū)位于陜西省神木縣,煤層平均厚度6.94 m,平均采深235 m,煤層傾角為1°~3°。工作面沿傾向線布設了25個監(jiān)測點記為TQ1~TQ25,點間距為25 m。礦區(qū)位置與觀測線位置如圖9所示。
各算法初值及參數(shù)范圍設置如表4所示。
為了防止試驗結(jié)果具有偶然性并比較算法穩(wěn)定性,利用6種算法分別進行10次概率積分法參數(shù)反演,取其平均值作為反演結(jié)果,如表5所示,實測值與預計值的擬合結(jié)果如圖10和表6所示。
從表6和圖10可知:各算法參數(shù)反演結(jié)果由于各自理論基礎與模型的不同而略有差異,實測值與預計值擬合情況較好;模矢法由于初值設置與真值相差較大,因此反演結(jié)果陷入局部解,準確性較差,該結(jié)果與模擬試驗結(jié)果相吻合;遺傳算法擬合中誤差最小,全局尋優(yōu)能力與穩(wěn)定性最強,與模擬試驗結(jié)果一致;運行效率方面模矢法有加速向最優(yōu)點移動的特性,因此尋優(yōu)效率最高,與模擬試驗結(jié)果一致;文化-粒子群算法的雙重演化機制避免了錯誤值的反復迭代試算,效率也較高,其他4種算法通過初始種群迭代尋優(yōu),效率較低,該結(jié)果與模擬試驗結(jié)論相吻合。
(1)利用Matlab語言編寫了模矢法、遺傳算法、文化-粒子群算法、粒子群算法、果蠅算法和蟻群算法的運算程序,通過模擬試驗和工程實例對比分析了上述幾種算法在概率積分法參數(shù)反演中的性能。試驗表明:雖然各算法的反演參數(shù)準確性均較高,但在求參過程中的表現(xiàn)各不相同,其中模矢法運行效率最高但易陷入局部解,遺傳算法全局尋優(yōu)能力和穩(wěn)定性最強,文化-粒子群算法各方面表現(xiàn)平均,運行效率方面較其他智能優(yōu)化算法高,粒子群算法具有一定的抗誤差干擾能力,但易早熟收斂,果蠅算法和蟻群算法在抗誤差性能方面表現(xiàn)較弱。
(2)當參數(shù)初值與真值相差較大時,模矢法反演準確性會大大降低。從算法反演準確性和運行效率方面綜合考慮,當參數(shù)真值或范圍已知時,模矢法為最優(yōu)選擇。當參數(shù)真值或范圍未知時,最優(yōu)算法的選擇依次為文化-粒子群算法、遺傳算法、果蠅算法、粒子群算法、蟻群算法和模矢法。