基于SBAS-InSAR與Prophet模型的引黃濟青沿線地表形變監測與預測

田雨情 劉國林 高騰飛 陶秋香

1 山東科技大學測繪與空間信息學院,青島市前灣港路579號,266590 2 山東省地質測繪院,濟南市二環東路11101號,250014

受區域地質構造活動、礦區和城市地下水開采等多種因素的影響,引黃濟青沿線干渠的渠堤及堤岸基礎存在不同程度的地表形變區。傳統監測地表形變的方法如水準測量、GPS等作業效率低、勞動強度大,只能獲得個別點的地表形變信息,難以精確獲得區域整體的形變分布特征[1]。InSAR及其改進技術具有高空間分辨率、范圍廣、全天候、全天時觀測等優勢,已成為地表形變監測領域的一種新的空間對地觀測方法[2-4]。

目前,基于InSAR技術對引黃濟青沿線地表形變進行監測的研究較少,沒有宏觀上獲取的引黃濟青沿線地表形變的位置、時空分布與演變過程的信息。基于此,本文首先利用時序SBAS-InSAR技術處理覆蓋研究區的38景Sentinel-1 SAR影像數據,獲得2019-02～2021-10引黃濟青沿線的年平均形變速率、各成像時刻的累積地表形變量、最終累積地表形變量等信息,提取特征點的地表形變時間序列,并結合實際情況對沿線地表形變的時空演化特征進行分析;然后,在地表形變明顯區域進行緩沖區處理,提取剖面線上高相干像元的地表形變值,并計算其形變梯度,分析引黃濟青干渠的穩定性;最后,利用Prophet模型對特征點的地表形變進行預測分析,并對模型預測的精確度和可靠性進行驗證。

1 研究區及數據來源

1.1 研究區概況

引黃濟青工程位于山東省中部,起點位于博興縣沉沙池出口,途經濱州、東營、濰坊、青島,最終引入青島即墨棘洪灘水庫,全長291.14 km,是一項長距離、跨流域的大型供水工程。本文研究區為引黃濟青路線中段(圖1矩形框)。

圖1 研究區示意圖Fig.1 Sketch map of study area

1.2 實驗數據

選取覆蓋研究區的38景Sentinel-1 SAR影像,時間跨度為2019-02-16～2021-10-03,影像極化方式為VV,波段為C波段。

此外,收集NASA提供的30 m的SRTM DEM高程數據和精密軌道數據,分別用于消除地形相位引起的誤差和軌道誤差對形變監測結果的影響。

2 原理與方法

2.1 SBAS-InSAR原理

假設已獲取時間依次為t0,…,ti,…,tN的N+1幅SAR影像數據,設定空間基線和時間基線閾值進行差分干涉,獲得M幅差分干涉圖,M滿足的條件為:

(1)

假設所有的差分干涉圖都經過正確的相位解纏,且解纏后的相位均被校正到某個形變趨勢穩定或者形變信息已知的高相干點像元上。記時刻ti(i=1,…,N)相對于初始時刻t0的差分干涉相位為φ(ti),φ(ti)為未知參數,通過數據處理得到的差分干涉相位記為δφ(tk)(k=1,…,M),則有2個時間序列:

φ(ti)=[φ(t1),…,φ(tN)]T

(2)

δφ(tk)=[δφ(ti),…,δφ(tM)]T

(3)

對于第k(k=1,…,M)幅差分干涉圖中的任一像元(x,r)有:

Δφk(x,r)=φ(tB,x,r)-φ(tA,x,r)≈

(4)

式中,λ為雷達波長,tA和tB分別為第k幅差分干涉圖對應的主、從影像的獲取時刻,d(tA,x,r)和d(tB,x,r)分別為tA和tB時刻像元(x,r)相對于初始時刻t0在LOS方向上的累積量形變。其中的φ(tA,x,r)可寫為:

(5)

對任一像元,可以列出線性形變模型:

Aφ=δΦ

(6)

式中,A是一個M×N維矩陣,可以直接由已知的差分干涉圖得到。

在實際情況中,N+1景SAR影像通常會被劃分到若干個短基線集中(即L≥2),此時,矩陣A的秩為N-L+1,即ATA是個奇異矩陣,因此式(6)有無窮多個解。SBAS-InSAR采用奇異值分解方法將多個基線集聯合,求其最小范數最小二乘解。

2.2 地表形變梯度

引黃濟青工程作為山東省的重點工程之一,其穩定性至關重要。本文使用形變梯度來描述引黃濟青工程地表形變在不同方向上的形變速率,可以展現出地表形變的不均勻性[5-6],而不均勻地表形變所帶來的危害高于均勻形變帶來的危害[7-8]。

在引黃濟青沿線等間距選取高相干像元,提取高相干像元的形變信息,計算2個高相干像元之間的形變梯度,即

(7)

式中,Kz,z+1是第z個和第z+1個像元之間的地表形變梯度,xz+1、xz分別為第z個和第z+1個像元的地表形變值,s為第z個和第z+1像元之間的距離。

2.3 Prophet模型

Prophet模型是一種結合時間序列分解和機器學習擬合的時序模型,其在建模過程中考慮了趨勢性、節假日效應、周期性以及外部變量等因素的影響,預測效果較好,相對于傳統時間序列模型有明顯優勢[9]。

Prophet模型可以表示為:

y(t)=g(t)+s(t)+h(t)+εt

(8)

式中,g(t)為趨勢項,表示時間序列非周期的變化趨勢,s(t)為周期項,表示時間序列中呈現周期性變化的部分,h(t)為節假日項,表示時間序列中那些潛在的具有非固定周期的節假日對預測值造成的影響,εt為誤差項,服從高斯分布。

基于邏輯回歸函數來模擬趨勢項:

g(t)=

(9)

式中,a(t)=[a1(t),a2(t),…,as(t)],δ=[δ1,δ2,…,δs]T,γ=[γ1,γ2,…,γs]T,C(t)為承載量,k為增長率,m為偏移量。

使用傅里葉級數來模擬時間序列的周期性,假設P表示時間序列的周期,則傅里葉級數的形式為:

(10)

式中,β=[a1,b1,a2,b2,…,aN,bN],為Prophet模型的初始化參數,其初始化參數β～N(0,σ2),σ越大,表示季節效應越明顯,σ越小,表示季節效應越不明顯。

假設有L個節假日,那么節假日效應模型為:

(11)

式中,Z(t)=[1{t∈Di},…,1{t∈DL}],k=[k1,…,kL]T。

2.4 數據處理流程

使用GAMMA軟件[10]處理SAR影像。根據綜合相干系數最小的原則[11],選取一景影像作為主影像,其余N景SAR影像作為從影像與主影像進行配準。對配準好的影像進行去斜、裁剪等處理;設置時空基線閾值,生成100個干涉像對,處理得到時空基線分布圖(圖2)。設置距離向和方位向為4∶1的多視處理來抑制噪聲,與已有DEM數據進行差分干涉處理;使用Goldstein濾波方法對各差分干涉圖進行濾波處理,得到M幅濾波增強后的時序差分干涉圖。使用最小費用流法(minimum cost flow, MCF)對濾波增強后的各差分干涉圖進行相位解纏,對SBAS-InSAR的結果進行后續處理,得到研究區形變信息。基于形變梯度對引黃濟青沿線地表的穩定性進行分析,然后利用Prophet模型對沿線若干特征點的地表形變進行預測。圖3給出了主要的數據處理流程與方法。

165 二仙湯抗骨質疏松有效組分對維甲酸致骨丟失大鼠的影響 張建花，沈 燚，何玉瓊，韓 婷，秦路平，張巧艷

圖2 時空基線分布Fig.2 Spatiotemporal baseline distribution

圖3 SBAS-InSAR數據處理流程Fig.3 SBAS-InSAR data processing flow chart

3 結果與分析

3.1 引黃濟青沿線地表形變分析

計算獲取2019-02-16～2021-10-03的地表形變信息,圖4為年平均形變速率圖,形變速率為負值表示地面沉降,正值表示地面抬升。圖5為引黃濟青沿線地表形變折線圖,橫坐標為從起始點到形變點的距離(以研究區左側為起始點)。

圖4 引黃濟青區域年平均形變速率Fig.4 Annual average deformation rate of the Yellow river to Qingdao diversion region

圖5 引黃濟青沿線地表形變速率折線圖Fig.5 Line chart of surface deformation rate along the Yellow river to Qingdao diversion route

從圖4可以看出,大部分區域年平均形變速率為正值,表現為地表抬升;稻莊鎮、大碼頭鎮、大王鎮、營里鎮、羊口鎮、稻田鎮等區域形變速率為負值,表現為地表沉降,沉降速率多在-60～0 mm/a之間,沉降最嚴重的廣饒縣最大形變速率達到-167 mm/a。從圖5看出,引黃濟青沿線的形變速率主要在-60～40 mm/a之間,形變速率波動與年平均形變速率圖形變區域較吻合。

造成引黃濟青沿線地表形變的因素可能包括城市地面沉降、開采沉陷、人類工程活動和數據處理誤差等。從圖4可以看出,地表形變集中在廣饒縣境內。查閱資料可知,造成廣饒縣地表形變的主要因素為地下水超采、城市建設等。中深層地下水長期超采形成地下水降落漏斗,加上城市建設,特別是高層建筑的興建,對其下部地層持續加載,造成上部松散地層固結壓縮,引起大范圍的地面沉降,進而對引黃濟青沿線造成影響。

在沉降中心選取4個特征點P1、P2、P3、P4(圖4),提取其地表形變時間序列,結果見圖6,其中,P1、P2、P3位于廣饒縣境內,P4位于壽光市營里鎮。

圖6 特征點累積形變時序折線圖Fig.6 Line chart of cumulative deformation time series of feature points

從圖6可以看出,各個特征點累積形變量基本處于持續沉降趨勢,沉降曲線波動較小,趨勢相似,表現出較強的線性特征。至2021-10-03四個特征點累積形變值分別達到-320 mm、-239 mm、-238 mm和-212 mm,其中,P1點形變量明顯大于其他3點,其地理位置在廣饒縣城區。廣饒縣地下水長期處于超采狀態,形成地下水降落漏斗,造成該區域出現地表沉降,代表區域包括廣饒縣城區、稻莊鎮和大王鎮等。

3.2 引黃濟青沿線地表形變穩定性分析

根據圖4,在地表形變明顯區域內對引黃濟青沿線2 km范圍進行緩沖區處理,獲得如圖7所示的緩沖區累積形變圖,其中,a、c為縱向剖面線,b、d為橫向剖面線。分別獲取橫縱剖面線累積形變值,以此為依據計算剖面線的形變梯度,結果見8,其中,圖8(a)為剖面線累積形變量,8(b)為剖面線形變梯度。

圖7 2 km緩沖區累積沉降圖Fig.7 Cumulative settlement map of 2 km buffer zone

圖8 累積形變值及形變梯度折線圖Fig.8 Line chart of cumulative deformation value and deformation gradient

對比分析圖8(a)、8(b)可以看出,剖面線a出現一次形變梯度值的驟升,與形變嚴重區域相吻合,剖面線c形變梯度值比較穩定;橫向剖面線b、d各像元點之間存在較大的沉降波動,相應的形變梯度值變化也較大,整體上地表形變嚴重區域和形變梯度較大區域有較大的重疊。由此可知,縱向剖面線的地表形變速率較穩定,不均勻形變較少;橫向剖面線存在較多的不均勻形變,可能是由引黃濟青沿線周圍的地質活動引起的,應多加關注。

引黃濟青沿線橫縱剖面線形變梯度穩定性分析表明,引黃濟青沿線整體處于一個比較穩定的狀態。

3.3 基于Prophet模型的引黃濟青地表形變預測

圖9 SBAS-InSAR觀測值與后8期Prophet模型預測值對比Fig.9 Comparison between SBAS-InSAR observation values and predicted values of Prophet model in the last 8 periods

從圖9中可以看出,4個特征點后8期數據的預測值與SBAS觀測值均具有較好的一致性,P1點預測結果與觀測值最吻合;P2、P3點預測值與觀測值之間差異相對大一些,但誤差值均小于15 mm;P4點預測結果與觀測值存在一些差異,但整體沉降趨勢一致。

從圖10看出,P1點沉降值逐步增加,沉降趨勢接近指數函數形式,到2022-04沉降值達到-375 mm;P2～P4點訓練樣本具有一定的周期波動性,預測數值同樣存在波動,沉降值不是單純的增加,其中,P2點出現2次小幅抬升,沉降值較2021-10增加30 mm,P3、P4點形變值上下波動,總體沉降值變化不大,分別增加10 mm和4 mm。

圖10 SBAS-InSAR觀測值與未來8期Prophet模型預測值對比Fig.10 Comparison between SBAS-InSAR observation values and predicted values of Prophet model in the future 8 periods

綜上,Prophet模型能夠在數據量和地表形變量較小時進行較好的模擬和預測。

4 結 語

本文利用SBAS-InSAR技術對2019-02-16～2021-10-03覆蓋引黃濟青沿線的38景Sentinel-1數據進行處理,獲得研究區大范圍、長時間的地表形變信息;基于地表形變結果,針對沉降嚴重區域進行穩定性分析,并利用Prophet模型預測特征點未來的形變趨勢,得到結論如下:

1)研究區域內大部分地區地表形變速率和累積形變量都較小,處于穩定狀態;大部分區域形變速率位于-60～40 mm/a之間,主要表現為地面抬升;沉降區主要分布在廣饒縣等地方,最大沉降速率達到-167 mm/a,最大累積沉降值為-366 mm。

2)地表形變嚴重區域與形變梯度較大區域高度重合,并且橫向剖面線形變梯度變化大于縱向剖面線,應給予更多關注。

3)Prophet模型在數據量和形變量較小的情況下預測效果較好。