杜子龍 儀夢婷 李志遠
(1. 山東科技大學 測繪科學與工程學院, 山東 青島 266590; 2. 中國海洋大學 文學與新聞傳播學院, 山東 青島 266100)
礦區地下開采引起的地表沉陷是一個復雜的四維空間問題,沉陷預測為采空區上方的建筑物,鐵路、水庫等采后支護與治理提供了理論依據,很多研究中注重于開采沉陷的最終結果,即靜態地表移動與變形預測。靜態預測側重于觀測的整體性,對開采沉陷的動態過程研究尚不充分,這對地表建筑物的保護是不利的。實際上,由于礦山開采導致的地表移動與變形是一個動態的過程,變形與移動越劇烈,地表建筑物的損壞就越嚴重[1]。地表的下沉具有三個過程:緩慢發展期、加速發展期、發展衰減期[2]。對地面建筑物等的保護需要考慮變形隨時間發展的全過程,目前常用動態預測地表下沉的時間函數有波蘭學者Knothe[3]于1952年提出的Knothe時間函數;Sroka與 Schober于1982年提出的Sroka-Schober[4]時間函數,Kowalski于1999年提出的廣義時間函數;Gonzalfz等在2007年研究使用的正態分布函數作為時間函數等;我國在采動過程中采用的預計函數多數為基于Knothe時間函數的改進預測模型,大量專家與學者針對Knothe時間函數無法反應下沉速度與加速度的缺陷提出了各種改進方法,常占強等提出的分段Knothe時間函數,提高了預測精度;張兵等改進了分段Knothe函數,使之居于更大的適用范圍,并提高了精度[5-8]。但Knothe時間函數觀測周期長,模型復雜,所需參數多。基于經驗函數的動態預計方法公式計算復雜,且其精度受限于開采實際經驗及巖移參數。基于Logistics模型的下沉預測時間函數能夠利用較少的觀測量提早預測開采沉陷的沉陷量,能夠反映地表下沉的全過程,對于采動過程中引發的地表移動與變形的預測和控制的研究有著積極意義。 Logistic模型能夠滿足單個點動態下沉預測的要求,模型簡單,參數少,應用方便。本文通過控制參與擬合的數據個數,一組使用72期的數據進行擬合,另一組使用81期的數據進行擬合,得到模型預測方程,與實測值的對比表明81期數據組的預計方程要優于72期數據組,通過預測方程對82期91期的沉降值進行預測,得到了較高的預測精度。
在當前的社會學,生物學,計量經濟學等方面,邏輯回歸模型都有廣泛的應用。Logistic模型在最初研究人口增長規律時由馬爾薩斯提出, 之后荷蘭數學家威赫爾斯特將其歸納為一般的數學公式,Logistic曲線在有限空間內的數值的增長會趨于一個穩定值。其發展具有三個過程:緩慢發展期、加速發展期、發展衰減期,它反映了事物的發生、發展、成熟并趨于穩定的過程[9]。
Logistic增長模型的微分方程為:
(1)
對該微分形式進行積分求解可得到該方程的解如下:
(2)
式(2)所表示的曲線即為Logistic模型,總體上呈現S型或反S型曲線(圖1),其特點是單調遞增(S型)或單調遞減(反S型)。由于該曲線為S型,具有地表下沉的相關特征,曲線經歷緩慢發展,加速發展,最終趨于穩定的三個階段,可以將該函數看為Logistic地表下沉時間函數,其中y代表采動過程中該點任意時間的下沉值,參數K為該點最終沉降量W0,參數A跟b為下沉影響因子。其表達式如下:
(3)

圖1 Logistic模型示意圖
Logistic模型曲線形態規律與采動引起的地表下沉規律相符合[10],運用Logistics模型可以預計地表下沉,通過預測的值和實測的值之間的對比來驗證模型。Logistic模型預計流程如圖2所示。

圖2 Logistic模型預計流程
Logistic曲線是以參數K為漸近線,連續的、單調遞減(增)的曲線,曲線隨時間T的變化速度由參數a和b確定。對于式(2)中的三個未知參數,本文將采用在倒數求和法的基礎上,使用非線性回歸的方法,估計模型參數。
(1)倒數求和法確定迭代初始值
為了提高精度,對W0,A,b初值的估計采用倒數求和法[11]。倒數求和法是將樣本均分為三等分,要求時間間隔一樣,設樣本為(ti,Wi),i=1,2,…,n;令r=int(n/3),當樣本數量N不能等分為三部分時,可以適當舍棄開始的1~2個數。
(4)
令D1=S1-S2,D2=S2-S3,則
(5)
(6)
(7)
(2)非線性回歸
為了提高擬合精度,迭代的方法采用Levenberg-Marquardt算法(簡稱L-M算法)。L-M算法是一種以最小二乘原理對多個參數進行擬合的方法,它是介于牛頓法與梯度下降法之間的一種非線性優化方法。L-M算法結合了最速下降法與線性方法,能有效處理冗余參數問題,保證在值域范圍內快速收斂[12]。用L-M算法求得各參數的最佳估計值可以作為Logistic模型參數的最小二乘無偏性估計值。
對某礦區觀測線上的某監測點經過91個周期的地表沉降監測,該監測點的累計沉降量達到38.7 cm。借助Logistic模型進行擬合,將A組使用72期的數據進行擬合并預測沉降,B組使用81期的數據進行模型驗證,并分別預測82期至91期的沉降值,本文采用計算機模擬分析參數。分別得到A、B組的預測函數與實測值的擬合方程及成果圖(圖3)及根據兩組擬合方程得到的82期至91期的沉降值預測曲線圖(圖4)。

圖3 某個地表沉降監測點實測值與A組與B組Logistic擬合結果圖

圖4 某個地表沉降監測點實測值與Logistic模型預測值圖
WA=-0.368 555/(1+14.790 1·
exp(-0.0787 091·t))
WB=-0.381 006/(1+13.702 2·
exp(-0.074 0703·t))
如圖3所示,擬合曲線WA、WB與該點地表下沉實測值在三個階段內的吻合度較高,曲線WB較WA擬合程度更高,原因是WB較WA選取了較多的實測數據,并且具有較多處于下沉第三個階段的實測數據,從預測結果來看,如圖4,WB也比WA更接近真實值,擬合程度高相應的預測精度也高,這說明在本例中,單點下沉符合Logistic下沉時間曲線,兩者具有一致的規律性。
根據統計結果顯示出B組采用的81期數據擬合方程要優于A組采用72期數據擬合方程,其誤差平方和及誤差均方都較A組更接近于0,表明其擬合程度要優于采用較少觀測點的擬合方程。B組參數估計成果表如表1,結果顯示運用L-M算法對三點法估計值進行修正,得到了較好的效果。B組擬合方程統計結果顯示采用該擬合方法的誤差平方和SSE為0.006 621 8,誤差均方MSE為0.000 087 1。通過A、B組Logistic預測方程分別得出82期到91期沉降預測值,并與實測沉降值進行比較,統計誤差,結果見表2。

表1 B組參數估計成果表
(1)表1、表2結果表明,采用81期數據的預測結果要優于采用72天數據的預測結果,兩組的預測差在1 cm左右,結果表明采用該方法時,選用的數據越多,擬合的程度就越高,預測的結果也會更準確。

表2 實測值與預測值對比 單位:m
(2)在采用72期的實測數據進行擬合預測的結果表明,采用較少的數據進行預測時預測值與實測值的誤差為2 cm左右,同樣具有比較高的精度,表明在應用該模型進行預測時,可以采用較少的數據來進行地表下沉預測,選取的地表沉降值數據組將對預測結果產生直接影響。
(3)B組實驗組采用該擬合方法的誤差平方和SSE為0.006 621 8,誤差均方MSE為0.000 087 1,接近于0,擬合精度較高,說明該方法對于動態預測具有一定的精度可靠性。
(4)由表2結果知B組預測函數與實測值的最大殘差出現在第91期,為0.012 067 3 m,最小殘差出現第85期,為0.007 882 1 m,平均誤差僅為0.01 m,具有較高的預測精度。
(1)通過Logistics模型進行參數擬合,運用L-M算法對參數進行修正,統計結果顯示效果較好,并得到了動態預測方程。
(2)從實測數據與擬合曲線圖可以看出,單點下沉關系曲線為反S型。A、B組誤差平方和SSE與均方MSE均接近于0,擬合精度較高。
(3)在大量的實測數據的基礎上,Logistic模型對于單點的地表沉降具有較好的擬合結果,預測精度較高,B組模型預測差最大為0.012 067 3 m,平均誤差為0.01 m。對礦區地表建筑物損壞的實時預防和維護具有實際的指導意義。
(4)對比信息表明采用該方法進行預測時,選用的數據越多,擬合程度越優,沉陷預測結果也更接近于真實值,選取不同的沉降值參與擬合會直接影響預測結果。
(5)Logistics模型簡單,參數少,可以適當減少觀測次數,對地表下沉的動態過程具有一定的使用價值。