田 雨，雷少剛，卞正富

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

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

1 研究區(qū)概況

勝利露天煤礦隸屬于神華北電勝利能源有限公司，其中勝利一號露天礦位于勝利煤田的中西部，錫林浩特市北郊5 km，本煤田位于內蒙古高原中部，氣候干燥，丘陵地形，地形標高970～1 326 m，露天礦共有4個排土場，其中3個外排土場:南排土場、北排土場、沿幫排土場。礦區(qū)至2013年實現(xiàn)全部內排，內排土場覆蓋范圍約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)技術由BERARDINO等于2002年提出[7]，其基本原理是將獲得的同一個區(qū)域的SAR數(shù)據(jù)按照空間基線組成多個集合，利用最小二乘方法分別得到每個集合的沉降時間序列，然后利用奇異值分解法(SVD)將多個集合聯(lián)合求解從而得到整個觀測時間周期的沉降時間序列[8]，經(jīng)眾多國內外學者驗證，它可以有效抑制由于時間、空間造成的失相干和大氣效應，尤其在低植被覆蓋區(qū)域可得到高精度的沉降時序數(shù)據(jù)。

2.2 數(shù)據(jù)準備

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模式下的單視復數(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 沉降結果及分析

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

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

圖3 沉降量Fig.3 Maps of settlement

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

2.4 下沉系數(shù)估計

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

(1)

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

依據(jù)內排土場排土時序圖，得A1～A8，B1～B8覆土時間見表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

因內排土場覆土層力學性質相近，覆土高度幾乎相同，所以可近似認為A1～A8，B1～B8為內排土場覆土工后不同時期的沉降狀態(tài)，其中A1接近為監(jiān)測周期開始時最新的覆土點，選擇A1～A8累積沉降量依據(jù)時間關系進行曲線擬合，估計A1點累積下沉量，結果如圖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，因此可求得觀測周期內，內排土場下沉系數(shù)約為0.639 cm/m。值得注意的是，A1點沉降監(jiān)測時間受限于原始數(shù)據(jù)，略晚于真實覆土時間，而覆土初期是沉降劇烈期，因此本文求出的下沉系數(shù)小于實際下沉系數(shù)，只能作為后續(xù)研究的參考。

3 沉降預測

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

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

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

3.2 相空間重構

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

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

對于時間序列{x(t1)，x(t2)，…，x(tn)}，當延遲時間為τ時，序列變?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各點的相鄰兩期累積沉降量差值構成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對應a2(i，m)的均值，即

(4)

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

(5)

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

(6)

(7)

若E2(m)的值恒為1，則表明所選時間序列為隨機時間序列，數(shù)據(jù)間無相關性;若其不恒為1，則表明時間序列數(shù)據(jù)間的相關性受嵌入維數(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ù)，結果見表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自適應濾波預測

(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自適應濾波器，可采用的自適應算法為時間正交自適應算法。

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

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

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

(1)平均絕對誤差

(12)

(2)平均相對誤差

(13)

(3)均方根誤差

(14)

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

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

4 結 論

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

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

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

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

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