基于SBAS-InSAR與GM(1,1)模型的張博線采空區地表形變監測

劉瑞祥, 陶秋香*, 劉曉朋, 韓宇, 李雪鵬, 肖怡欣

(1.山東科技大學測繪與空間信息學院, 青島 266590; 2.中鐵工程設計咨詢集團有限公司濟南設計院, 濟南 250022; 3.山東省青島市勘查測繪研究院, 青島 266033)

當今經濟社會的快速發展需要大量的煤炭資源作為支撐,但在開發地下資源的過程中,形成了大量地下采空區,使原有的土質結構遭到破壞,很容易造成地面塌陷等人為災害[1]。而擬建的張博線鐵路沿線(即淄博市張店區至博山區路段)恰好經過了大量采空區,若采空區上覆地表尚未穩定,存在較嚴重的地表形變,就會對鐵路建設以及運營造成安全隱患,給人們的出行帶來危險。

盡管水準測量、全球導航衛星系統(global navigation satellite system, GNSS)測量等技術手段都能夠進行精確的地表形變監測,但卻具有一定的局限性,例如,人力、物力投入較多,難以實現大規模的覆蓋和普查[2]。近年來,隨著空間對地觀測技術的發展,合成孔徑雷達差分干涉測量(differential interferometric synthetic aperture radar, D-InSAR)技術[3-4]和小基線集合成孔徑雷達差分干涉(small baseline subset interferometric synthetic aperture radar, SBAS-InSAR)技術[5]的發展,以其高效率、低成本、大范圍的優點,極大地彌補了傳統測量的不足,為地表形變監測提供了新的測量方法,之后的研究員利用這些技術取得了一系列的研究成果[6-11]。GM(1,1)灰度預測模型被提出后,為研究人員對地表形變的預測提供了一種極為簡便可行的方法,石曉宇等[12]、劉哲[13]、郭在潔等[14]分別對GM(1,1)在不同情形下的應用進行了研究。

現利用SBAS-InSAR技術對43景Sentinel-1A SAR影像數據進行處理,獲取擬建的張博線鐵路沿線的地表形變信息;以張博線鐵路沿線為中心線,對其左右兩側1 km緩沖區內總體地表形變情況、重點形變區域及其里程點的地表形變情況等進行分析;利用GM(1,1)灰度預測模型對地表特征點監測結果進行模擬和預測,以研究GM(1,1)在采空區地表形變監測中預測的可行性。

1 研究區與數據源

1.1 研究區概況

張博線位于山東省淄博市境內,北起膠濟線淄博站,沿張南路南行至南定鎮,轉向西南,跨越省道S703及國道G205后沿G205至淄川區,向西跨越孝婦河,繼續沿G205直至設計終點博山站。線路經過區域有淄博市張店區、淄川區及博山區,線路全長39.233 km。沿線村鎮密集。圖1為張博線地理位置及其經過的主要采空區。

圖1 研究區地理位置Fig.1 Geographical location of study area

1.2 數據源

Sentinel-1SAR衛星為載有C波段的合成孔徑雷達,具備干涉寬帶(interferometric wide swath,IW)等多種工作模式,是歐洲航天局哥白尼計劃中的極其重要的地球觀測衛星。該衛星具有水平發射水平接收(VV)、水平發射垂直接收(VH)等不同的極化方式,能提供連續的、全天候的SAR影像,有著高重訪頻率、高覆蓋能力等無可比擬的的優點,以及極好的時效性和可靠性[15]。現選取覆蓋張博線,時間跨度為2019年1月5日—2022年1月1日的43景C波段、VV極化的哨兵1A升軌數據。表1為43景SAR影像的獲取日期。

表1 Sentinel-1A SAR影像數據情況Table 1 The data information of Sentinel-1A SAR images

2 原理與方法

2.1 SBAS-InSAR技術原理

Berardino等[16]在2002年提出的SBAS-InSAR技術克服了大量傳統InSAR的諸多弊端,有效地提高了地表形變監測的精度。

設研究區域有N+1幅雷達影像,它們獲取的時間分別為t0、t1、…、tn,現在選擇其中一幅影像作為公共主影像,任何影像都能與其他至少一幅影像進行干涉,可以得到M幅干涉圖,M滿足如下關系式。

(1)

利用在tA和tB時刻獲得的兩景已經去除地形相位的SAR影像(tA

δφi(x,r)=φ(tB,r,x)-φ(tA,r,x)

(2)

式(2)中:φ(tA,r,x)和φ(tB,r,x)分別為像元(r,x)在tA與tB時刻相對于參考時刻t0的雷達視線向的累計形變相位;λ為雷達的波長;d(tA,r,x)和d(tB,r,x)分別為tA與tB時刻相對于t0的雷達視線向的累計形變量,d(t0,r,x)≡0。有

(3)

假設IE=[ IE1IE2… IEM]、IS=[ IS1IS2… ISM]分別為干涉數據處理時按照時間順序排列的主、輔影像序列,并滿足條件為

IEj>ISj,j=1,2,…,M

(4)

那么所有差分干涉圖相位就可以組成的觀測方程為

δφi=φ(tIEi)-φ(tISi)

(5)

將式(5)轉換為線性方程可得

Aφ=δφ

(6)

式(6)中:A為M×N維矩陣,當所有的SAR影像都被分在一組,且M≥N,矩陣A的秩是N,通過對式(6)采用最小二乘法求解,可得

(7)

2.2 GM(1,1)模型

GM(1,1)模型是一種用于研究數據量與信息較少的不確定性問題的理論方法,是鄧聚龍教授于1982年提出的一種灰色系統理論,利用該理論可以利用較少的數據來獲得較好的預測效果[17]。原理如下。

設有原始觀測數列為

X(0)=[X(0)(1),X(0)(2),…,X(0)(n)]

(8)

將原始觀測數列進行一次累加,得到新的數列為

X(1)=[X(1)(1),X(1)(2),…,X(1)(n)]

(9)

對數列X(1)建立一階微分方程可得

(10)

式(10)中:a為發展灰數,反映X(1)及X(0)的發展趨勢;u為內生控制灰數,反映了數據間變化的關系。

k=2,3,…,n

(11)

令

Z(1)(k)=μX(1)(k)+(1-μ)X(1)(k-1),

k=2,3,…,n

(12)

式(12)中:Z(1)(k)被稱為式(10)的背景值;μ∈[0,1],μ為權重系數。

假定μ取值為0.5,則有

(13)

將式(10)離散化,則有

X(1)(k)-X(1)(k-1)+aZ(1)(k)=u,

k=2,3,…,n

(14)

利用最小二乘法求解式(14)可得

(15)

求得a和u后,繼續求解式(10),得

(16)

將式(16)離散化得

(17)

(18)

代入式(17)得

k=0,1,…,n-1

(19)

對式(19)做累減還原便可得到X(0)的灰色預測模型為

k=2,3,…,n-1

(20)

2.3 數據處理流程

數據處理采用的平臺是世界著名的瑞士 GAMMA 遙感公司開發的專門用于處理SAR數據的商業軟件:GAMMA,采用的方法是SBAS-InSAR技術,充分利用現有43景sentinel-1A數據,得到最優結果。首先,利用連續的常規差分方案,通過設定的時空基線指標使43景數據進行自由配對組合的方式,并選取其中2020年1月24日獲取的SAR影像作為公共主影像,然后從剩下的42個單主影像組成的干涉對中選取時空基線指標較好的像對,聯合外部參照 DEM 處理為雷達坐標系下的模擬解纏地形相位,之后再進行常規差分處理, 獲取區域內形變概略位置,以及形變區的演變過程,據此提取強度閾值、相位相干性閾值等經驗參數;然后進行SBAS-InSAR 處理,通過時空基線指標進行自由配對組合,以時間基線小于96 d,空間基線小于 194.9 m 為時空臨界基線,選取146個組合對,進行差分干涉處理,得到差分干涉圖;對差分干涉圖進行濾波處理,去除大氣延遲誤差、軌道誤差、殘余DEM誤差等;對濾波后的差分干涉圖進行相位解纏;利用解纏后的干涉圖進行精密基線估計;然后利用估計的精密基線再次進行“差分干涉處理→濾波處理→差分干涉圖解纏”處理;對第二次獲取的解纏后的差分干涉圖進行瀏覽,剔除依然存在軌道誤差、大氣延遲誤差的差分干涉圖,形成最優組合的解纏后的差分干涉圖序列;最后進行地表形變反演和歸一化處理,得到研究區的地表形變時間序列;對求取的形變結果進行克呂格插值方法分析,基于參考點進行整體修正,得到最終形變結果。具體的數據處理流程如圖2所示。

圖2 數據處理流程Fig.2 Data processing flow

3 實驗結果分析

3.1 鐵路沿線沉降分析

利用圖2所示的數據處理流程和方法對表1的43景Sentine-1A SAR數據進行處理,獲取2019年1月—2022年1月張博線鐵路沿線的地表形變信息,對張博線鐵路沿線左右500 m緩沖區內的地表形變進行分析,得到如圖3所示的鐵路沿線形變速率圖和如圖4所示的1 km緩沖區內最終累計形變圖。

圖3 鐵路沿線1 km緩沖區內形變速率Fig.3 The deformation rate in the 1 km buffer zone along the railway

分析圖3和圖4可以看出,整個鐵路沿線在監測周期中的地表形變以沉降為主,最大沉降速率與累計沉降量分別為-7.9 mm/a與-23.8 mm;小部分區域存在微小抬升,最大抬升速率與累計抬升量為2.7 mm/a與7.7 mm。其大部分區域年平均形變速率在-2.0～2.0 mm/a,累積形變量集中在 -5.0～5.0 mm,整個區域相對穩定;結合礦區范圍來看,各礦區域范圍內都存在沉降區,且沉降范圍與沉降量各不相同。除此之外,在不存在采空區的十里鋪村西南側存在范圍較大的持續沉降區域,應進行實地調查以查明原因。

為了更好地分析鐵路沿線的形變區域進行分析,將形變范圍分成了3個等級,大于等于0 mm、小于等于0 mm、小于等于-5 mm共3段,如圖5所示。

圖5 監測周期內3個區間內累計地表形變區域分布Fig.5 The regional distribution of cumulative surface deformation in 3 intervals during the monitoring period

分析圖5,存在沉降的區域主要位于GDK1～GDK11、GDK16～GDK21、GDK28～GDK30、GDK31～GDK34、GDK37～GDK39,沉降區域的面積為 18.95 km2;沉降較嚴重(沉降量多于-5.0 mm)的區域主要位于GDK16～GDK19區域,在 GDK23～GDK24、GDK8～GDK9附近也有少許,通過進行面積統計發現,沉降小于等于-5.0 mm的區域的面積為3.88 km2,小于整個監測區域的1/10,整體相對穩定;GDK3～GDK5、GDK11～GDK16、GDK21～GDK24、GDK25～GDK28、GDK34～GDK37等附近區域未沉降或略有抬升的區域面積為21.22 km2。

分別提取鐵路中線以及左右兩側各500 m在監測周期內的累計形變量,制作剖面圖如圖6所示。

圖6 鐵路沿線中心及其兩側各500 m沿線形變速率剖面圖Fig.6 The deformation rate profile along the center of the railway and 500 m on both sides

分析圖6,在2019年1月—2022年1月期間,鐵路沿線的形變速率整體都較小,最大不超過 -4.0 mm/a,從中線可以看出,形變速率主要集中在-1.0～1.0 mm/a,其中在16～18 km,沉降速率超過了-1.0 mm/a,最大達到-1.7 mm/a;從形變穩定性上來看,中線整體波動較小,但是在11～23 km 處,形變出現了較大波動,形變速率由 -0.2 mm/a 變化到-1.7 mm/a,再由-1.7 mm/a變化到-0.2 mm/a,變化幅度相對于其他路段較大;對于左側500 m,相對于中線來說,鐵路沿線的形變速率和形變波動幅度,都有了一定程度的增大,在6～9 km處出現了最大形變速率,為-2.6 mm/a,在16～18 km,沉降速率均超過了-2.0 mm/a;在形變穩定性上,在6～10 km和23～26 km處,形變呈現出W形起伏,在11～23 km處,形變呈現出V形起伏,且起伏程度相對較大,但考慮到該區域存在河道,應進行實地勘察,排除不穩定沉降對鐵路線的威脅;對于右側500 m,從整體來看,大部分地區形變速率相對較為穩定,形變速率在-2～1.3 mm/a, 在16.4 km處出現了最大沉降速率,為-4.0 mm/a;在15～21 km處(光正實業西北部河道內),形變呈現出W形起伏,且起伏狀態較大,穩定較差,須進行實地排查。

3.2 GM(1,1)模型在地表形變監測中的應用

選擇沉降較為嚴重的區域選擇P1、P2、P3三點,以及出現抬升的P4點,如圖7所示。提取此4個點在43個時期內的累計沉降值,以前40期數據為初始數據建立GM(1,1)模型,并用該模型對后三期進行預測,與實際數據相比較。

圖7 選取測試點的位置Fig.7 The location of the selected test point

預測結果如圖8～圖11所示。其中,P1點的發展灰數a=-0.037 8;P2點的發展灰數a=-0.046 9;P3點的發展灰數a=-0.050 1;P4點的發展灰數a=-0.036 2。

圖8 P1點GM(1,1)預測模型與監測值對比Fig.8 The GM (1,1) prediction model of P1 point is compared with the monitoring value

圖9 P2點GM(1,1)預測模型與監測值對比Fig.9 The GM (1,1) prediction model of P2 point is compared with the monitoring value

圖10 P3點GM(1,1)預測模型與監測值對比Fig.10 The GM (1,1) prediction model of P3 point is compared with the monitoring value

由文獻[18]可知:

(1)當-a<0.3時,GM(1,1)可用于中長期預測。

(2)當0.3<-a<0.5時,GM(1,1)可用于短期預測,中長期預測慎用。

(3)當0.5<-a<0.8時,用GM(1,1)作短期預測應十分謹慎。

(4)當0.8 <-a<1時,應采用參差修正GM(1,1)模型。

(5)當-a>1時,不宜采用GM(1,1)模型。

即GM(1,1)模型在本次地表形變監測中可以用于中長期預測。

4 結論

利用43景Sentinel-1A數據,基于SBAS-InSAR技術獲取了2019年1月5日—2022年1月1日淄博市擬建張博線鐵路沿線地表形變信息,對其分析可得出以下主要結論。

(1)2019年1月—2022年1月期間,沉降區域的面積在逐年增加,統計發現三年來最終呈現沉降的區域面積為22.83 km2,占總面積的51.8%,但沉降量級不大,大部分區域年平均形變速率在-2.0～2.0 mm/a,整個區域比較穩定;沉降最嚴重的區域位于GDK23+550左側232.2 m,出現了一個明顯的沉降中心,最大沉降速率達到-7.9 mm/a, 最大沉降量達到-23.8 mm,且呈現出線性沉降趨勢,對其沉降情況需要進一步關注。

(2)大部分里程點的形變較為穩定,形變速率在-1.0 ～ 1.0 mm/a,累計形變量在-3.0～3.0 mm。形變嚴重的點為GDK16、GDK17、GDK18,沉降值分別達到-9.0、-11.0、-11.0 mm,需要進一步關注。

(3)鐵路中線大部分地區處于穩定狀態,鐵路左側500 m和右側500 m沿線形變速率和波動幅度要略大于鐵路中線,需要進一步關注,為鐵路的安全建設與運行提供保障。

(4)GM(1,1)模型在地表形變監測中擬合度較高,表現較好,可用于中長期的預測。