開采沉陷中GPS高程擬合代替四等水準測量方法研究

張曉東

(山西三元福達煤業有限公司,山西 長治 046000)

為保障煤礦安全生產,解決好建筑物下、鐵路下和水體下壓煤開采問題,《煤礦測量規程》要求,首采面及重要建構筑物必須建立地表移動觀測站或專門觀測站,觀測站的高程精度須達到四等水準以上精度要求。目前,按進行觀測站的形變監測時所使用的儀器可劃分為傳統手段和新技術。其中,傳統手段主要包括:水準、全站儀和GPS等,新技術主要包括機載LiDAR、InSAR等[1]。由于目前規范中尚未規定新技術應用于地表移動觀測站監測時的具體要求,因此在應用新技術時也會建立觀測站進行驗證。但是,傳統的水準測量和三角高程測量適用于平原和植被不茂密的山區,且效率低、工作強度大;對于植被茂密的山區一般采用GPS觀測,GPS數據處理時多應用卡爾曼濾波、神經網絡等數學模型,這些方法處理的均為某一點的多次GPS數據,且要求數據量大。然而,對于單一工作面觀測而言,視工作面長度的不同,一般僅需觀測8～16期,同一觀測站的觀測次數少,從而導致誤差分布不一定符合正態分布,使得使用常規GPS 數據處理方法時適用性差。因此,針對山區煤炭開采沉陷監測時,應用GPS高程擬合代替四等水準測量的方法研究成為亟待解決的問題。

目前,制約著GPS代替四等水準測量主要有兩個:一是GPS測量的是大地高,水準測量的是正常高,兩者間存在高程異常差值;二是GPS本身的高程觀測精度低于四等水準。針對這兩個主要問題,本文提出了基于觀測線地形拐點的高程異常擬合方法和基于下沉曲線的GPS數據處理方法,實現了在山區開采沉陷中采用GPS高程擬合代替常規四等水準測量。

1 項目概況

研究案例位于山西省長治市武鄉縣,工作面走向長度1 200 m,傾向長度220 m,煤層傾角12°,平均煤厚3.2 m。工作面上方為中低山嶺,溝谷發育,有2條“V”字型沖溝與工作面斜交,坡度較大,最大高差可達190 m,且植被茂密,灌木叢生,現場如圖1所示。

圖1 現場影像Fig.1 On-site images

為研究山區高程異常分布規律,同時驗證GPS處理精度,在工作面上方布設了1條傾向觀測線和半條走向觀測線,相鄰地表移動觀測站間距約為25 m。其中,首次和末次觀測時同時采用三角高程和GPS進行觀測,中間8期數據僅采用GPS進行測量。根據相關規范、文獻和現場監測成果,在山區采用三角高程可達到四等水準的同等精度[2]。

2 基于地形拐點的高程異常擬合

目前,常用的高程異常擬合方法主要有重力法和解析法兩種,其中采用重力法擬合精度較低[3],解析法是應用多面函數、多項式等數學模型擬合高程異常的空間分布,因而對已知點的分布要求較高,而地表移動觀測站一般沿走向和傾向主斷面布置成“十”或“T”型[4],因此將觀測點作為已知點進行全區域擬合是不合理的。針對這一問題,本文提出了一種基于觀測線地形拐點的高程異常擬合方法。

2.1 觀測線上高程異常分布規律

地表移動觀測站一般布設在主斷面上,雖然實際布設時考慮到地形等影響可能略有偏差,但總體呈現為直線或近似直線形狀,因此可將高程異常擬合問題從曲面擬合簡化為曲線擬合,圖2和圖3分別給出了走向和傾向觀測線上高程異常與地形的無因次曲線。

圖2 走向觀測線高程異常與地形的無因次曲線Fig.2 Elevation anomaly and terrain dimentionless curves toward observation line

圖3 傾向觀測線高程異常與地形的無因次曲線Fig.3 Elevation anomaly and terrain dimentionless curves leaning to observation line

從圖中可看出:①高程異常隨地形的變化而變化,整體上具有正相關性;②高程異常的拐點一般出現在地形的拐點位置,但地形的拐點并不全是高程異常的拐點;③在直線型高程異常擬合時,選取地形拐點作為高程異常擬合的關鍵點更合理。

由于傾向觀測線高程異常分布較為簡單,以下僅對走向觀測線高程異常進行擬合。結合地形起伏情況,選取觀測線起點、終點和中間5個地形拐點作為已知點,分別采用三次樣條曲線和多項式曲線擬合整條觀測線的高程異常曲線,其余點用于擬合精度的檢核。

2.2 三次樣條擬合

設區間[a,b]上有n個節點,a=x1

圖4 三次樣條高程異常曲線Fig.4 Elevation anomaly curve of cubic spline

2.3 多項式擬合

多項式擬合是基于最小二乘原理采用n次多項式去擬合包含訓練樣本點在內的所有觀測點[6]。本次樣本點數為7個,最高可進行6次多項式擬合,為對比不同多項式階數的擬合效果,分別進行4次、5次和6次多項式擬合,擬合結果如圖5所示。

圖5 多項式高程異常曲線Fig.5 Polynomial elevation anomaly curves

2.4 擬合精度分析

根據本案例中選取的7個樣本點和29個檢查點,分別進行上述三次樣條擬合和4次、5次、6次多項式擬合,根據擬合結果得出擬合精度如表1所示。

從表1中可知:①擬合精度最高為5次多項式擬合和三次樣條擬合,4次和6次多項式擬合精度最差,說明多項式階次越高并不代表著擬合效果越好;②在實際觀測中,若僅觀測地形拐點,無法確定多項式階次取何值時的擬合效果最好,因此采用三次樣條擬合比采用多項式擬合更合理;③在實際觀測時僅需采用三角高程或水準觀測地形拐點處的高程,極大地提高了觀測效率。

表1 擬合精度統計表Table 1 Fitting accuracy

3 基于下沉剖面函數-EMD的GPS數據處理方法

針對開采沉陷中用GPS 代替四等水準測量時所存在的兩個問題,在上一節研究中指出,可在僅采用三角高程或水準觀測地形拐點處的高程的情況下,采用三次樣條擬合法獲取全觀測線的高程異常值,但仍需解決GPS自身觀測精度低于四等水準的問題。

針對這一問題,本文以地表穩定后走向下沉曲線上36個觀測點的作為研究對象,分別采用常規的EMD、小波分析以及改進的下沉剖面函數3種方法對其GPS監測成果進行處理分析,并根據三角高程的觀測數據對比分析不同方法的處理效果。

3.1 EMD方法

經驗模式分解(empirical mode decomposition,EMD)方法是美國學者HUANG提出的自適應信號分析方法[7]。EMD方法認為,任何復雜非平穩信號都是由簡單的非正弦函數分量信號組成的,一維信號x(t)可分解為多本征函數imf和分解余量r,如公式(1)所示。

將GPS觀測的下沉值看作是由開采引起的下沉信號和GPS觀測的誤差信號構成的序列,圖6給出了下沉值分解的本征函數曲線,由于GPS觀測的下沉曲線主要受采動影響,故將本征函數imfi作為GPS觀測誤差予以剔除,其余項為采動引起。

(1)

圖6 EMD分解圖Fig.6 Decomposition diagram of EMD

圖7顯示了EMD擬合得到的下沉曲線。從圖中可看出:斜率變化小的下沉曲線前半段和后半段擬合效果較好,中間段(400～600 m)偏離觀測值較遠。

圖7 EMD擬合下沉曲線Fig.7 Fitted subsidence curvesof EMD

3.2 小波變換方法

小波變換的過程為:首先將原始信號進行多尺度劃分為高頻和低頻小波系數,并逐層分解低頻小波系數,然后對每層的高頻小波系數進行閾值處理,最后將小波系數進行重構,即濾除噪聲的信號[8]。

將GPS觀測的下沉值看作是由開采引起的下沉信號和GPS觀測的誤差信號構成的序列,GPS誤差可看成是高頻項,此處僅進行一層的小波分解。

1)將GPS實測下沉值通過小波變換分解為a1和b1,如圖8所示。

2)設定閾值,根據誤差傳播定律,h=h末-h首,計算出下沉值的限差56 mm,將其作為閾值。

3)對分解項進行重構,得到下沉曲線的估計值,如圖9所示。

從圖中可看出:經小波變換處理后下沉曲線整體擬合較好,但在下沉曲線以300 m為界的前段和后段,d1項(即誤差項)的大小明顯不同,這與測量誤差的隨機分布不符。

圖8 小波分解圖Fig.8 Wavelet decomposition diagram

圖9 小波變換后下沉曲線Fig.9 Subsidence curves after wavelet transform

3.3 下沉剖面函數-EMD方法

GPS實測下沉曲線可看成是由規律性下沉曲線、非規律性下沉曲線和觀測誤差曲線組成,規律性下沉曲線是滿足下沉剖面函數的曲線,占主要成分;非規律性下沉曲線主要由地形變化、特殊地質構造、煤層傾角變化等因素引起。與非規律性下沉曲線相比,GPS觀測誤差存在于每個觀測點且符號具有隨機性,即具有高頻特征,觀測誤差的絕對值是有限的,經誤差傳播定律推導,本案例中GPS實測中誤差為28 mm。

由于前面EMD和小波變換2種方法是直接對實測下沉曲線進行處理,導致分解的高頻曲線(即誤差曲線)中含有部分非規律性下沉曲線,使得誤差濾除不理想。針對這一問題,本文提出了一種下沉剖面函數-EMD的誤差濾除方法。首先依據剖面函數法和山區下沉盆地主斷面地表下沉預計公式迭代計算出規律性下沉曲線,將殘余曲線(由非規律性下沉曲線和誤差曲線組成)應用EMD方法進行分解,根據2種曲線的特點,分離出誤差曲線。

剖面函數法是以某些函數來表示各種開采條件下主斷面移動變形分布情況,與概率積分法相比,剖面函數法更適用于表示主斷面移動變形情況。目前我國使用最多的剖面函數為負指數函數法,是由煤炭科學研究總院唐山分院根據我國大量實測資料建立的,其下沉表達式如下[9]。

(2)

式中Wm為最大下沉值,m;a為陡度參數;b為形狀參數;c為位置參數;h為平均采深,m。

由于山區地表變形與平原有較大的差異性,根據相關文獻,山區下沉盆地主斷面地表下沉預計公式如下[10]。

W=W1+b2W1tg2α.

(3)

其中,坡度影響系數函數

b2(x)=1/[a(x/r)+t].

(4)

式中,W1為類似地質采礦條件下平地下沉值;α為地面傾角;a、t為待定系數;r為主要影響半徑,m。

公式(2)中W(x)是不受地形影響的下沉值,即公式(3)的W1,將公式(2)帶入公式(3),并將式(3)中的a、t分別替換為d、f,即可得到改進的下沉剖面函數,如公式(5)所示。

(5)

將所有的GPS實測數據應用概率積分法及最小二乘法計算出Wm和r,并作為初始值帶入公式(5),根據收集的工作面資料確定平均采深h和每個觀測點的地面傾角α,帶入公式(5),應用最小二乘法求解待定系數a、b、c、d、f。將擬合結果作為觀測值再次求解Wm和r,直至待定系數沒有明顯變化為止。下沉剖面函數法擬合出的規律性下沉曲線如圖10所示。

圖10 剖面函數法擬合的規律性下沉曲線Fig.10 Regular subsidence curves fitted by section function method

將GPS實測下沉曲線減去規律性下沉曲線,即為殘余曲線(由非規律性下沉曲線和誤差曲線組成),應用EMD方法進行分解,圖11給出了EMD分解的本征函數曲線,由于誤差曲線具有高頻特征,且各觀測點中誤差為28 mm,故將IMF1作為GPS觀測誤差予以剔除,其余項為非規律性下沉曲線。圖12給出了下沉剖面函數-EMD方法的擬合效果,從圖中可看出:經下沉剖面函數-EMD方法處理后下沉曲線與三角高程觀測的下沉曲線擬合較好。

圖11 EMD分解圖Fig.11 Decomposition diagram of EMD

圖12 下沉剖面函數-EMD方法擬合下沉曲線Fig.12 Subsidence curves fitted by subsidence profile function-EMD method

3.4 擬合精度分析

將36個觀測點的三角高程觀測數據作為真值,分別對比GPS觀測數據、三種方法的擬合數據,統計結果如表2所示。

從表2中可看出:

1)GPS直接觀測數據中誤差為21 mm,不能滿足四等水準的精度要求。

2)EMD和小波分析方法能提高GPS觀測數據的精度,但由于規律性下沉占主要成分,分解出的誤差項含有部分非規律性下沉,導致精度提高有限。

3)下沉剖面函數-EMD方法的精度要好于小波變換和EMD方法,根據《工程測量》中山區水準線路閉合差的規定,按3測站觀測一個測點,可推算出測點中誤差最大為10.4 mm,滿足四等水準的要求。

表2 擬合精度統計表Table 2 Fitting accuracy

4 結論

1)通過分析山區高程異常分布規律,提出了基于觀測線地形拐點的高程異常擬合方法,對比發現,三次樣條方法適用性更強,中誤差可將至3 mm,在保證精度的同時減小了三角高程的工作量;

2)根據實際觀測情況,文中提出了基于下沉曲線的下沉剖面函數-EMD的GPS數據處理方法,通過對比EMD、小波變換和下沉剖面函數-EMD3種方法的最終擬合精度可知,下沉剖面函數-EMD方法精度最高,中誤差為8 mm。

3)通過文中提出的方法進行GPS數據處理,研究案例中按每3測站觀測一個測點可推算出測點中誤差最大為10.4 mm,可達到四等水準的精度要求,可為類似山區開采沉陷的工程應用中提供借鑒。