洞庭湖軟土區域時序InSAR形變與環境物理參數聯合估計方法

朱 珺,朱凌杰,2,邢學敏,3,張 銳,鮑 亮,張騰飛,鮑皓丹,3

1. 長沙理工大學交通運輸工程學院,湖南 長沙 410114; 2. 中南大學地球科學與信息物理學院,湖南 長沙 410083; 3. 洞庭湖生態環境遙感監測湖南省重點實驗室,湖南 長沙 410114

洞庭湖位于長江中游荊江南岸,跨湘、鄂兩省,是我國第二大淡水湖,湖區所處的洞庭湖生態經濟區作為國家級發展片區,是中部地區崛起的重要戰略支撐點。洞庭湖周邊水系發達,且有湘江貫穿湖泊南北,區域內土質以軟土為主,主要由淤泥、淤泥質黏性土、軟塑狀亞黏土、亞砂土及松散狀粉細砂組成[1]。由于軟土含水量大、可壓縮性高、強度較低、結構松散等一系列不良土質特點,在軟土區域修建的基礎設施更易發生危險性沉降[2-3],累計形變甚至會導致各種構筑物及基礎設施嚴重變形和損毀,是社會生產生活的嚴重安全隱患。因此對該區分布密集的基礎設施進行長期持續穩定性監測,并剖析環境物理因素對時序形變特征的影響具有十分重要的意義。

針對基礎設施的形變監測,傳統的水準測量方法費時費力,成本高昂且難以實現大尺度監測[4-5]。永久散射體差分干涉測量技術(permanent scatterer interferometric synthetic aperture radar,PSInSAR)通過提取測區內長時間保持穩定散射特性的點(permanent scatterers,PS)進行建模分析,解算形變參數,能夠高時效大尺度地獲取測區時序形變結果,同時其在理論上具備亞毫米級的形變反演精度,在城市大型基礎設施形變監測中極具潛力[6-9]。

利用PSInSAR技術進行時序形變反演的關鍵是對高相干點形變分量與形變參數之間的時序函數關系進行建模。準確且符合物理規律的形變模型不僅可以合理地描述形變隨時間演化的規律,提高形變估計的精度,更可反演相應的環境物理參數,為后續形變結果的解譯及災害防治提供參考[10]。目前在PSInSAR技術中,針對城市基礎設施變形監測的形變建模普遍采用簡單的線性模型,即假定干涉組合內形變速率隨時間呈線性變化趨勢。但是,軟土區域內基礎設施形變在時序上具有顯著非線性特征。已有研究表明,溫度變化對基礎設施變形的影響比運行荷載或結構損傷的影響更顯著[11-12]。X波段的TerraSAR數據對熱膨脹有高靈敏度,特別是當建筑設施以金屬成分為主時[11,13],由于溫度變化導致的材料膨脹或收縮會導致InSAR時間序列位移的強烈季節性變化[14]。因此,簡單的線性假設無法揭示這一實際變形規律,不但會影響變形序列的反演精度,還難以獲取可供災害防治分析的環境物理參數。如果將民用基礎設施的熱膨脹特性引入InSAR形變建模和解譯中,將會獲取到更可靠的城市基礎設施變形趨勢信息,為本文研究提供了啟示。

針對上述問題,本文提出一種基于熱膨脹雙曲線模型的時序InSAR形變與環境物理參數聯合估計方法。該方法在PS基線網絡的相位建模環節以軟土區域形變的先驗模型——雙曲線模型為基礎模型,將熱膨脹系數和降水參數融入其中,即分別建立受熱膨脹和外界環境降水影響的相位分量模型。用上述模型取代原線性速率模型,解算時不但可以獲取每個PS點的變形值,還可同時獲取模型中相應的環境物理參數,實現洞庭湖區域基礎設施總形變序列與環境物理參數的聯合反演。獲得的形變序列與環境物理參數可為基礎設施的安全穩定分析和災害防治提供可靠的數據支撐。

1 時序InSAR模型構建

1.1 PSInSAR傳統線性速率模型

PS點提取時,主要是通過對覆蓋同一地區的多幅SAR影像的相位和幅度信息進行統計分析,篩選出具有較高后向散射特性的目標點。為有效減弱大氣延遲相位、噪聲相位的影響,處理中對篩選出來的PS點構建基線網絡,即將任意相鄰兩PS點連接為一條PS基線,然后對每一條PS基線進行時序建模。

對于PS網絡中任意一條基線(分別連接PS點i和j),對這兩點的相位值作差,可建立如下模型[15]

(1)

(2)

1.2 熱膨脹雙曲線模型

(3)

雙曲線模型是一種反映時序非線性關系的函數模型,可以表征沉降速率隨時間表現為先快后慢的變化特征。符合軟土形變隨時間演化的非線性規律,因此被廣泛用于路基路面工程中的軟土次固結沉降預測[17]。雙曲線模型的基本形式為

(4)

(5)

式中,Tm為干涉圖m的從影像獲取時刻。對等號右邊進行泰勒級數展開,取前兩項后,式(5)可表示為[18]

(6)

(7)

(8)

(9)

式中,Tempm為研究區域在SAR影像過境時刻的環境溫度;Precm為研究區域當月降水量;Th和Pr分別為熱膨脹系數和降水參數,為待求模型未知參數。

(10)

式中,a、b為變形參數;Th和Pr為環境物理參數?；谑?10)可對變形和環境物理參數聯合求解。

(11)

式(11)顧及了湖區軟土形變隨時間呈非線性變化的特征,且考慮了基礎設施受熱膨脹和外界環境降水的影響,具有更為明確的物理意義。

2 模型的未知參數估計

2.1 基于LAMDBA-SVD的時間維參數估計

將式(11)構建的熱膨脹雙曲線模型代入式(1)中,則任一基線邊上兩PS點相位差可寫為

(12)

完成對式(12)的時間維參數估計后,為衡量基線參數估值質量,采用PS基線邊的時序相關系數作為質量指標[4]?；€邊的時序相關系數可表示為

(13)

2.2 基于雅可比迭代的空間維參數估計

在獲取了形變參數增量的估值之后,需要以此為觀測值,恢復出每個PS點的絕對參數值,這一過程即為空間維參數估計。在此,利用間接平差的方法對整個PS基線網絡來實現空間維解纏[22]。在求解過程中,需要在研究區域中選取一穩定的參考點以解決觀測方程的法方程系數陣秩虧問題。由于對PS基線網絡構建的觀測方程,其法方程的系數陣B是一個大型稀疏矩陣,如果直接利用間接平差法進行求逆操作需要耗費大量計算時間,且求解結果不可靠。因此,本文采用雅可比迭代法實現法方程系數陣的迭代求逆。雅可比迭代是求解線性方程組中一種常用的迭代方法[23],主要用于求解帶有大型稀疏矩陣的線性方程組,主要包括以下幾個步驟。

(1) 對觀測方程的系數矩陣進行LDU分解:B=D+L+U,其中D為對角陣,L為下三角陣,而U為上三角陣。

(2) 分解后原方程可轉化為

X=B0X+f

(14)

式中,B0=-D-1(L+U);f=D-1b。由于D為對角陣,對其對角線元素求倒數即可得出D-1,因此極大地降低了運算的壓力。

(3) 給X初值后可采用式(15)進行迭代

(15)

在時間維和空間維參數估計完成之后,即獲取了研究區所有PS點的形變參數值,對這些形變參數在時間域上積分,即得到PS點上的低通形變分量。在PS基線網絡布設時,設置基線邊距離均不超過800 m,以盡最大可能減弱大氣延遲相位的影響。通過對殘余相位進行時空濾波以分離出殘余相位中殘留的高通形變分量部分,再與低通形變疊加,得到整個研究區域的時序形變場。算法流程如圖1所示。

圖1 基于熱膨脹雙曲線模型的洞庭湖軟土區域時序InSAR形變估計處理流程Fig.1 Processing flow of InSAR soft soil area deformation estimation algorithm based on thermal expansion hyperbolic model

3 試驗研究與分析

3.1 測區概況

洞庭湖周邊區域是典型的軟土地區。其中洞庭湖大橋位于岳陽市北門渡口下游,是湘北干線上跨越洞庭湖口的一座特大型橋梁,于1996年12月19日動工興建,2000年12月26日通車,運營總長為5 784.5 m,是連接君山區和岳陽樓區的城市主干快速路[24]。試驗中采用的TerraSAR-X影像空間覆蓋范圍和研究區域范圍如圖2所示。岳陽是一個擁有580萬人口的城市,居住區眾多,交通網絡發達。每逢汛期,洞庭湖周邊居民區都存在著堆積變形、地表塌陷甚至邊坡滑坡的安全隱患。因此,對該地區進行長期的時空變形監測,對預防交通安全隱患具有重要意義。

圖2 研究區概況Fig.2 Study area overview

3.2 試驗處理流程

試驗收集了2011年12月至2013年4月的24幅TerraSAR-X雷達衛星遙感影像。影像的前期干涉處理利用ENVI SARScape 5.2的雷達干涉處理模塊,外部高程數據為30 m分辨率的SRTM-DEM(shuttle radar topography mission digital elevation model)數據。為了保證影像原始分辨率,將距離向和方位向的多視比設置為1∶1。干涉組合采用傳統PSInSAR技術的單一主影像干涉模式,主影像選擇2012年9月17日的影像,參數影像具體見表1。前期干涉處理的時空基線的閾值分別設置為130 m和300 d,通過平地相位的去除、濾波、去地形相位、軌道精化等處理后,共生成23幅差分干涉圖。然后利用Matlab編程實現PS點選取、PS基線網絡布設、時序InSAR建模及時空參數估計。在PS點識別時,采用三重指數閾值串行的方法,即振幅離差閾值、相干性閾值和強度值三重閾值均滿足要求的點為候選PS點。在實際處理時,采用基于振幅離散指數、相干系數的平均值和強度值的三重指數法來挑選高相干點,分別設置三重指數閾值為0.46、0.80和0.68[25-27],同時手動去掉了河流里少數由采砂船等造成的高強度值的PS候選點。在PS基線網絡布設時采用Delaunay三角網方法。圖3為篩選前后的對比圖。對每一條PS基線開展時間維參數估計,利用式(13)進行基線邊質量評估,設置時間相關系數閾值為0.8,并迭代刪除質量差的基線和孤立PS候選點。如圖3所示,經過篩選后,測區共有4420個PS點和9224條基線。對整個PS網絡進行空間維參數估計后,可獲取PS點上的各參數估值,反演得到低通形變分量,最后對殘余相位進行時間維高通濾波和空間維低通濾波,以獲取殘余相位中的高通形變分量,與低通形變分量累加后即可獲取高相干點上時序形變結果。

表1 所用SAR影像的基本參數

圖3 PS點與PS基線迭代篩選前后對比Fig.3 Comparison of PS points and PS baselines before and after the iterative screening

3.3 試驗結果與分析

圖4(a)顯示了測區內熱膨脹系數Th的分布結果,圖4(b)為PS點上熱膨脹系數的定量分布統計結果?？梢园l現熱膨脹效應對于此測區地表形變的影響非常明顯。通過統計分析,58%的PS點其熱膨脹系數值分布在[-0.1, 0.1] mm/℃內,86%的PS點在[-0.2, 0.2] mm/℃內,95%的PS點在[-0.3, 0.3] mm/℃內。從顏色空間分布可以看出,熱膨脹系數值在不同區域存在顯著差異,從北至南呈現數值降低的趨勢,且數值有正有負。這一現象的原因是在熱膨脹特性參數恒定的情況下(當不同的基礎設施由類似的材料制成時,其熱膨脹特性參數視為恒定),熱膨脹系數值與建筑物的高度成正比[28-29]。通過研究區域歷史谷歌影像可知,在2012—2013年期間,研究區域由北向南建筑物高度整體呈由高向低的趨勢,與本文獲取的熱膨脹系數分布整體保持一致。數值有正有負的原因如圖5所示,熱膨脹導致的監測對象上某個點的總形變是水平形變和垂直形變的矢量之和,當這個點沿視線向遠離SAR傳感器時,該點的水平形變表現為負值。因此當建筑比較高時,其水平方向膨脹在雷達視線向分量較小,而豎向的分量則較大,使得熱膨脹引起的總形變為正[28-29]。為詳細分析熱膨脹系數分布的特征,本文分別在研究區域內的北部、中部和南部3個典型區域開展局部分析。以圖4中A、B、C幾個典型區域為例,由局部放大圖(圖4(c)、(d))可以看出,A、B區域都以高樓為主,因此其熱膨脹系數值均表現為正。其中A區內最大熱膨脹系數達0.31 mm/℃,B區的熱膨脹系數為0.1 mm/℃。而C區的情況則相反(4(e))。C區主要以學校、老住宅區、長條型道路為主,建筑物高度相對較低,水平方向膨脹在雷達視線向負向的分量,遠遠大于其豎直向的熱膨脹抬升,使得到的熱膨脹參值以負值為主,熱膨脹系數最小為-0.23 mm/℃。

圖4 熱膨脹參數分布Fig.4 Distribution of thermal expansion parameters

圖5 單點熱膨脹形變的幾何圖示Fig.5 Geometric illustration of single point thermal expansion deformation

圖6(a)、圖6(b)分別為雙曲線形狀參數va和vb的結果。由圖6(a)—(b)可知,雙曲線參數va在整體空間分布上并無顯著趨勢,其原因主要是va為雙曲線形狀參數,與各PS點空間分布并不相關;而vb為直線的斜率,與各PS點的沉降量相關,因此在空間上有明顯的趨勢性。經統計發現,va、vb為負值的占比分別為78%和84%,表明整個測區的形變以沉降趨勢為主。圖6(c)為降水參數Pr分布圖。由圖6(c)可知,圖中靠近水域的1、2區域降水參數偏小,都在1.6×10-3左右,而處在內陸的3區域則達到了2.8×10-3,呈現出從水域附近到內陸逐漸增大的趨勢。圖6(d)為獲取的高程改正結果,主要集中分布在-10～10 m。

圖6 模型參數估計結果Fig.6 Results of model parameter estimation

需要特別指出的是,由于洞庭湖大橋并非修建在軟土地基之上,不能利用本文的熱膨脹雙曲線模型。因此,對大橋區域的點進行掩膜處理,將這一部分點采用周期模型作為基礎模型,替換式(7)中的雙曲線模型分量部分進行形變建模[30],如式(16)所示

(16)

圖7為獲取的研究區域2011年12月28日至2013年4月3日的時間序列形變結果。從空間分布上來看,沉降較為明顯的區域一般分布在洞庭湖東岸(D區域),由湖岸向內部城區沉降區域逐漸減少。洞庭湖東岸沿線一帶沉降較為明顯,最大累積沉降達13.3 mm,而城區中部區域沉降量較小,最大累積沉降僅為4.7 mm。從時序變化上看,研究區域沉降速率在時序上呈現出明顯的先快后慢的特點,符合軟土沉降的非線性規律。2011年12月28日至2012年11月11日,洞庭湖東岸沿線一帶沉降趨于穩定,而從2012年11月11日至2013年4月3日,沉降趨勢逐漸增大,其中最大沉降發生在2013年4月3日,局部(D區內)累積最大沉降超過15 mm。除此之外,測區右下角(E區域)從2013年1月5日開始逐漸下沉,直到2013年4月3日,最大沉降量達到12 mm。通過查找歷史資料,發現測區右下角在2013年前后開展了廠房、建筑物拆除等基建工作,工程沉降或許是導致此處沉降的主要原因。試驗中還探測到了少部分抬升區域,如東風湖沿湖區域(F區域)分布了部分抬升點,抬升量主要集中在2～6 mm。由于該區位于元江凹陷與阜山隆起之間的地質構造邊緣地帶[1],因此推斷該抬升與岳陽市所處的地理位置,江漢洞庭盆地的地質構造、阜山隆起的影響有關。

圖7 研究區的時序形變Fig.7 Time-series deformation of the study area

為了更清晰地顯示變形的時序演化過程,分別提取了兩個特征點開展時序形變分析。如圖8所示,PS1位于岳陽市巴陵廣場瞻岳門城樓樓頂東南角處,PS2位于岳陽市岳陽樓區建設北路與外灘一路交叉口以南50 m處。圖9為兩個特征點上各形變分量的時間序列演化情況。PS點上各形變分量的占比如表2所示?？梢钥闯?無論是PS1還是PS2,總形變均以雙曲線分量為主,占比分別為81.5%和78.6%;熱膨脹及降水相關形變分量均小于2 mm。PS1的累積沉降量達9.5 mm,發生在2013年1月16日,而PS2累積沉降達7.4 mm。

表2 PS點上各形變分量占比

圖8 PS點位置及實地照片Fig.8 PS points location and field photos

圖9 PS1、PS2各形變分量時序形變結果Fig.9 Time-series deformation of each component on PS1 and PS2

分別提取了洞庭湖大橋的3個特征點(見圖8中PS3、PS4和PS5),其時間序列形變如圖10所示。由于PS3與PS5位于河兩岸的橋墩附近,形變結果接近,相比于位于橋梁中間的PS4形變整體浮動偏大。由圖10中可以看出,在3個特征點上的總時序形變都呈現出明顯的周期性變化,其中周期形變分量都達到了80%以上,最大可達-13 mm,而熱膨脹分量形變都在±2.5 mm之內。

圖10 PS3、PS4和PS5各形變分量時序結果Fig.10 Time-series deformation of each component on PS3, PS4 and PS5

根據文獻[13],熱膨脹特性參數能夠反映觀測對象的物理性質,用簡單的函數估計橋梁材料的熱膨脹特性參數,其表達式為

K=Th/H

(17)

式中,K為熱膨脹特性參數;H表示橋梁主梁相對于橋墩底部的高度;熱膨脹系數Th可由本文算法獲取。根據本文的估算,橋梁材料的平均熱膨脹特性參數為12.1×10-6/℃。根據收集的現場施工設計資料,洞庭湖大橋主梁采用C50預應力混凝土連續肋板結構。該混凝土的熱膨脹特性參數約為6×10-6～13×10-6/℃。因此,本文的估算結果與橋梁材料的物理性質有很好的一致性[31]。對該橋3個特征點的觀測結果顯示,在觀測周期內,該橋累積最大沉降值為8 mm,發生在2012年8月26日。表明在整個觀測期間,整個橋梁基本處于穩定狀態。

為了揭示熱膨脹分量與環境溫度之間的相關性,分別對3個特征點上的熱膨脹形變分量開展了相關性分析。圖11為各特征點熱膨脹形變分量與環境氣溫的相關性分析結果。由圖11可以看出熱膨脹形變分量與環境氣溫之間具有明顯的相關性。通過計算可得,3個特征點與溫度的平均相關系數分別為0.92、0.93和0.91。

圖11 熱膨脹分量與環境氣溫Fig.11 The thermal expansion-related deformation component and the air temperature

表3 熱膨脹雙曲線模型參數的不確定度

圖12 模型參數敏感性分析結果Fig.12 Sensitivity analysis results of model parameters

4 精度分析與討論

由于研究區域沒有可用的外部形變測量數據,為了驗證算法的精度,分別利用殘余相位和時間相關系數評估本文方法的可靠性。根據文獻[33]所述,殘余相位的大小反映了模型與真實形變的擬合程度,即可以用來表征模型的建模精度。干涉圖中殘余相位越小,表征建模精度越高。因式(7)相位模型是基于每一條PS基線建立,需將每一基線邊的殘余相位恢復至每一個PS點上,再分別對23個干涉對中所有PS點的殘余相位進行標準差(standard deviation,STD)的計算,再與傳統線性速率模型估計的結果進行對比。對比結果如圖13所示。由圖13可以明顯看出,熱膨脹雙曲線模型的殘余相位明顯小于線性模型,說明熱膨脹雙曲線模型在建模本測區的軟土時序形變中更加適用。熱膨脹雙曲線模型在所有干涉對上的殘余相位均小于0.60 rad,STD為0.40 rad,而線性模型殘余相位的STD為0.63 rad,提升了36.5%。

圖13 各干涉圖的殘余相位對比Fig.13 Residual phase comparison of each interferogram

由文獻[34—35]可知,時間相關系數也可以作為模型建模精度的又一指標。對應PS點上的時間相關系數越高,則表明其建模精度越高。因此,本文計算了所有PS點在不同模型下的平均時間相關系數,結果如圖14所示,可以明顯看出熱膨脹雙曲線模型的平均時間相關系數(整個圖都是深黃色)整體顯著高于線性模型結果。熱膨脹雙曲線模型在所有PS點上的平均時間相干系數STD為0.92,而線性模型為0.63。

圖14 模型的平均時間相關系數對比Fig.14 Comparison of the average temporal coherence of the models

根據前述試驗結果及精度分析,現針對本文方法進行討論如下。

(1) 本文提出將雙曲線模型作為PSInSAR處理基礎模型,并融入熱膨脹和降水參數分量,取代原線性速率模型,在理論上顧及了軟土形變隨時間呈非線性變化的物理特征、湖區基礎設施的熱膨脹效應和外界環境降水的影響因素。相比于傳統線性速率模型,本文方法殘余相位的顯著提升也證實了方法的可靠性。本文方法可利用InSAR相位方程組估計出各PS點的熱膨脹系數值,由此可以計算出基礎設施的熱膨脹特性參數,提供了一種利用InSAR相位觀測值獲取材料熱膨脹參數的手段。

(2) 由于研究區域內存在著時空異質性,不同區域、不同時期軟黏土區分布面積,土質成分等特征并不完全相同。例如,洞庭湖大橋并不是全部修建在軟土地基之上,利用本文的模型建模仍存在局限性。因此在試驗處理時,將反映周期變化特征的周期模型作為基礎模型,用于洞庭湖大橋上的點的形變建模,以更好地建模橋梁隨熱膨脹變化的周期特性。試驗結果表明周期模型能更好地反映大橋這一特殊設施的實際形變隨時間演化的規律。而本文的雙曲線模型對建設在軟黏土區域、高程相對較低的基礎設施(如路基路面)則更有優勢。

(3) 試驗中獲取的數據集為2011年12月至2013年4月,由于過境時間的限制未跨越兩年的完整周期,且在數據集期間沒有相匹配的地面水準測量數據作為外部驗證數據。對于基礎設施的時序形變特征的研究,如果可以用長達兩年以上的數據進行分析則更具說服力。除利用殘余相位來驗證本文模型建模效果的提升之外,為彌補外部水準數據的缺失,搜集了洞庭湖區域早期沉降歷史資料以進行驗證。通過文獻[36]得知,在研究時段內,洞庭湖區沉降量在1 cm以內,這與本文研究得到的沉降量一致。

5 結 論

本文提出一種軟土區域時序InSAR形變與環境物理參數聯合估計方法。根據軟土區的變形特性,選擇雙曲線模型作為基礎形變模型,并引入熱膨脹和降水參數分量,改進了傳統的PSInSAR線性速率形變模型,并可同時獲取區域特征點的熱膨脹系數,提供了一種獲取材料熱膨脹參數的遙感手段。將洞庭湖區域基礎設施分布密集區作為試驗區域,設計了包含PS點選點、基線網絡構建、InSAR時序建模、時空維參數估計、時序總形變的生成這一完整算法流程。試驗利用2011年12月至2013年4月的24幅TerraSAR-X雷達衛星遙感影像,獲取了測區熱膨脹雙曲線模型的各參數值,并反演得到了測區最終的時序形變場。試驗結果表明,截至2013年4月3日,洞庭湖東岸沿線一帶沉降較為明顯,最大累積沉降達13.3 mm,而城區中部區域沉降量較小,最大累積沉降僅為4.7 mm。根據生成的熱膨脹系數圖,估計了洞庭湖大橋的熱膨脹特性參數值,與橋梁混凝土的物理性質一致。為了驗證結果的可靠性,分析了23個干涉對的殘余相位,得出本文方法的殘余相位相比傳統線性速率模型提升了36.5%。通過搜集洞庭湖區域早期沉降歷史資料得知,在研究時段內,洞庭湖城區沉降量總體分布在1 cm以內,與本文研究得到的沉降量級一致。