郭昌放，武 祥，楊 真，陳一鼎，馬留柱，馬中原

(1.中國礦業(yè)大學(xué) 人工智能研究院，江蘇 徐州 221116；2.中國礦業(yè)大學(xué) 礦業(yè)工程學(xué)院，江蘇 徐州 221116；3. 昆明煤炭科學(xué)研究所，云南 昆明 650051；4. 中煤科工集團(tuán)西安研究院有限公司，陜西 西安 710000；5. 中國礦業(yè)大學(xué)(北京) 能源與礦業(yè)工程學(xué)院，北京 100083)

工作面地質(zhì)異常的精準(zhǔn)預(yù)測是實(shí)現(xiàn)煤礦智能精準(zhǔn)開采的重要基礎(chǔ)和前提，對于推動(dòng)煤炭資源的安全高效開采具有重要意義[1-2]。煤礦工作面電磁波CT(Computerized Tomography)技術(shù)作為一種非接觸式探測方法，具有施工便捷、效率高、分辨能力強(qiáng)、成本較低、無破壞性等優(yōu)點(diǎn)，已成為目前煤礦工作面內(nèi)地質(zhì)異常探測的主要方法之一[3]。

煤礦工作面電磁波層析成像技術(shù)來自于醫(yī)學(xué)上常用的CT技術(shù)，屬于圖像重構(gòu)技術(shù)的應(yīng)用之一。在利用電磁波CT技術(shù)對工作面內(nèi)地質(zhì)異常進(jìn)行層析反演計(jì)算時(shí)，需要求解的是一個(gè)病態(tài)矩陣方程。針對該矩陣方程的求解方法，最為經(jīng)典的是GORDEN[4]提出的代數(shù)重建算法(ART)、GILBERT[5]提出的聯(lián)合迭代重建算法(SIRT)及其優(yōu)化算法[6-7]。然而這些常見的層析反演算法大都采用單一的初始模型，通過反復(fù)向改進(jìn)的鄰近解集移動(dòng)，直到滿足允許條件以求得最終解。其反演結(jié)果很大程度上取決于單一初始模型的選擇，而且容易陷入局部最優(yōu)解。隨著計(jì)算機(jī)計(jì)算能力的提高以及人工智能技術(shù)的快速發(fā)展，先進(jìn)的智能計(jì)算方法已經(jīng)被廣泛應(yīng)用于各種復(fù)雜模型問題求解過程中。

自1960年以來，蒙特卡洛法、BP神經(jīng)網(wǎng)絡(luò)、遺傳算法、粒子群等智能算法通過模仿生物的生存和進(jìn)化規(guī)則而逐步發(fā)展起來，并為復(fù)雜多維度問題的處理提供了新的思路和方法。由于這些算法對初始模型的依賴程度較低，其強(qiáng)大的全局搜索能力和特殊的進(jìn)化搜索機(jī)制引起了各領(lǐng)域?qū)W者的廣泛關(guān)注[8-9]。1985年以來，這些智能算法開始應(yīng)用于地球物理反演問題中，特別是在地震反演領(lǐng)域取得較多的進(jìn)展[10-11]。隨后，智能計(jì)算方法也逐漸在大地電磁反演領(lǐng)域得到推廣和應(yīng)用。然而智能算法在工作面電磁波CT層析反演領(lǐng)域的應(yīng)用和研究較少。針對目標(biāo)函數(shù)極值的求解，蒙特卡洛法在樣本數(shù)量比較少的情況下，會朝著局部最優(yōu)解方向收斂，但在求解參數(shù)維度較大的情況下缺乏實(shí)用性[12]。神經(jīng)網(wǎng)絡(luò)算法在實(shí)際應(yīng)用過程中難以解決多參量的優(yōu)化求解問題。遺傳算法是一種具有高效、并行、全局搜索能力的智能算法，在進(jìn)化搜索的過程中，能夠通過之前積累的進(jìn)化知識自主調(diào)整進(jìn)化方向，求得目標(biāo)函數(shù)的最優(yōu)解。遺傳算法作為一種實(shí)用、高效、魯棒性強(qiáng)的優(yōu)化算法，由于其本身的開放特性，具有較好的改進(jìn)潛能，發(fā)展非常迅速，是目前應(yīng)用較為廣泛的智能尋優(yōu)搜索算法[13-15]。因此，筆者將遺傳算法及其優(yōu)化算法引入工作面電磁波CT地質(zhì)異常反演重構(gòu)過程，對于工作面內(nèi)地質(zhì)異常的精準(zhǔn)識別具有重要意義。

工作面掘進(jìn)和回采期間時(shí)刻都在揭露地質(zhì)信息，地質(zhì)數(shù)據(jù)的產(chǎn)生在時(shí)間和空間上是持續(xù)發(fā)展和動(dòng)態(tài)變化的過程。在工作面電磁波CT層析反演的過程中，如果能夠充分結(jié)合工作面在掘進(jìn)及回采過程中揭露的掘進(jìn)數(shù)據(jù)、鉆探數(shù)據(jù)以及回采數(shù)據(jù)等多源先驗(yàn)地質(zhì)信息作為約束，對于提高工作面內(nèi)地質(zhì)異常反演的精準(zhǔn)性具有重要意義[16]。筆者建立了基于多種群自適應(yīng)遺傳算法和多源數(shù)據(jù)約束的地質(zhì)異常層析反演模型，通過函數(shù)測試和數(shù)值模擬對提出方法的可行性和有效性進(jìn)行了驗(yàn)證，并成功應(yīng)用于工程探測資料處理，為工作面內(nèi)地質(zhì)異常的精準(zhǔn)反演和預(yù)測提供了一種新思路。

1 工作面層析反演模型

電磁波CT探測方法主要是利用不同介質(zhì)對電磁波吸收能力的差異，根據(jù)采集到電磁波信號的強(qiáng)度對傳播區(qū)域內(nèi)的異常介質(zhì)進(jìn)行預(yù)測。電磁波在同一介質(zhì)中傳播的過程中，其輻射強(qiáng)度的高低主要與發(fā)射功率以及天線的組合形式有關(guān)。因此，在選擇天線的形式時(shí)需要綜合考慮探測所需的最佳匹配形式，從而使得輻射的能量達(dá)到最大化。在井下探測時(shí)，受到工作面現(xiàn)場施工條件的限制，通?？紤]方向因子較簡單的天線，因此選擇磁偶極子天線激發(fā)產(chǎn)生電磁波，如圖1所示。

圖1 磁偶極子天線輻射場Fig.1 Radiation field of magnetic dipole antenna

如圖1所示，當(dāng)電磁波輻射源在各向同性、均勻的介質(zhì)中傳播時(shí)，若觀測點(diǎn)P到輻射源點(diǎn)A的距離為r(觀測點(diǎn)在輻射場內(nèi))，則點(diǎn)P處的電磁波場強(qiáng)E為

(1)

其中，E為傳播介質(zhì)中距離發(fā)射點(diǎn)r位置的實(shí)測場強(qiáng)，μV；E0為電磁波發(fā)射端產(chǎn)生的初始場強(qiáng)，μV；β為傳播介質(zhì)的吸收系數(shù)， dB/m；R為發(fā)射點(diǎn)距接收點(diǎn)之間的距離，m；sinθ為方向性因子；θ為發(fā)射端天線軸與電磁波采集點(diǎn)方向之間的夾角。對于工作面電磁波CT探測的現(xiàn)場空間條件下，θ的數(shù)值可近似為90°，則場強(qiáng)E的計(jì)算公式變?yōu)?/p>

(2)

在煤礦工作面進(jìn)行電磁波CT探測時(shí)，接收機(jī)采集到的場強(qiáng)值都是H(dB)。因此，式(2)經(jīng)過分貝轉(zhuǎn)換[17]可得

H=H0-βr-lnr

(3)

式中，H為實(shí)測記錄的場強(qiáng)值，dB；H0為發(fā)射初始場強(qiáng)值，dB。

電磁波層析反演計(jì)算時(shí)，首先將煤礦工作面目標(biāo)探測區(qū)域進(jìn)行網(wǎng)格離散，假設(shè)共劃分成B個(gè)網(wǎng)格。其中，每個(gè)網(wǎng)格稱為一個(gè)像素，并用xj代表第j個(gè)網(wǎng)格的真實(shí)吸收系數(shù)，j=1,2,3,…,B。

假設(shè)射線的總條數(shù)為N，則從發(fā)射端穿過探測區(qū)域到達(dá)接收端的任一條射線yi(i=1,2,3,…,N)，在第j個(gè)網(wǎng)格的截距為dij，則在第i條射線路徑上有

(4)

由式(3)可知

Hi=H0-yi-lnri

(5)

即

(6)

式中，Hi為第i條射線采集到的電磁波場強(qiáng)數(shù)據(jù)；ri為第i條射線的長度。

工作面電磁波CT探測過程是在多個(gè)發(fā)射點(diǎn)發(fā)射電磁波，其中針對每一個(gè)發(fā)射點(diǎn)激發(fā)產(chǎn)生的電磁波在多個(gè)接收點(diǎn)進(jìn)行數(shù)據(jù)采集，則式(4)可進(jìn)一步推廣為

(7)

即

D×X=Y

(8)

式中，D為A×B階系數(shù)矩陣，即每條射線通過每個(gè)網(wǎng)格的截距，其中射線條數(shù)為N，網(wǎng)格數(shù)目為B；X為B×1階吸收系數(shù)矩陣；Y為A×1階實(shí)測常數(shù)矩陣。

由于煤礦井下工作面現(xiàn)場探測條件的限制，射線的條數(shù)一般遠(yuǎn)小于網(wǎng)格的數(shù)目，因此需要求解的是一個(gè)具有病態(tài)特征的稀疏矩陣的解。由于方程(7)沒有唯一解，而目前常見的ART,SIRT等傳統(tǒng)層析反演算法大多依賴于問題的線性化形式，采用單一的初始模型，通過反復(fù)向改進(jìn)的鄰近解集移動(dòng)，反演的結(jié)果很大程度上取決于初始值的選擇。針對該問題，將式(7)矩陣求解問題轉(zhuǎn)化為一個(gè)泛函極值求解問題，并建立目標(biāo)函數(shù)反演模型，進(jìn)而通過具有較強(qiáng)全局搜索能力的智能算法對目標(biāo)函數(shù)進(jìn)行求解。

在電磁波層析反演過程中，目標(biāo)函數(shù)的定義為觀測數(shù)據(jù)與理論數(shù)據(jù)之間的差值。首先，根據(jù)式(5)得出電磁波在傳播過程中場強(qiáng)的實(shí)測損耗值；然后，對探測范圍內(nèi)介質(zhì)網(wǎng)格賦予隨機(jī)的初始吸收系數(shù)，進(jìn)而根據(jù)式(4)，計(jì)算得到電磁波在傳播過程中場強(qiáng)的理論損耗值；最后，通過智能算法反復(fù)更新每個(gè)網(wǎng)格的吸收系數(shù)值，直到目標(biāo)函數(shù)值取得最小或迭代次數(shù)達(dá)到設(shè)定上限時(shí)，即滿足停止條件，輸出當(dāng)前時(shí)刻每個(gè)網(wǎng)格的吸收系數(shù)數(shù)值作為反演結(jié)果。具體目標(biāo)函數(shù)定義為

(9)

因此，地質(zhì)異常層析反演模型可定義為：求X′∈C，使得

f(X′)=minf(X)

(10)

2 多種群自適應(yīng)遺傳算法反演方法

2.1 基本遺傳算法

基本遺傳算法(Single Genetic Algorithm，SGA)是一種借鑒生物進(jìn)化過程中適者生存的自然選擇規(guī)律，從而演變形成的一種具有較強(qiáng)全局搜索能力的智能算法。SGA作為一種全局優(yōu)化的隨機(jī)搜索智能算法，其搜索性能的優(yōu)劣主要決定于帶有概率性質(zhì)的交叉算子和變異算子[18]。由于SGA往往采用單一固定的遺傳參數(shù)，因此遺傳參數(shù)的選擇尤為重要。應(yīng)用SGA解決不同的具體問題時(shí)，不同行業(yè)的學(xué)者提出了不同遺傳參數(shù)的取值方法[19]，但到目前為止，仍然沒有一個(gè)統(tǒng)一的參數(shù)選擇規(guī)則。當(dāng)面對一個(gè)新的優(yōu)化問題時(shí)，需要經(jīng)過大量的測試才能從中選擇出最優(yōu)的結(jié)果。然而這不僅需要消耗大量的時(shí)間，而且得到的結(jié)果也并非全局最優(yōu)解。

因此，利用SGA智能算法對建立的電磁波層析反演目標(biāo)函數(shù)模型求解時(shí)，為了盡可能避免早熟收斂現(xiàn)象的出現(xiàn)，筆者在SGA的基礎(chǔ)上，通過引入多種群協(xié)同進(jìn)化算子以及自適應(yīng)調(diào)整算子對其進(jìn)行優(yōu)化，在保證種群多樣性的同時(shí)(全局搜索性能)，盡可能提高搜索的精準(zhǔn)性(局部搜索性能)，以獲得更加精確的層析反演結(jié)果。

2.2 基本遺傳算法的優(yōu)化

2.2.1多種群遺傳算法

多種群遺傳算法(Multiple Population Genetic Algorithm，MPGA)是在SGA的基礎(chǔ)上，通過具有不同遺傳參數(shù)的多個(gè)種群協(xié)同進(jìn)化搜索的方式，使得每個(gè)種群朝著不同方向進(jìn)化，從而具有更大范圍和更加精細(xì)化的進(jìn)化搜索能力。同時(shí)，通過移民算子實(shí)現(xiàn)了多個(gè)種群協(xié)同進(jìn)化過程中優(yōu)質(zhì)信息的交流，而且最終獲得的結(jié)果是由多個(gè)種群協(xié)同進(jìn)化綜合搜索所得。MPGA通過多個(gè)具有不同遺傳參數(shù)的種群協(xié)同進(jìn)化搜索，使得SGA的全局搜索性能得到一定提高。

2.2.2自適應(yīng)遺傳算法

交叉概率Pc和變異概率Pm的選擇是影響SGA進(jìn)化搜索性能的關(guān)鍵指標(biāo)。Pc的取值越大，越容易產(chǎn)生新個(gè)體。但是如果Pc的取值過大時(shí)，種群中的個(gè)體特別是高質(zhì)量的個(gè)體越容易被破壞，從而導(dǎo)致結(jié)果難以收斂。如果Pc的取值過小時(shí)，新個(gè)體就越不容易產(chǎn)生，從而使得整個(gè)搜索進(jìn)化過程變得更加緩慢，甚至停滯不前。另一方面，對于變異概率Pm的取值來講，如果Pm的取值太小，就不容易產(chǎn)生新的個(gè)體基因組合結(jié)構(gòu)。但是，如果Pm的取值太大，那么遺傳算法就等同于一種隨機(jī)搜索算法。經(jīng)過多年理論和實(shí)踐發(fā)現(xiàn)，SGA在進(jìn)化到不同的階段需要不同的參數(shù)設(shè)置[20-21]。自適應(yīng)遺傳算法(Adaptive Genetic Algorithm，AGA)能夠使得種群在進(jìn)化過程中對應(yīng)交叉概率Pc和變異概率Pm的取值隨著適應(yīng)度值的不同而自適應(yīng)的發(fā)生變化。對于適應(yīng)度值較大的個(gè)體(距離目標(biāo)函數(shù)最優(yōu)解較近的個(gè)體)，選擇較低的Pc和Pm取值，從而使得該個(gè)體得以保護(hù)，并進(jìn)入下一代中；而對于那些適應(yīng)度值較小的個(gè)體(距離目標(biāo)函數(shù)最優(yōu)解較遠(yuǎn)的個(gè)體)，往往選擇較高的Pc和Pm取值，從而使得這些個(gè)體更容易被破壞并產(chǎn)生新的個(gè)體。

因此，AGA能夠很好地解決SGA交叉概率Pc和變異概率Pm取值問題，在SGA搜索進(jìn)化過程中，能夠時(shí)刻保持種群的多樣性，同時(shí)也有利于保證種群的收斂性能，進(jìn)一步提高了SGA的局部搜索能力。

2.2.3多種群自適應(yīng)遺傳算法

MPGA智能算法具有較強(qiáng)的全面搜索能力，而AGA智能算法能夠?qū)崟r(shí)有效的調(diào)整種群的進(jìn)化方向，因此將這2種思想進(jìn)行整合，通過設(shè)計(jì)多種群自適應(yīng)遺傳算法(Multi-population Adaptive Genetic Algorithm,MAGA)能夠顯著提高SGA的搜索性能。另一方面，通過設(shè)置精英個(gè)體(所有種群中適應(yīng)度值最高的個(gè)體)并將其單獨(dú)保留，不參與遺傳操作，防止其結(jié)構(gòu)被破壞，保證了算法的全局收斂能力。其具體實(shí)現(xiàn)流程如圖2所示。

圖2 多種群自適應(yīng)遺傳算法進(jìn)化流程Fig.2 Evolutionary flow chart of MAGA

(1)參數(shù)初始化，其中，MP為種群的數(shù)量；NP為種群中個(gè)體的數(shù)量；Nmax為種群最大遺傳代數(shù)；Ggap為代溝；K為種群二進(jìn)制編碼位數(shù)；Pc1為多種群交叉概率上限；Pc2為多種群交叉概率下限；Pm1為多種群變異概率上限；Pm2為多種群變異概率下限；Ngen為當(dāng)前進(jìn)化代數(shù)。

(2)初始化多個(gè)種群Porigin(Ngen=0)，則

Porigin={Pop(Ngen,1),Pop(Ngen,2),…, Pop(Ngen,p),…,Pop(Ngen,MP)}

(11)

其中，Pop(Ngen,p)為在第Ngen代進(jìn)化過程中第p個(gè)種群。

Pop(Ngen,p)={P(Ngen,p,1),P(Ngen,p,2),…,P(Ngen,p,n),…,P(Ngen,p,NP)}

(12)

其中，P(Ngen,p,n)為在第Ngen代進(jìn)化過程中第p個(gè)種群中第n個(gè)個(gè)體。

(13)

其中，xj(Ngen,p,n)為在第Ngen代進(jìn)化過程中第p個(gè)種群中第n個(gè)個(gè)體中第j個(gè)變量。

(3)個(gè)體評價(jià)，計(jì)算每個(gè)種群中每個(gè)個(gè)體的適應(yīng)度值F(Ngen,p,n)，其中，F(xiàn)(Ngen,p,n)為進(jìn)化到第Ngen代時(shí)，第p個(gè)種群中第n個(gè)個(gè)體的適應(yīng)度值。

Ngen=1,2,3,…,Naxg；p=1,2,3,…,MP；n=1,2,3，…，NP。

(4)根據(jù)設(shè)置的Ggap進(jìn)行選擇操作。

(5)交叉操作。其中，第p個(gè)種群在第Ngen代進(jìn)化過程中，交叉概率Pc(Ngen,p)的自適應(yīng)變化范圍為[Pc2(p),Pc1(p)]，則

Pc1(p)=Pc2+(Pc1-Pc2)δ

(14)

Pc2(p)=Pc1(p)-0.3

(15)

各種群交叉概率定義為

(16)

其中，δ為0～1的隨機(jī)數(shù)；Fc(Ngen,p)為第Ngen代進(jìn)化過程中第p個(gè)種群中要進(jìn)行交叉操作的2個(gè)個(gè)體中較大的適應(yīng)度值；Favg(Ngen,p)為第Ngen代進(jìn)化過程中第p個(gè)種群的平均適應(yīng)度值；Fmax(Ngen,p)為第Ngen代進(jìn)化過程中第p個(gè)種群中最大適應(yīng)度值。

(6)變異操作。所有種群中第p個(gè)種群在第Ngen代進(jìn)化過程中，變異概率Pm(Ngen,p)的自適應(yīng)變化范圍為[Pm2(p),Pm1(p)]，則

Pm1(p)=Pm2+(Pm1-Pm2)δ

(17)

Pm2(p)=Pm1(p)/10

(18)

各種群變異概率定義為

(19)

其中，F(xiàn)m(Ngen,p)為第Ngen代進(jìn)化過程中第p個(gè)種群中變異個(gè)體的適應(yīng)度值。

(7)移民算子。不同種群之間都互相獨(dú)立，主要通過移民算子實(shí)現(xiàn)種群間的信息互換。具體的互換規(guī)則為：將目標(biāo)種群中質(zhì)量最差的個(gè)體替換為源種群中質(zhì)量最好的個(gè)體。

(8)人工選擇操作。每一代進(jìn)化結(jié)束后，由人工算子選擇出每個(gè)種群中的最優(yōu)個(gè)體，然后通過對比選擇出所有種群中最優(yōu)個(gè)體進(jìn)行單獨(dú)存儲。由于精英個(gè)體不進(jìn)行選擇、交叉和變異等遺傳操作，以確保不同群體在進(jìn)化過程中產(chǎn)生的最優(yōu)個(gè)體不被破壞或丟失。

(9)判斷是否達(dá)到終止條件。如果未達(dá)到，則返回到第(3)步；如果達(dá)到，輸出結(jié)果。

2.3 算法測試

為了測試MAGA智能算法的進(jìn)化搜索能力，首先通過一個(gè)典型的含有多個(gè)極值的一元函數(shù)g(x)為例進(jìn)行測試，如圖3所示。分別通過SGA和MAGA在[1,2]內(nèi)尋找以下目標(biāo)函數(shù)的最小值：

(20)

由圖3可以看出，該目標(biāo)函數(shù)含有多個(gè)極小值。為了求得該目標(biāo)函數(shù)g(x)的全局最小值，并檢驗(yàn)MAGA的搜索性能，分別通過SGA,MPGA,AGA和MAGA按照表1中的參數(shù)設(shè)置，進(jìn)行重復(fù)50次重復(fù)搜索計(jì)算，得到的進(jìn)化過程和進(jìn)化結(jié)果如圖4和表2所示。

圖3 目標(biāo)函數(shù)分布示意Fig.3 Distribution diagram of objective function

由圖4和表2可以看出，上述4種不同尋優(yōu)算法經(jīng)過50次重復(fù)計(jì)算發(fā)現(xiàn)，SGA算法出現(xiàn)16次早熟收斂現(xiàn)象，平均尋優(yōu)結(jié)果最高為-0.792；MPGA和AGA算法分別出現(xiàn)5次和3次早熟收斂現(xiàn)象，AGA出現(xiàn)3次早熟收斂現(xiàn)象；而MAGA未出現(xiàn)早熟收斂現(xiàn)象，平均尋優(yōu)結(jié)果最低為-0.869。

表1 參數(shù)設(shè)置

圖4 SGA,MPGA,AGA和MAGA智能算法50次 重復(fù)進(jìn)化搜索結(jié)果Fig.4 Search results of SGA,MPGA,AGA and MAGA intelligent algorithms with 50 times of repeated evolution

表2 50次重復(fù)計(jì)算平均尋優(yōu)結(jié)果

因此，通過對具有多極值特征的一元函數(shù)g(x)的求解可知，在SGA算法的基礎(chǔ)上，MAGA算法通過多種群協(xié)同進(jìn)化算子和遺傳參數(shù)自適應(yīng)調(diào)節(jié)算子能夠更穩(wěn)定的搜索到更為精確的全局最優(yōu)解，有效提高了SGA算法的全局和局部搜索能力。

3 多源數(shù)據(jù)約束模型

為了進(jìn)一步提高工作面內(nèi)地質(zhì)異常反演預(yù)測的精準(zhǔn)性，在MAGA層析反演智能算法的基礎(chǔ)上提出了范圍約束、就近約束和平均值約束3種約束模型，使得巷道、鉆探以及回采揭露的先驗(yàn)約束數(shù)據(jù)融入到MAGA智能算法的層析反演搜索計(jì)算過程中。

3.1 范圍約束

在煤礦工作面已揭露區(qū)域內(nèi)，對采用MAGA智能算法進(jìn)行層析反演計(jì)算得到的每個(gè)網(wǎng)格的吸收系數(shù)值進(jìn)行約束處理。其中，在無地質(zhì)異常揭露的區(qū)域內(nèi)，若x′j小于當(dāng)前通過反演得到的所有網(wǎng)格的平均吸收系數(shù)值xavg，則x′j并不發(fā)生改變；若該區(qū)域中x′j大于xavg，則通過在變量的下限Ll與xavg之間隨機(jī)產(chǎn)生一個(gè)數(shù)值來代替x′j。

對于在有地質(zhì)異常揭露的區(qū)域內(nèi)，若x′j大于xavg，則x′j并不發(fā)生改變；若x′j小于xavg，則通過在變量的上限Lu與xavg之間隨機(jī)產(chǎn)生一個(gè)數(shù)值來代替x′j。

在范圍約束條件下，基于MAGA智能算法的層析反演數(shù)學(xué)模型為

(21)

(22)

其中，f為目標(biāo)函數(shù)；x′j為探測區(qū)域第j個(gè)網(wǎng)格對應(yīng)的反演吸收系數(shù)值；xavg為探測區(qū)域所有網(wǎng)格對應(yīng)反演的層析反演吸收系數(shù)平均值；δ為0～1的隨機(jī)數(shù)；Q1為正常區(qū)域(無異常構(gòu)造揭露區(qū)域)內(nèi)大于xavg的網(wǎng)格序號集合；Q2為異常區(qū)域(有異常構(gòu)造揭露區(qū)域)內(nèi)小于xavg的網(wǎng)格序號集合。

3.2 就近約束

在煤礦工作面已揭露區(qū)域范圍內(nèi)，對采用MAGA智能算法進(jìn)行層析反演計(jì)算得到的每個(gè)網(wǎng)格的吸收系數(shù)值x′j進(jìn)行約束處理。其中，在無地質(zhì)異常揭露的區(qū)域內(nèi)，若x′j小于xavg，則x′j并不改變；若x′j大于xavg，則x′j替換為距離該網(wǎng)格最近的反演吸收系數(shù)值小于xavg的網(wǎng)格所對應(yīng)反演的吸收系數(shù)值xc。

對于有地質(zhì)異常揭露的區(qū)域內(nèi)，若x′j>xavg，則x′j并不發(fā)生改變；若x′j

在就近約束條件下，基于MAGA智能算法的層析反演數(shù)學(xué)模型如式(23)所示：

(23)

其中，Q3為反演區(qū)域內(nèi)小于當(dāng)前所有網(wǎng)格對應(yīng)反演的平均吸收系數(shù)值xavg的網(wǎng)格序號集合；Q4為反演區(qū)域內(nèi)大于當(dāng)前所有網(wǎng)格對應(yīng)反演的平均吸收系數(shù)值xavg的網(wǎng)格序號集合。

3.3 平均值約束

在煤礦工作面已揭露區(qū)域范圍內(nèi)，對采用MAGA智能算法進(jìn)行層析反演計(jì)算所得到的每個(gè)網(wǎng)格的吸收系數(shù)值x′j進(jìn)行約束處理。其中，在無地質(zhì)異常揭露的區(qū)域內(nèi)，若小于xavg，則x′j并不改變；若x′j大于xavg，則該網(wǎng)格的吸收系數(shù)值x′j用xavg進(jìn)行替代。

在地質(zhì)異常揭露的區(qū)域內(nèi)，若x′j大于xavg，則x′j并不發(fā)生改變；若x′j小于xavg，則該網(wǎng)格的反演吸收系數(shù)值x′j用xavg進(jìn)行替代。

在平均值約束條件下，基于MAGA智能算法的層析反演數(shù)學(xué)模型為

(24)

將上述3種約束模型融入到MAGA智能算法的層析反演進(jìn)化過程中，以控制其進(jìn)化方向，進(jìn)一步提高工作面內(nèi)地質(zhì)異常層析反演的精準(zhǔn)性。其中已揭露的多源先驗(yàn)約束數(shù)據(jù)主要有巷道揭露數(shù)據(jù)、鉆探揭露數(shù)據(jù)以及回采過程中揭露數(shù)據(jù)。

4 數(shù)值模擬

為了進(jìn)一步驗(yàn)證基于多源數(shù)據(jù)約束的MAGA智能算法在煤礦工作面內(nèi)部地質(zhì)異常層析反演應(yīng)用中的性能和效果，根據(jù)煤礦工作面狹長的特點(diǎn)，建立長600 m和寬200 m的工作面模型，如圖5所示。在該工作面模型中存在2處常見的地質(zhì)異常區(qū)域，由于通過電磁波CT技術(shù)對工作面地質(zhì)異常進(jìn)行反演計(jì)算的基礎(chǔ)是煤層與地質(zhì)異常介質(zhì)對電磁波吸收能力存在一定的相對差異，同時(shí)為了更加直觀地驗(yàn)證本文提出預(yù)測方法的可行性，將異常區(qū)域的電磁波吸收系數(shù)設(shè)置為β1=1.0 dB/m，其余正常區(qū)域的電磁波吸收系數(shù)設(shè)置為β2=0.5 dB/m。

以電磁波CT探測方法為基礎(chǔ)，利用MAGA智能算法并分別結(jié)合范圍約束、就近約束和平均值約束3種約束模型，對該工作面模型內(nèi)的異常區(qū)域進(jìn)行層析反演，最終通過工作面模型區(qū)域內(nèi)的電磁波吸收系數(shù)擬合分布結(jié)果，進(jìn)而推斷出異常區(qū)域的分布情況。

圖5 模型構(gòu)造分布示意Fig.5 Distribution diagram of model construction

在觀測系統(tǒng)布置中，發(fā)射點(diǎn)之間的間距為30 m，對應(yīng)接收點(diǎn)的間距為10 m，射線總數(shù)量N=420條。在層析反演計(jì)算之前，首先對探測區(qū)域進(jìn)行網(wǎng)格離散化處理，每個(gè)網(wǎng)格的尺寸為10 m×10 m，探測區(qū)域共計(jì)劃分成B=1 200個(gè)網(wǎng)格，網(wǎng)格的吸收系數(shù)值即為需要反演求解的目標(biāo)。其中，MAGA智能算法在反演過程中參數(shù)設(shè)置見表3。

將巷道、鉆探以及回采過程中揭露的多源地質(zhì)數(shù)據(jù)通過建立的約束模型融入到MAGA智能算法反演計(jì)算過程中，其中多源先驗(yàn)約束信息分布如圖6所示。其中，鉆探間距為50 m，正常鉆孔深度為150 m，兩巷接替進(jìn)行，對于揭露構(gòu)造的區(qū)域增加鉆孔密度，如圖6(b)所示。巷道揭露數(shù)據(jù)主要體現(xiàn)在巷道兩幫及開切眼周邊的網(wǎng)格吸收系數(shù)值；鉆探及回采揭露數(shù)據(jù)主要體現(xiàn)在所遍歷的網(wǎng)格吸收系數(shù)值。

表3 工作面層析反演模型多種群自適應(yīng)遺傳算法參數(shù)設(shè)置

經(jīng)過500代的進(jìn)化搜索，在MAGA智能算法的基礎(chǔ)上，通過3種不同約束模型分別得到的目標(biāo)函數(shù)f的進(jìn)化過程如圖7所示，最終的收斂結(jié)果見表4。

圖7 MAGA智能算法在不同約束模型條件下隨揭露先驗(yàn)數(shù)據(jù)增加反演得到的目標(biāo)函數(shù)f進(jìn)化曲線Fig.7 Evolution curves of objective function f obtained by MAGA based on different constraint models with the increase of exposed prior data

表4 MAGA智能算法在不同約束模型條件下隨揭露先驗(yàn)數(shù)據(jù) 增加反演得到的目標(biāo)函數(shù)f收斂結(jié)果(Ngen=500)

從上結(jié)果中可以發(fā)現(xiàn)，無約束條件下MAGA智能算法經(jīng)過500代的搜索進(jìn)化最終的目標(biāo)函數(shù)收斂于6.014 dB/m。隨著揭露先驗(yàn)約束數(shù)據(jù)的增加，范圍約束和就近約束模型對約束區(qū)域網(wǎng)格吸收系數(shù)的更新將會產(chǎn)生較大幅度的改變，范圍約束和就近約束模型對應(yīng)反演的最終目標(biāo)函數(shù)值整體呈上升趨勢，并高于無約束條件下MAGA智能算法對應(yīng)的收斂結(jié)果，最終數(shù)值分別為11.08 dB/m(范圍約束)和8.272 dB/m(就近約束)。

而平均值約束模型對應(yīng)反演的最終目標(biāo)函數(shù)值隨著約束數(shù)據(jù)的增加逐漸降低，始終低于無約束條件下MAGA智能算法對應(yīng)的收斂結(jié)果，最終收斂數(shù)值為2.998 dB/m。因此，在平均值約束模型下，MAGA智能算法對地質(zhì)異常反演目標(biāo)函數(shù)最優(yōu)解的搜索性能更好。

經(jīng)過500代的搜索進(jìn)化，基于MAGA智能算法和平均值約束的反演模型隨著外部揭露多源地質(zhì)數(shù)據(jù)的增加，得到的層析成像結(jié)果如圖8所示。

從基于平均值約束的MAGA智能算反演得到的層析成像結(jié)果中可以發(fā)現(xiàn)，隨著外部揭露多源地質(zhì)信息的增加，層析成像結(jié)果中異常區(qū)的范圍更加集中，輪廓也更加明顯，與模型設(shè)置的異常區(qū)范圍也越來越吻合。因此，通過綜合考慮工作面層析反演算法本身的優(yōu)化以及對周圍揭露多源地質(zhì)信息的充分利用，有效提高了工作面內(nèi)地質(zhì)異常層析反演的精準(zhǔn)性，為工作面的透明化建設(shè)奠定了基礎(chǔ)。

圖8 基于平均值約束的MAGA智能算法隨揭露先驗(yàn)約束數(shù)據(jù)增加反演得到的層析成像結(jié)果Fig.8 Tomography results obtained by MAGA based on mean constrained model with the increase of exposed prior data

5 工程應(yīng)用

為了驗(yàn)證所提出方法的可行性，以山西大同礦區(qū)8208工作面探測為例，對工作面內(nèi)的地質(zhì)異常體進(jìn)行反演和預(yù)測。該工作面煤層結(jié)構(gòu)復(fù)雜，厚度不穩(wěn)定，煤厚為4.29～9.44 m，含有1～4層夾矸，夾矸厚度為0.10～0.40 m，煤層大致走向南北，傾向東，傾角為2°～10°。煤層整體賦存形式為背、向斜組成的褶曲構(gòu)造。該工作面的走向長1 000 m，傾向?qū)?50 m，在工作面的兩巷共揭露了7條斷層和3個(gè)陷落柱，其分布如圖9所示。

本次電磁波CT勘探采用WKT-E型坑透儀，發(fā)射點(diǎn)間距為50 m，接收點(diǎn)間距為10 m。觀測射線呈扇形分布，每個(gè)發(fā)射點(diǎn)對應(yīng)的射線覆蓋扇形區(qū)角度為37°，共計(jì)418條射線，如圖10所示。

圖9 8208工作面地質(zhì)概況示意Fig.9 Geological sketch map of 8208 working face

圖10 8208工作面探測觀測系統(tǒng)Fig.10 Layout of the observing systems for 8208 working face

在對8208工作面采集到的電磁波CT數(shù)據(jù)進(jìn)行層析反演之前，首先對工作面探測區(qū)域進(jìn)行網(wǎng)格離散化處理。為了保證每條射線都能分布在網(wǎng)格的內(nèi)部，定義每個(gè)網(wǎng)格的尺寸為10.3 m×10.3 m，網(wǎng)格數(shù)量共計(jì)為1 470個(gè)。

在工作面兩巷和鉆孔揭露的先驗(yàn)數(shù)據(jù)基礎(chǔ)上，基于平均值約束和MAGA智能算法經(jīng)過500代的搜索進(jìn)化，目標(biāo)函數(shù)值最終收斂于2.605 dB/m，如圖11所示，最終得到的層析反演結(jié)果如圖12所示。

由圖12可以看出，共包含2處異常區(qū)，左側(cè)異常區(qū)域主要分布在距離開切眼180～320 m，右側(cè)異常區(qū)主要分布在距離開切眼745～930 m，與周圍揭露構(gòu)造較為吻合。從層析成像反演結(jié)果中可以看出，異常區(qū)的成像效果較為明顯。從左右兩側(cè)異常區(qū)的分布中可初步判斷，2處異常區(qū)附近揭露的地質(zhì)構(gòu)造很大可能貫穿了整個(gè)工作面。因此，工作面在從左往右回采過程中，當(dāng)回采至距離開切眼180 m附近時(shí)，異常構(gòu)造將會對回采工作產(chǎn)生較大的影響，建議在320 m附近重新開切眼，繼續(xù)進(jìn)行工作面的回采，最終建議終采線設(shè)置在745 m附近。

圖11 基于巷道和鉆孔揭露先驗(yàn)數(shù)據(jù)約束的MAGA 智能算法反演得到的目標(biāo)函數(shù)值進(jìn)化曲線Fig.11 Evolution curves obtained by MAGA based on the prior data of roadways and boreholes

圖12 8208工作面層析成像結(jié)果Fig.12 Tomographic result of 8208 working face

8208工作面在回采過程中，從開切眼位置(0 m)回采至185 m位置時(shí)，受構(gòu)造影響，煤層變薄，回采工作被迫停止，并在325 m的位置重新開切眼，與MAGA智能算法層析反演推斷結(jié)果基本一致。

6 結(jié) 論

(1)將工作面內(nèi)地質(zhì)異常反演過程中的矩陣求解問題轉(zhuǎn)換為泛函極值問題，并建立了基于智能算法求解的電磁波反演目標(biāo)函數(shù)模型，為工作面內(nèi)地質(zhì)異常的反演提供了一種新思路。

(2)設(shè)計(jì)了多種群協(xié)同進(jìn)化算子和遺傳參數(shù)自適應(yīng)調(diào)節(jié)算子對SGA進(jìn)行優(yōu)化，并建立了基于多種群自適應(yīng)遺傳算法的地質(zhì)異常反演模型。通過函數(shù)測試發(fā)現(xiàn)MAGA智能算法有效增強(qiáng)了SGA的全局和局部搜索能力，保證了種群的多樣性，提高了目標(biāo)函數(shù)搜索的準(zhǔn)確性。

(3)提出了范圍約束、就近約束和平均值約束3種電磁波層析反演先驗(yàn)數(shù)據(jù)約束模型，建立了基于多源先驗(yàn)數(shù)據(jù)約束的MAGA智能算法反演模型；通過數(shù)值模擬和工程應(yīng)用揭示了平均值約束模型能夠充分利用先驗(yàn)約束數(shù)據(jù)，在反演計(jì)算過程中有效調(diào)整MAGA智能算法的進(jìn)化方向，進(jìn)一步提高了工作面地質(zhì)異常層析反演的精準(zhǔn)性。