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

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

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

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

1 研究區(qū)及數(shù)據(jù)來源

1.1 研究區(qū)概況

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

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

1.2 實驗數(shù)據(jù)

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

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

2 原理與方法

2.1 SBAS-InSAR原理

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

(1)

假設(shè)所有的差分干涉圖都經(jīng)過正確的相位解纏,且解纏后的相位均被校正到某個形變趨勢穩(wěn)定或者形變信息已知的高相干點像元上。記時刻ti(i=1,…,N)相對于初始時刻t0的差分干涉相位為φ(ti),φ(ti)為未知參數(shù),通過數(shù)據(jù)處理得到的差分干涉相位記為δφ(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幅差分干涉圖對應(yīng)的主、從影像的獲取時刻,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采用奇異值分解方法將多個基線集聯(lián)合,求其最小范數(shù)最小二乘解。

2.2 地表形變梯度

引黃濟青工程作為山東省的重點工程之一,其穩(wěn)定性至關(guān)重要。本文使用形變梯度來描述引黃濟青工程地表形變在不同方向上的形變速率,可以展現(xiàn)出地表形變的不均勻性[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模型是一種結(jié)合時間序列分解和機器學習擬合的時序模型,其在建模過程中考慮了趨勢性、節(jié)假日效應(yīng)、周期性以及外部變量等因素的影響,預測效果較好,相對于傳統(tǒng)時間序列模型有明顯優(yōu)勢[9]。

Prophet模型可以表示為:

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

(8)

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

基于邏輯回歸函數(shù)來模擬趨勢項:

g(t)=

(9)

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

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

(10)

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

假設(shè)有L個節(jié)假日,那么節(jié)假日效應(yīng)模型為:

(11)

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

2.4 數(shù)據(jù)處理流程

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

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

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

圖3 SBAS-InSAR數(shù)據(jù)處理流程Fig.3 SBAS-InSAR data processing flow chart

3 結(jié)果與分析

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

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

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

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

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

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

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

3.2 引黃濟青沿線地表形變穩(wěn)定性分析

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

圖7 2 km緩沖區(qū)累積沉降圖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出現(xiàn)一次形變梯度值的驟升,與形變嚴重區(qū)域相吻合,剖面線c形變梯度值比較穩(wěn)定;橫向剖面線b、d各像元點之間存在較大的沉降波動,相應(yīng)的形變梯度值變化也較大,整體上地表形變嚴重區(qū)域和形變梯度較大區(qū)域有較大的重疊。由此可知,縱向剖面線的地表形變速率較穩(wěn)定,不均勻形變較少;橫向剖面線存在較多的不均勻形變,可能是由引黃濟青沿線周圍的地質(zhì)活動引起的,應(yīng)多加關(guān)注。

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

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期數(shù)據(jù)的預測值與SBAS觀測值均具有較好的一致性,P1點預測結(jié)果與觀測值最吻合;P2、P3點預測值與觀測值之間差異相對大一些,但誤差值均小于15 mm;P4點預測結(jié)果與觀測值存在一些差異,但整體沉降趨勢一致。

從圖10看出,P1點沉降值逐步增加,沉降趨勢接近指數(shù)函數(shù)形式,到2022-04沉降值達到-375 mm;P2～P4點訓練樣本具有一定的周期波動性,預測數(shù)值同樣存在波動,沉降值不是單純的增加,其中,P2點出現(xiàn)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模型能夠在數(shù)據(jù)量和地表形變量較小時進行較好的模擬和預測。

4 結(jié) 語

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

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

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

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