基于鉆柱狀態估計的坑道回轉鉆進智能優化方法

姚寧平，魏宏超，張金寶，陸承達，李 浩，姚亞峰，柯友剛，張幼振

(1.煤炭科學研究總院，北京 100013；2.中煤科工西安研究院(集團)有限公司，陜西 西安 710077；3.中國地質大學(武漢) 自動化學院，湖北 武漢 430074)

坑道近水平回轉鉆進是保障煤炭資源安全高效開采的重要手段之一。在鉆進過程中，司鉆調整操作參數時往往依賴自身經驗，難以實現最優鉆進效率同時節省鉆進成本。為平衡鉆進效率、鉆進成本、鉆進安全等多個指標，研究鉆進過程多目標優化，給出操作參數的最優推薦值，是保證鉆進性能的關鍵。

由于鉆進過程充滿復雜性和耦合性，單目標操作參數優化在實際應用中往往存在局限。為平衡鉆進過程的經濟成本，學者們利用機理方法建立了評價鉆速、鉆頭壽命及鉆進效率的多種鉆進多參數多目標優化模型[1-2]；還有學者運用混合蝙蝠算法、支持向量機等數據驅動方法建立了鉆速預測模型[3-4]。基于這些目標評價模型，學者們運用粒子群算法[5]、改進蟻群算法[6]、非支配排序遺傳(NSGA-II)算法[7]等智能優化算法得到了鉆進過程操作參數的最優推薦結果。利用隨鉆測量技術獲知孔底鉆進數據為進行操作參數優化提供了便利[8]。但由于隨鉆測量設備成本較高且多以軌跡參數測量為主，其在坑道回轉鉆進中應用較少。為分析鉆頭的運動狀態和鉆柱系統的穩定性，有學者根據鉆柱系統結構建立了鉆柱動力學模型[9-10]。還有學者將鉆柱動力學模型、目標函數和操作變量相結合，提出了一種基于模型的操作參數優化方法[11]。因此，相比于較高成本的孔底信息測量，建立鉆柱動力學模型進行孔底信息估計，是分析孔底鉆柱運動狀態的一個重要方式。

筆者以坑道近水平回轉鉆進為研究對象，建立了軸向和扭轉維度的鉆柱動力學模型，設計全維狀態觀測器，建立孔口-孔底鉆柱運動狀態映射關系；提出了機械鉆速和鉆頭磨損的優化目標評價模型，基于孔口數據和孔底狀態估計，利用NSGA-II 算法優化轉速和給進壓力這兩個參數；最后利用淮南某煤礦實鉆數據進行驗證與分析，以驗證所提方法的有效性。

1 煤礦井下坑道回轉鉆進過程與特性

在煤礦井下鉆進過程中，坑道鉆機的液壓系統在扭轉維度上通過液壓泵將液壓油輸送至液壓馬達驅動動力頭回轉，進而帶動鉆桿旋轉運動；在軸向維度上通過輸送液壓油至液壓缸，使液壓缸推動鉆桿給進或起拔。鉆頭不斷旋轉和推進，從而實現碎巖鉆進。在實際鉆機操作中，司鉆往往通過遙控器下發控制指令調整動力頭轉速和給進系統給進壓力，以適應不同地質條件和鉆孔要求[12]。

隨著煤礦坑道鉆進不斷深入，地層分布愈加復雜，可能存在褶皺、斷層等地質構造，容易造成卡鉆埋鉆事故，安全鉆進難度大。為了確保煤礦坑道鉆進的安全和效率，需要對鉆機操作參數進行嚴格的監控和調整，以保持煤礦鉆進過程穩步進行。煤礦井下坑道鉆進過程的復雜性、非線性使得鉆進參數受到多種因素的影響，研究鉆進參數多目標優化方法能夠提高鉆進效率和成孔質量。

2 鉆進參數多目標優化

煤礦井下鉆進參數多目標優化通過考慮多個目標，綜合評估不同參數組合的優劣，從而選擇最優的鉆進參數，這樣可以提高鉆進效率、降低成本、減少風險，從而提高煤礦的生產效益和安全性。針對本文的研究問題，設計的技術路線如圖1 所示。

圖1 鉆進參數多目標優化技術路線Fig.1 Roadmap of multi-objective optimization technology for drilling parameters

圖 2 孔內狀態觀測器原理Fig.2 Principle of the in-hole state observer

首先將孔口信息輸入到鉆柱運動模型獲取狀態信息，進而通過非線性狀態觀測器得到孔內狀態估計信息；然后利用獲取的孔內信息，通過非線性多元回歸方法辨識機械鉆速模型的煤巖層系數、門限鉆壓和轉速指數3 個系數。最后利用NSGA-II 算法進行多目標優化求解，得到轉速和鉆壓等操作參數，并對優化的結果進行分析。

2.1 基于孔口數據驅動的鉆柱狀態估計

針對煤礦坑道近水平回轉鉆進中孔內數據難以獲得、操作參數優化缺少依據等問題，進行基于軸向和扭轉鉆柱模型的非線性狀態觀測器設計。

所使用的動力學模型采用集中質量法對鉆柱系統各部分動力學特性進行建模。其中，軸向鉆柱動力學模型采用下式表示：

式中：jal為各集中質量單元的質量，l取1，2，3，下標a為軸向鉆柱模型；zl、、分別為各單元的軸向位移、軸向速度和軸向加速度，l取1，2，3；hal為各單元軸向阻尼，l取1，2，3；kal為各單元軸向剛度，l取1，2；W0為給進力；Wb為軸向鉆頭-煤巖相互作用力。

選擇 [z1,z2,z3]T作為狀態變量可構成模型系數矩陣Ja、Ha、Ka，表示為：

扭轉鉆柱動力學模型采用下式表示：

式中：jtl、htl、ktl分別為鉆柱系統各單元的轉動慣量、扭轉阻尼、扭轉剛度，jtl和htl中l取1，2，3，ktl中l取1，2，下標t 為扭轉鉆柱模型；θl、、分別為鉆柱系統各單元的轉動位移、轉速和扭轉加速度；T0為動力頭輸入扭矩；Tb為鉆頭-煤巖相互作用帶來的負載扭矩。選擇[θ1,θ2,θ3]T作為狀態變量可構成模型系數矩陣Jt、Ht、Kt，表示為：

為后續進行狀態變量觀測，需要重構狀態變量，選擇[,,,z2-z1,z3-z2]T作為軸向鉆柱模型的狀態變量xa，選擇[,,,θ2-θ1,θ3-θ2]T為扭轉鉆柱模型的狀態變量xt。

由于軸向鉆柱模型和扭轉鉆柱模型的狀態空間方程結構相同，結合所重構狀態向量，通過矩陣代換，采用以下通式對軸向和扭轉鉆柱動力學模型的狀態空間方程進行描述：

式中：i=a為軸向鉆柱模型，i=t為扭轉鉆柱模型；xi為鉆柱模型的狀態變量；yi為輸出變量；Ai為狀態轉移矩陣；Bi為輸入矩陣；Ri為鉆頭-煤巖作用輸入矩陣；Ci為輸出矩陣；n為模型的自由度；Ji、Hi、Ki為動力學系數矩陣，由已建立的軸向、扭轉動力學模型系數確定，具體形式可參考文獻[13]；ui為鉆柱模型輸入(給進力或扭矩)；wi(xi)為鉆柱模型鉆頭-煤巖作用，||wi(xi)||

后續對狀態觀測器進行李雅普諾夫穩定性分析時，會用到以下幾個引理。

引理1：當V(t)是連續函數，w(t)為指數有界函數時，如果存在α>0 和b>0，那么對于所有t>t0有以下不等式成立：

式中：t為時間；t0為初始時間；b為標量。

進一步可知：

引理2：當系統具有指數衰減率α>0 并存在指數最終有界時，若有常數β≥1，則：

式中：κ為||x(t)||2的最終上界；η2為給定標量，代表θ 所存在的范圍大小。

狀態空間方程的優點是能夠更直觀地描述系統的行為，也更容易進行分析和控制。它可以用來描述線性和非線性系統，也可以用來描述時變和時不變系統。狀態空間方程在控制工程、信號處理、機器人控制、電力系統等領域都有廣泛的應用。

狀態觀測器是一種用于估計系統狀態的工具。它的作用是通過測量系統的輸出和輸入，以及系統的動態模型，來估計系統的狀態變量，包括未被直接測量的狀態變量?；诳卓谝阎畔?，通過狀態觀測器估計孔內信息，狀態觀測器原理如圖2 所示，圖中s表示積分。

式(5)中引入反饋矩陣L，并加入觀測器估計值與系統的被估計量的誤差yi-，構造如下狀態觀測器：

基于建立的模型狀態空間方程，輸入給進力或扭矩，利用可測量參數(動力頭轉速和動力頭鉆速)，通過觀測器估計煤礦井下孔底所需參數(鉆頭轉速和鉆頭鉆速)。

引入狀態估計誤差量ri=xi-來分析觀測器的收斂性及收斂速度，式(5)與式(11)相減可得：

式(12)中引入狀態估計誤差量ri可得：

已知||wi(xi)||

針對的狀態估計誤差ri，構造李雅普諾夫函數V：

式中：P為正定對稱矩陣。對式(14)求導，并在等式兩邊加上2αV可得：

其中，α>0且為實數。

由上述條件可知，存在φi(xi)Tεφi(xi)<，其中ε為設置大于0 的中間變量，后續可根據線性矩陣不等式計算得出，代入式(15)可得：

設σ(t)=[ri,φi(xi)]T，基于線性矩陣不等式性質設置反饋增益矩陣L=，μ為代求常數，代入式(16)可得：

式 中：λmin(P) 為P的最小特征值；為一個不小于λmax(P)的常數，由式(18)可得：

將給進力和扭轉分別作為軸向鉆柱和扭轉系統模型的輸入，通過觀測器式(11)獲取狀態變量的值，即可獲取鉆頭鉆速估計，為下一步多目標優化問題提供孔內信息。

2.2 鉆進參數優化方法

由于煤礦環境的復雜性和多樣性，在實際煤礦鉆進工程中，需要把如機械鉆速、鉆頭磨損等目標綜合考慮。傳統的多目標優化問題是轉化成單目標優化問題來解決，這種方法受人為影響，優化結果具有較大誤差。隨著智能優化算法的發展，國內外學者將這些算法應用到鉆進參數優化中，本文嘗試用NSGA-II 算法來解決多目標優化問題。

在煤礦井下實際鉆進中，鉆進速度和鉆頭磨損對煤礦生產效率、安全和環保等方面都有著重要的影響，希望機械鉆速和鉆頭磨損都能達到最優的指標。以最大鉆速和最小鉆頭磨損為優化目標，建立能定量評估鉆進效率同時能反映鉆進參數對鉆進過程影響的多目標優化模型。修正的楊格機械鉆速模型[14]是一種典型的物理模型，綜合考慮了多種因素對鉆速的影響：

式中：v為機械鉆速，m/h；kd為地層可鉆性系數；Cp為壓差影響系數；CH為水力凈化系數；p1為鉆壓，MPa；p2為門限鉆壓，MPa；N為轉速，r/min；λ為轉速指數；C2為切削齒磨損系數；h為鉆頭切削齒相對磨損量。其中CP和CH的值在本文中設定為1，kd、p2、λ的值與實際鉆進現場工況有關，通過回歸分析得到這些值。為便于回歸分析求取模型參數，將鉆速模型進行簡化：

式中：A1=kd為煤巖層系數；A2=p1為門限鉆壓；A3=λ為轉速指數。通過回歸分析可得到A1、A2和A3。其中，鉆速模型中的1/(1+C2h)是切削齒磨損對鉆速的影響，該影響在下文的鉆頭磨損模型中已經體現，故在鉆速模型中將這部分忽略。該模型是在楊格鉆速模型基礎上，綜合考慮鉆壓、轉速、鉆井液性能以及水力因素對鉆速的影響，并結合鮑格因等人的研究成果歸納建立[15]。

本文中的鉆頭磨損主要考慮切削齒磨損，切削齒磨損的速度模型如下：

式中：D1、D2為鉆壓影響系數；Q1、Q2為轉速影響系數；Af為地層研磨性系數；C1為切削齒磨損減慢系數。相關系數的具體值通過查閱相關資料在下文的實驗部分給出。

鉆頭轉速和鉆壓是影響鉆進效率和質量的兩個重要因素，同時也是鉆進過程中易于控制的兩個變量。以鉆頭轉速N和鉆壓W為決策變量，綜合考慮機械鉆速和鉆頭磨損，建立多目標優化模型，優化問題可用下式的數學模型表示：

回歸分析是用于研究自變量和因變量之間關系的一種統計學方法，利用孔口數據結合狀態觀測器估計的孔底信息，可以通過回歸分析求得機械鉆速模型簡化形式(式(21))中的系數A1、A2和A3。具體的步驟如下。

(1) 收集數據：收集的數據包括狀態觀測器觀測的孔底鉆速，孔口轉速和鉆壓。

(2) 數據處理：使用均值濾波對初始數據進行預處理。

(3) 估計參數：利用非線性多元回歸估計鉆速模型中的參數。

(4) 檢驗模型：對模型的擬合程度進行判斷，檢驗回歸模型是否合適，包括檢驗模型的顯著性和擬合程度等。

上述優化問題(式(23))提出了兩個優化目標，需要利用有效的算法尋得最優解。NSGA-II 算法在解決雙目標優化問題時表現出良好的效果，已廣泛應用于實際工程問題中，該方法適用于本文的研究。

NSGA-II 主要克服了NSGA 算法的一些缺點。在NSGA-II 中，使用了快速非支配排序算法，針對多個優化目標無法同時達到最優，將種群根據相應的支配關系劃分為不同級別的Pareto 解，算法的目標是找到一組盡可能接近Pareto 最優域的解，即對算法的收斂性有要求，同時要保證這組解的多樣性。NSGA-II同時采用精英保留策略，父代優秀的個體被保留下來的概率大大增加。NSGA-II 算法采用擁擠度比較算子作為種群中個體間的比較標準，克服了NSGA 算法需要人為指定共享參數的缺陷。算法流程如下。

(1)初始化種群：隨機生成一定數量的個體作為初始種群，這里的個體是一組機械鉆速和鉆頭磨損的值。

(2)計算適應度：對于每個個體，計算其在多個目標函數下的適應度值。

(3)非支配排序：根據個體之間的支配關系，將種群中的個體劃分為不同的等級，即非支配解、一級支配解、二級支配解等。

(4)計算擁擠度：對于每個等級中的個體，計算其在目標空間中的擁擠度，即其周圍個體的密度。

(5)選擇操作：根據非支配排序和擁擠度計算，選擇一定數量的個體作為下一代種群。

(6)進化操作：通過交叉和變異等進化操作，生成新的個體，加入下一代種群中。

(7)終止條件判斷：判斷是否滿足終止條件，如達到最大迭代次數或目標函數值達到一定精度等。

(8)輸出結果：輸出最終的帕累托前沿解集。

其中種群大小、迭代次數、交叉和變異的概率，根據實際問題進行設置。從得到的Pareto 解中，根據現場的實際情況和鉆進需求，挑選出滿意解作為鉆進參數的優化值，為決策者在實際工程問題中提供有效參考。

3 仿真分析與實驗

為了驗證上述方法的有效性和可行性，基于狀態觀測器所獲得的數據進行仿真實驗，對相關鉆進參數進行優化，上述算法基于Python 3.8 實現。

3.1 回歸分析與優化結果

通過非線性多元回歸對機械鉆速模型式(18)進行擬合，所得的3 個參數值為A1=12.939 2，A2=8.447 9，A3=-0.327 7。通過回歸分析常用評價指標對結果進行分析見表1。

表1 回歸分析評價參數Table 1 Regression analysis evaluation parameters

均方誤差是觀測值與真值偏差的平方和與觀測次數的比值，均方誤差的值越小，說明預測模型描述實驗數據具有更好的精確度。平均絕對誤差能更好地反映預測值誤差的實際情況。R2介于0 到1 之間，可以理解為擬合程度的好壞，R2為0.988 1，可以認為模型能解釋98%的輸出變化。

通過非線性多元回歸確定鉆速模型的系數后，利用機械鉆速模型得出的鉆速和狀態觀測器估計的孔底鉆速值對比如圖3 所示。機械鉆速模型計算的鉆速和狀態觀測器估計的鉆速相差較小，可認為非線性多元回歸擬合的鉆速方程效果較好，具有實用價值。

圖3 計算鉆速與孔底估計鉆速對比Fig.3 Comparison between ROP calculated by ROP model and downhole ROP estimated by state observer

鉆頭磨損方程中的相關系數可根據實際工況和所使用的鉆頭情況來確定[16]。D1=0.014 8，D2=6.38，Q1=1.5，Q2=6.53×10-5，C1=5，Af=0.002 28。優化算法中的相關參數設置根據經驗來確定，種群大小為100，進化代數為120，圖4-圖6 表示不同遺傳代數解的分布情況。

圖4 未優化時初始解的分布情況Fig.4 Distribution of initial solution without optimization

圖5 循環次數為60 時解的分布情況Fig.5 Distribution of solution under 60 cycles

圖6 優化完成時解的分布情況Fig.6 Distribution of solutions under completed optimization

從圖4-圖6 可以看出，初始時，種群散亂地分布在解的空間中，當種群進化到60 代時，種群個體已經分布得較為均勻，當完成最終的進化時，種群個體均勻地分布在解空間中，獲得的Pareto 曲線較為平滑，最終解的分布情況證實了NSGA-II 算法的有效性，選取3 組Pareto 最優解分析對比，見表2。

表2 3 組Pareto 最優解對比Table 2 Comparison of three sets of Pareto optimal solutions

由表2 可知，當鉆速增大時，鉆頭磨損速度加快，可見兩者是相互制約的關系，一個目標性能的改善將引起另一個目標性能的降低。在實際煤礦鉆進的現場中，對鉆速有一定的要求，首先考慮滿足鉆速在一定的范圍內，然后去考慮鉆頭磨損速度，滿意解的選取還是需要根據現場的實際工況來確定。

3.2 實測數據實驗分析

為了進一步說明鉆柱狀態估計對優化結果的提升，利用淮南某煤礦實測數據進行了實驗分析。所使用的數據為淮南顧北煤礦采集到的井下鉆進數據，此時鉆進深度為18.75 m，所使用的鉆桿長度為0.75 m。在下文的實驗中，鉆頭磨損模型相關系數與上文中進行仿真數據使用一致，機械鉆速模型中的參數通過實測數據進行回歸擬合得到。文獻[17]僅利用孔口的實測數據，未使用狀態觀測器獲取孔底信息，并利用NSGAII 算法進行仿真優化。本文方法利用狀態觀測器估計的孔底信息，并使用NSGA-II 算法進行優化。將兩種方法得到的結果同時與現場實際操作參數進行對比，并計算預計鉆速，實驗結果見表2 和表3。由于機械鉆速的增大，另一個優化目標鉆頭磨損速度也會增大，表2 中優化后的目標值選擇方案是保持鉆頭磨損速度變化較小，讓鉆速較大的增大，對比情況見表3。

表3 實際煤礦鉆進數據與優化后的數據對比Table 3 Comparison between actual coal mine drilling data and optimized data

分析表3 可得出如下結論：使用文獻[17]中的方法，基于實測孔口數據進行優化，優化后的鉆速為22.64 mm/s，相比于實際采集數據中的鉆速19.68 mm/s，提升了15.04%；基于鉆柱狀態估計孔底信息進行優化后，鉆速為26.07 mm/s，相比于實際采集的數據鉆速提升了32.47%。表明構造狀態觀測器估計孔底信息，能讓優化結果提升更高，證實了本文所提方法的有效性。

4 結論

a.分析了煤礦井下鉆進過程存在的優化問題，綜合考慮機械鉆速模型和鉆頭磨損模型這兩個優化目標，通過NSGA-II 算法求解優化問題，實現了最優鉆壓和轉速的匹配，用實測數據驗證了該方法的有效性。

b.通過建立鉆柱動力學模型和設計全維狀態觀測器，建立了孔口-孔底鉆柱運動狀態映射關系，實現了孔底不可測量數據的估計。實驗結果表明，基于鉆柱狀態估計孔底信息進行優化后的鉆速優于僅利用孔口的實測數據優化的鉆速。

c.提出的鉆速模型計算鉆速與狀態觀測器估計孔底鉆速擬合程度好，鉆柱狀態估計對優化結果的提升顯著，但還需要通過實踐驗證優化結果的準確性。