田 雨，雷少剛，卞正富

(1.中國礦業(yè)大學(xué) 國土資源研究所，江蘇 徐州 221116; 2.礦山生態(tài)修復(fù)教育部工程研究中心，江蘇 徐州 221116)

排土場是露天礦存放廢棄石材土料的地方，是露天礦整個生產(chǎn)周期中不可或缺的一項永久性工程建筑，排土場可分為內(nèi)排土場和外排土場，其中，內(nèi)排土場是采礦過程中將剝離的巖土直接排棄到采空區(qū)進而形成的堆積體。內(nèi)排土技術(shù)運距短，可以極大地降低排土成本，相比于外排土場壓占面積更小，同時內(nèi)排土場經(jīng)平整恢復(fù)，可為復(fù)墾和工程建設(shè)提供可觀的土地，成為礦區(qū)實現(xiàn)可持續(xù)發(fā)展的關(guān)鍵區(qū)域之一。但排土場因其填方土石的自重作用、降雨徑流等因素影響，堆積成形后隨時間推移不可避免的發(fā)生變形。目前學(xué)者們對于排土場的研究主要集中在邊坡穩(wěn)定性上[1-3]，對于排土場主體的變形沉降研究較少，而內(nèi)排土場土層沉降壓實過程中，土壤質(zhì)地和滲透性發(fā)生變化，有效水分和養(yǎng)分隨之改變，不利于植被及微生物生長[4-5]，工程建設(shè)更是對其主體變形的速率、累積量等都提出了嚴(yán)格要求，本文旨在通過時序SAR數(shù)據(jù)高精度監(jiān)測內(nèi)排土場沉降，探究其沉降規(guī)律，并依據(jù)監(jiān)測數(shù)據(jù)采用混沌時間序列相空間重構(gòu)理論結(jié)合二階Volterra自適應(yīng)濾波建立預(yù)測模型，進行內(nèi)排土場沉降動態(tài)預(yù)計，為內(nèi)排土場的充填復(fù)墾[6]，工程利用等提供一定參考。

1 研究區(qū)概況

勝利露天煤礦隸屬于神華北電勝利能源有限公司，其中勝利一號露天礦位于勝利煤田的中西部，錫林浩特市北郊5 km，本煤田位于內(nèi)蒙古高原中部，氣候干燥，丘陵地形，地形標(biāo)高970～1 326 m，露天礦共有4個排土場，其中3個外排土場:南排土場、北排土場、沿幫排土場。礦區(qū)至2013年實現(xiàn)全部內(nèi)排，內(nèi)排土場覆蓋范圍約2.3 km×1.1 km，其上植被較少，適宜利用SAR數(shù)據(jù)進行沉降監(jiān)測，如圖1所示。

圖1 研究區(qū)示意Fig.1 Schematic map of the study area

2 沉降監(jiān)測

2.1 SBAS算法

小基線集(SBAS)技術(shù)由BERARDINO等于2002年提出[7]，其基本原理是將獲得的同一個區(qū)域的SAR數(shù)據(jù)按照空間基線組成多個集合，利用最小二乘方法分別得到每個集合的沉降時間序列，然后利用奇異值分解法(SVD)將多個集合聯(lián)合求解從而得到整個觀測時間周期的沉降時間序列[8]，經(jīng)眾多國內(nèi)外學(xué)者驗證，它可以有效抑制由于時間、空間造成的失相干和大氣效應(yīng)，尤其在低植被覆蓋區(qū)域可得到高精度的沉降時序數(shù)據(jù)。

2.2 數(shù)據(jù)準(zhǔn)備

Sentine-1 A于2014-04-03發(fā)射，搭載了C波段SAR傳感器，時間分辨率為12 d，共有4種數(shù)據(jù)獲取模式:SM(Stripmap Model)，IW(Interferometric Wide Swath)，EW(Extra-Wide Swath)和WM(Wave Model)[9]。本文選用Sentine-1 A IW模式下的單視復(fù)數(shù)據(jù)(SLC)，極化方式為VV和VH，其空間分辨率為5 m×20 m，數(shù)據(jù)時間跨度為2016-09-29—2018-02-03，除缺失2017-06-20一期數(shù)據(jù)外，其余數(shù)據(jù)時間間隔均為12 d，共計41期。

本文選取時間基線閾值為60 d，時空基線如圖2所示。

圖2 干涉圖的時空基線Fig.2 Spatial-temporal baselines of the interfergrams

2.3 沉降結(jié)果及分析

通過生成的干涉圖，經(jīng)相位解纏、軌道精煉和重去平、時間高通濾波、空間低通濾波等得到只含形變信息的差分干涉圖，進而SVD求解得到時間形變序列[10]。

圖3為處理獲得的2016-09-29—2018-02-03地面累積沉降量，由圖3可以看出，勝利一號露天礦內(nèi)排土場的累積沉降量和沉降面積隨時間推移不斷增大，總體呈現(xiàn)西部累積沉降大，東部累積沉降小的規(guī)律，截至2018-02-03(參考2016-09-29)，最大累積沉降量達394 mm。

圖3 沉降量Fig.3 Maps of settlement

為進一步研究內(nèi)排土場沉降的時空分布規(guī)律，將2018-02-03累積沉降量按照50 mm等間隔劃分為Ⅰ～Ⅷ共8個區(qū)域，并從累積沉降量最大區(qū)域向外延伸兩條直線A和B，每個區(qū)域內(nèi)均在直線上隨機取一點，如圖4所示。

圖4 線和點位置Fig.4 Location of lines and points

分別統(tǒng)計沿直線A,B的累積沉降量和A1～A8點，B1～B8點處的沉降時間序列。

圖6 A1～A8沉降時間序列Fig.6 Settlement time series of A1-A8

圖7 B1～B8沉降時間序列Fig.7 Settlement time series of B1-B8

圖5為觀測時間周期內(nèi)累積沉降量沿直線A和B的剖面圖，可以看出A,B兩條曲線相似度很高，內(nèi)排土場累積沉降量空間上呈現(xiàn)明顯的半漏斗狀，總體上累積沉降量與到礦坑距離成反比，符合內(nèi)排土場排棄規(guī)律。圖6，7分別代表直線A和B在每個區(qū)域取的隨機點隨時間的沉降規(guī)律，可看出同一個區(qū)域內(nèi)兩點的累積沉降量-日期曲線近乎相同，兩圖均反映出在Ⅷ區(qū)域，沉降速率最后趨近于0，表明已達穩(wěn)定狀態(tài)，可初步判斷Ⅷ區(qū)域已基本滿足了土地復(fù)墾及建設(shè)簡單構(gòu)筑物的基本要求，Ⅲ～Ⅶ區(qū)域隨時間推移沉降速率有所減緩，但仍未穩(wěn)定，還需一定時間的固結(jié)壓實才可以考慮工程利用，Ⅰ,Ⅱ區(qū)域觀測周期內(nèi)沉降速率大且基本保持不變，沉降劇烈，處于沉降活躍期，主要是由于其為觀測周期起始時覆土，比其他區(qū)域較晚，且臨近礦坑，受施工、雨水沖刷等影響邊坡部分土壤松動滑落所致，存在滑坡、泥石流風(fēng)險，是后期沉降監(jiān)測的重點區(qū)域。

2.4 下沉系數(shù)估計

礦區(qū)傳統(tǒng)的下沉系數(shù)是井工礦開采沉陷和建筑物下采煤進行地表移動變形預(yù)計時的一個重要參數(shù)。但露天礦內(nèi)排土場覆土工后物料也會在自重等作用下逐漸壓實沉降，期間地表難免發(fā)生移動變形，而內(nèi)排土場的再利用如工程建設(shè)和土地復(fù)墾等都對土體穩(wěn)定性提出了嚴(yán)格要求，因此露天礦內(nèi)排土場主體穩(wěn)定性必須考慮下沉系數(shù)，本文定義內(nèi)排土場下沉系數(shù)為排土工后最終沉降量與初始覆土高度的比值:

(1)

式中，q為內(nèi)排土場下沉系數(shù);H0為內(nèi)排土場原始覆土高度，m;Hcm為內(nèi)排土場最終沉降量，cm。

依據(jù)內(nèi)排土場排土?xí)r序圖，得A1～A8，B1～B8覆土?xí)r間見表1。

表1A1～A8，B1～B8覆土日期
Table 1 Covering soil times ofA1-A8andB1-B8

點號覆土日期(年-月)點號覆土日期(年-月)A12016-10B12016-09A22016-08B22016-06A32016-05B32016-05A42016-04B42016-04A52016-03B52016-03A62015-03B62015-03A72014-09B72013-12A82013-07B82013-07

因內(nèi)排土場覆土層力學(xué)性質(zhì)相近，覆土高度幾乎相同，所以可近似認(rèn)為A1～A8，B1～B8為內(nèi)排土場覆土工后不同時期的沉降狀態(tài)，其中A1接近為監(jiān)測周期開始時最新的覆土點，選擇A1～A8累積沉降量依據(jù)時間關(guān)系進行曲線擬合，估計A1點累積下沉量，結(jié)果如圖8所示。

圖8 A1點累積下沉擬合曲線Fig.8 Cumulative settlement fitting curve of A1

由圖8可知，累積下沉量隨時間呈現(xiàn)指數(shù)分布，R2高達0.991，A1最終下沉量約為662 mm，下沉量達到相對穩(wěn)定約需要110期監(jiān)測，即1 320 d。A1點排土場原始覆土高度為103.58 m，因此可求得觀測周期內(nèi)，內(nèi)排土場下沉系數(shù)約為0.639 cm/m。值得注意的是，A1點沉降監(jiān)測時間受限于原始數(shù)據(jù)，略晚于真實覆土?xí)r間，而覆土初期是沉降劇烈期，因此本文求出的下沉系數(shù)小于實際下沉系數(shù)，只能作為后續(xù)研究的參考。

3 沉降預(yù)測

3.1 混沌系統(tǒng)的識別

混沌運動可認(rèn)為是確定性非線性動力系統(tǒng)中獨特的重復(fù)性運動，具有不同于一般確定性運動的幾何和統(tǒng)計特征[11]。混沌系統(tǒng)與高維隨機運動的主要區(qū)別在于:混沌系統(tǒng)通過相位重構(gòu)，其軌跡可預(yù)測，而隨機運動則不能[12]。

時間序列可分為周期性序列、準(zhǔn)周期性序列和混沌序列[13]。對于指定的時間序列需要進行性質(zhì)鑒別?；煦鐣r間序列的判定方法有很多，目前為止使用最為廣泛的Lyapunov指數(shù)法[14-15]。Lyapunov指數(shù)反映了混沌動力學(xué)系統(tǒng)對初始條件的敏感性，最大的Lyapunov指數(shù)是否大于0可作為判定動力學(xué)系統(tǒng)混沌性的一個依據(jù)?，F(xiàn)有的Lyapunov指數(shù)計算方法有定義法、Wolf方法、P-范數(shù)方法、小數(shù)據(jù)量法以及奇異值分解方法等[14-15]。

3.2 相空間重構(gòu)

根據(jù)荷蘭數(shù)學(xué)家TANKENS[16]提出的延遲嵌入定理，即只要將采集到的時間序列嵌入到延遲時間為τ的m維相空間中，則所重構(gòu)的m維狀態(tài)空間與原混沌動力系統(tǒng)的幾何特征等價，它們具有相同的拓?fù)浣Y(jié)構(gòu)。因此，選取合適的嵌入維數(shù)m和延遲時間τ是相空間重構(gòu)的關(guān)鍵步驟。

本文選取互信息量法求延遲時間τ。Swinny和Frasert[17]提出的互信息法是一種通過度量2個變量之間的隨機關(guān)聯(lián)程度來判斷系統(tǒng)非線性相關(guān)性的方法，計算較為準(zhǔn)確，在相空間重構(gòu)延遲時間的確定中應(yīng)用較為廣泛。

對于時間序列{x(t1)，x(t2)，…，x(tn)}，當(dāng)延遲時間為τ時，序列變?yōu)閧xi+τ，i=1，2，…，n}，則互信息函數(shù)為

(2)

其中，P(xk)為xk在時間序列{xi+τ，i=1，2，…，n}中出現(xiàn)的概率;P(xkτ)為xk+τ出現(xiàn)的概率；P(xk，xk+τ)為xk和xk+τ共同出現(xiàn)的概率，在使用互信息量法計算最佳時間延遲時一般采用互信息函數(shù)I(τ)的第1個極小值點作為最優(yōu)時間延遲。

選取A1,B3,A5,B7各點的相鄰兩期累積沉降量差值構(gòu)成4組長度為40的時間序列，因缺失一景SAR數(shù)據(jù)無法組成等時間間隔序列，經(jīng)研究采用3次樣條插值對數(shù)據(jù)進行加密，從而得到4組以1 d為間隔，長度為492的時間序列數(shù)據(jù)。圖9是4組時間序列數(shù)據(jù)的互信息，由圖9可以看出A1,B3,A5,B7的時間延遲分別為4,3，3,3。

圖9 互信息法選擇時延Fig.9 Determine delay time through mutual information method

對于嵌入維數(shù)m的選取利用Cao Liangyue[18]提出的改進虛假臨近點法，即Cao方法，定義:

(3)

定義E(m)為所有i對應(yīng)a2(i，m)的均值，即

(4)

為了進一步研究E(m)隨m的變化情況，定義

(5)

若時間序列存在混沌吸引子，則當(dāng)嵌入維m大于某個值m0時，E1(m)逐漸穩(wěn)定，那么m0+1可以作為系統(tǒng)的最佳嵌入維數(shù)。但在實際生活中，時間序列數(shù)據(jù)長度有限，有時難以準(zhǔn)確判斷E1(m)為否已達到收斂狀態(tài)，因此，額外定義

(6)

(7)

若E2(m)的值恒為1，則表明所選時間序列為隨機時間序列，數(shù)據(jù)間無相關(guān)性;若其不恒為1，則表明時間序列數(shù)據(jù)間的相關(guān)性受嵌入維數(shù)m影響[19]，因此可以依據(jù)E1(m)和E2(m)的值共同確定出最佳嵌入維數(shù)。

圖10為4組時間序列數(shù)據(jù)的E1(m)和E2(m)值，由圖10可以確定A1,B3,A5,B7的最佳嵌入維數(shù)分別為6,5，6,6。

圖10 Cao方法選擇嵌入維Fig.10 Determine embedding dimension through Cao method

得到A1,B3,A5,B7四組時間序列的最佳延遲時間τ和嵌入維數(shù)m后，采用小數(shù)據(jù)量法分別計算其最大Lyapunov指數(shù)，結(jié)果見表2。由表2可以看出，4組時間序列的最大Lyapunov指數(shù)均大于0，所以可以判斷A1,B3,A5,B7四組時間序列均具有混沌特征。

表2 最大Lyapunov指數(shù)
Table 2 Largest Lyapunov exponent

時間序列A1B3A5B7最大Lyapunov指數(shù)0.020 90.024 10.015 70.032 6

3.3 Volterra自適應(yīng)濾波預(yù)測

(8)

式中，hp(m1,m2,…,mp)為p階Volterra核;m為濾波器的輸入維數(shù);N1=N2=m≥2D+1。

文獻[20]指出，要實現(xiàn)一個混沌時間序列完全描述出原動力系統(tǒng)的動力特征，至少需要m≥2D+1個變量，所以可取N1=N2=m≥2D+1。

定義U(n)和H(n)分別為濾波器的輸入矢量和系數(shù)矢量:

U(n)=[1,x(n),x(n-1),…,x(n-m+1),

x2(n),x(n)x(n-1),…,x2(n-m+1)]T

(9)

H(n)=[h0,h(0),h(1),…,h(m-1),h2(0,0),

h2(0,1),…,h2(m-1,m-1)]T

(10)

因此，式(7)可表示為

(11)

對于二階Volterra自適應(yīng)濾波器，可采用的自適應(yīng)算法為時間正交自適應(yīng)算法。

本文分別使用A1,B3,A5,B7四組時間序列的前472個數(shù)據(jù)為訓(xùn)練樣本，后20個樣本為測試樣本，Volterra模型沉降預(yù)測結(jié)果如圖11所示。

由圖11可以看出，二階Volterra自適應(yīng)濾波預(yù)測結(jié)果在短期內(nèi)可以很好地反映出真實值的變化規(guī)律，預(yù)測絕對誤差在可接受范圍內(nèi)波動，預(yù)測精度較高，但隨預(yù)測步數(shù)的增加，預(yù)測效果逐漸下降。

為進一步評價沉降預(yù)測結(jié)果的優(yōu)劣，采用下述3個指標(biāo)進行度量:

(1)平均絕對誤差

(12)

(2)平均相對誤差

(13)

(3)均方根誤差

(14)

其中，n為預(yù)測樣本數(shù);xd(i)和x(i)分別為i時刻對應(yīng)的預(yù)測沉降值和真實沉降值。經(jīng)計算，預(yù)測模型各指標(biāo)見表3。

表3 Volterra自適應(yīng)預(yù)測預(yù)測誤差
Table 3 Error of Volterra self-adaptive prediction%

時間序列前10步MAEMAPERMSE后10步MAEMAPERMSEA10.590.060.7616.471.7020.86B34.651.005.7722.123.8529.64A51.590.672.2054.0714.1063.69B73.572.223.9141.5811.7953.14

基于表3的統(tǒng)計結(jié)果，得出二階Volterra自適應(yīng)濾波對4組時間序列的前10步預(yù)測結(jié)果精度都很高，平均絕對誤差(MAE)、平均相對誤差(MAPE)和均方根誤差(RMSE)均在6%以下，實現(xiàn)了沉降時間序列的高精度預(yù)測，分析其原因在于二階Volterra自適應(yīng)濾波避免引入其他干擾因素數(shù)據(jù)，僅從沉降數(shù)據(jù)本身出發(fā)，綜合利用了訓(xùn)練樣本中的線性與非線性因素。而隨著預(yù)測步數(shù)的增加，如后10步的預(yù)測，其MAE,MAPE和RMSE相比于前10步大幅上升，表明預(yù)測結(jié)果精度的大幅下降。這證明二階Volterra自適應(yīng)濾波僅適用于沉降時間序列的短期預(yù)測，將其應(yīng)用于長期沉降預(yù)測結(jié)果不可靠。

4 結(jié) 論

(1)引入小基線集(SBAS)技術(shù)利用sentinel-1 A數(shù)據(jù)監(jiān)測內(nèi)排土場沉降，因內(nèi)排土場植被覆蓋度極低，sentinel-1 A數(shù)據(jù)時間分辨率高，多余觀測多，所得沉降數(shù)據(jù)精度和可信度高，為露天礦內(nèi)排土場沉降監(jiān)測提供了新思路。

(2)內(nèi)排土場沉降剖面呈現(xiàn)明顯的半漏斗狀，總體上累積沉降量與到礦坑的距離成反比，通過累積沉降量將內(nèi)排土場劃分為了8個區(qū)域，通過對8個區(qū)域上隨機兩點的沉降速率分析可得，Ⅰ,Ⅱ區(qū)域處于沉降活躍期，存在滑坡、泥石流風(fēng)險，是后期沉降監(jiān)測的重點區(qū)域，Ⅲ～Ⅶ區(qū)域步入穩(wěn)定過渡期，Ⅷ區(qū)域已基本穩(wěn)定，初步判斷該區(qū)已基本滿足了土地復(fù)墾及建設(shè)簡單構(gòu)筑物的基本要求。

(3)依據(jù)A1～A8的覆土?xí)r間差，對A1累積沉降量進行曲線擬合，得A1最終下沉量約為662 mm，下沉量達到相對穩(wěn)定約需1 320 d，定義內(nèi)排土場下沉系數(shù)為地表最終沉降量與初始覆土高度的比值，求得觀測周期內(nèi)，內(nèi)排土場下沉系數(shù)約為0.639 cm/m。值得注意的是，A1點沉降監(jiān)測時間受限于原始數(shù)據(jù)，略晚于真實覆土?xí)r間，而覆土初期是沉降劇烈期，因此本文求出的下沉系數(shù)小于實際下沉系數(shù)，只能作為后續(xù)研究的參考。

(4)針對SBAS技術(shù)得到的沉降量時間序列，從中任選4點，分別求取其最佳延遲時間τ和嵌入維數(shù)m，經(jīng)檢驗4組序列數(shù)據(jù)均具有混沌特征。

(5)運用混沌理論中的相空間重構(gòu)理論結(jié)合二階Volterra自適應(yīng)濾波對沉降量進行預(yù)測，預(yù)測結(jié)果在短期內(nèi)與真實值吻合較好，但隨預(yù)測步數(shù)增加，預(yù)測誤差逐步增加。表明二階Volterra自適應(yīng)濾波可對SBAS獲取的一維沉降觀測數(shù)據(jù)進行有效地短期預(yù)測，而其對沉降的長期預(yù)測結(jié)果不可靠。