張曉東

(山西三元福達(dá)煤業(yè)有限公司,山西 長(zhǎng)治 046000)

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

目前,制約著GPS代替四等水準(zhǔn)測(cè)量主要有兩個(gè):一是GPS測(cè)量的是大地高,水準(zhǔn)測(cè)量的是正常高,兩者間存在高程異常差值;二是GPS本身的高程觀測(cè)精度低于四等水準(zhǔn)。針對(duì)這兩個(gè)主要問題,本文提出了基于觀測(cè)線地形拐點(diǎn)的高程異常擬合方法和基于下沉曲線的GPS數(shù)據(jù)處理方法,實(shí)現(xiàn)了在山區(qū)開采沉陷中采用GPS高程擬合代替常規(guī)四等水準(zhǔn)測(cè)量。

1 項(xiàng)目概況

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

圖1 現(xiàn)場(chǎng)影像Fig.1 On-site images

為研究山區(qū)高程異常分布規(guī)律,同時(shí)驗(yàn)證GPS處理精度,在工作面上方布設(shè)了1條傾向觀測(cè)線和半條走向觀測(cè)線,相鄰地表移動(dòng)觀測(cè)站間距約為25 m。其中,首次和末次觀測(cè)時(shí)同時(shí)采用三角高程和GPS進(jìn)行觀測(cè),中間8期數(shù)據(jù)僅采用GPS進(jìn)行測(cè)量。根據(jù)相關(guān)規(guī)范、文獻(xiàn)和現(xiàn)場(chǎng)監(jiān)測(cè)成果,在山區(qū)采用三角高程可達(dá)到四等水準(zhǔn)的同等精度[2]。

2 基于地形拐點(diǎn)的高程異常擬合

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

2.1 觀測(cè)線上高程異常分布規(guī)律

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

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

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

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

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

2.2 三次樣條擬合

設(shè)區(qū)間[a,b]上有n個(gè)節(jié)點(diǎn),a=x1

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

2.3 多項(xiàng)式擬合

多項(xiàng)式擬合是基于最小二乘原理采用n次多項(xiàng)式去擬合包含訓(xùn)練樣本點(diǎn)在內(nèi)的所有觀測(cè)點(diǎn)[6]。本次樣本點(diǎn)數(shù)為7個(gè),最高可進(jìn)行6次多項(xiàng)式擬合,為對(duì)比不同多項(xiàng)式階數(shù)的擬合效果,分別進(jìn)行4次、5次和6次多項(xiàng)式擬合,擬合結(jié)果如圖5所示。

圖5 多項(xiàng)式高程異常曲線Fig.5 Polynomial elevation anomaly curves

2.4 擬合精度分析

根據(jù)本案例中選取的7個(gè)樣本點(diǎn)和29個(gè)檢查點(diǎn),分別進(jìn)行上述三次樣條擬合和4次、5次、6次多項(xiàng)式擬合,根據(jù)擬合結(jié)果得出擬合精度如表1所示。

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

表1 擬合精度統(tǒng)計(jì)表Table 1 Fitting accuracy

3 基于下沉剖面函數(shù)-EMD的GPS數(shù)據(jù)處理方法

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

針對(duì)這一問題,本文以地表穩(wěn)定后走向下沉曲線上36個(gè)觀測(cè)點(diǎn)的作為研究對(duì)象,分別采用常規(guī)的EMD、小波分析以及改進(jìn)的下沉剖面函數(shù)3種方法對(duì)其GPS監(jiān)測(cè)成果進(jìn)行處理分析,并根據(jù)三角高程的觀測(cè)數(shù)據(jù)對(duì)比分析不同方法的處理效果。

3.1 EMD方法

經(jīng)驗(yàn)?zāi)Ｊ椒纸?empirical mode decomposition,EMD)方法是美國(guó)學(xué)者HUANG提出的自適應(yīng)信號(hào)分析方法[7]。EMD方法認(rèn)為,任何復(fù)雜非平穩(wěn)信號(hào)都是由簡(jiǎn)單的非正弦函數(shù)分量信號(hào)組成的,一維信號(hào)x(t)可分解為多本征函數(shù)imf和分解余量r,如公式(1)所示。

將GPS觀測(cè)的下沉值看作是由開采引起的下沉信號(hào)和GPS觀測(cè)的誤差信號(hào)構(gòu)成的序列,圖6給出了下沉值分解的本征函數(shù)曲線,由于GPS觀測(cè)的下沉曲線主要受采動(dòng)影響,故將本征函數(shù)imfi作為GPS觀測(cè)誤差予以剔除,其余項(xiàng)為采動(dòng)引起。

(1)

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

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

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

3.2 小波變換方法

小波變換的過(guò)程為:首先將原始信號(hào)進(jìn)行多尺度劃分為高頻和低頻小波系數(shù),并逐層分解低頻小波系數(shù),然后對(duì)每層的高頻小波系數(shù)進(jìn)行閾值處理,最后將小波系數(shù)進(jìn)行重構(gòu),即濾除噪聲的信號(hào)[8]。

將GPS觀測(cè)的下沉值看作是由開采引起的下沉信號(hào)和GPS觀測(cè)的誤差信號(hào)構(gòu)成的序列,GPS誤差可看成是高頻項(xiàng),此處僅進(jìn)行一層的小波分解。

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

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

3)對(duì)分解項(xiàng)進(jìn)行重構(gòu),得到下沉曲線的估計(jì)值,如圖9所示。

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

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

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

3.3 下沉剖面函數(shù)-EMD方法

GPS實(shí)測(cè)下沉曲線可看成是由規(guī)律性下沉曲線、非規(guī)律性下沉曲線和觀測(cè)誤差曲線組成,規(guī)律性下沉曲線是滿足下沉剖面函數(shù)的曲線,占主要成分;非規(guī)律性下沉曲線主要由地形變化、特殊地質(zhì)構(gòu)造、煤層傾角變化等因素引起。與非規(guī)律性下沉曲線相比,GPS觀測(cè)誤差存在于每個(gè)觀測(cè)點(diǎn)且符號(hào)具有隨機(jī)性,即具有高頻特征,觀測(cè)誤差的絕對(duì)值是有限的,經(jīng)誤差傳播定律推導(dǎo),本案例中GPS實(shí)測(cè)中誤差為28 mm。

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

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

(2)

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

由于山區(qū)地表變形與平原有較大的差異性,根據(jù)相關(guān)文獻(xiàn),山區(qū)下沉盆地主斷面地表下沉預(yù)計(jì)公式如下[10]。

W=W1+b2W1tg2α.

(3)

其中,坡度影響系數(shù)函數(shù)

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

(4)

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

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

(5)

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

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

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

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

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

3.4 擬合精度分析

將36個(gè)觀測(cè)點(diǎn)的三角高程觀測(cè)數(shù)據(jù)作為真值,分別對(duì)比GPS觀測(cè)數(shù)據(jù)、三種方法的擬合數(shù)據(jù),統(tǒng)計(jì)結(jié)果如表2所示。

從表2中可看出:

1)GPS直接觀測(cè)數(shù)據(jù)中誤差為21 mm,不能滿足四等水準(zhǔn)的精度要求。

2)EMD和小波分析方法能提高GPS觀測(cè)數(shù)據(jù)的精度,但由于規(guī)律性下沉占主要成分,分解出的誤差項(xiàng)含有部分非規(guī)律性下沉,導(dǎo)致精度提高有限。

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

表2 擬合精度統(tǒng)計(jì)表Table 2 Fitting accuracy

4 結(jié)論

1)通過(guò)分析山區(qū)高程異常分布規(guī)律,提出了基于觀測(cè)線地形拐點(diǎn)的高程異常擬合方法,對(duì)比發(fā)現(xiàn),三次樣條方法適用性更強(qiáng),中誤差可將至3 mm,在保證精度的同時(shí)減小了三角高程的工作量;

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

3)通過(guò)文中提出的方法進(jìn)行GPS數(shù)據(jù)處理,研究案例中按每3測(cè)站觀測(cè)一個(gè)測(cè)點(diǎn)可推算出測(cè)點(diǎn)中誤差最大為10.4 mm,可達(dá)到四等水準(zhǔn)的精度要求,可為類似山區(qū)開采沉陷的工程應(yīng)用中提供借鑒。