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

1. 長沙理工大學(xué)交通運(yùn)輸工程學(xué)院,湖南 長沙 410114; 2. 中南大學(xué)地球科學(xué)與信息物理學(xué)院,湖南 長沙 410083; 3. 洞庭湖生態(tài)環(huán)境遙感監(jiān)測湖南省重點(diǎn)實(shí)驗(yàn)室,湖南 長沙 410114

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

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

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

針對上述問題,本文提出一種基于熱膨脹雙曲線模型的時序InSAR形變與環(huán)境物理參數(shù)聯(lián)合估計(jì)方法。該方法在PS基線網(wǎng)絡(luò)的相位建模環(huán)節(jié)以軟土區(qū)域形變的先驗(yàn)?zāi)Ｐ汀p曲線模型為基礎(chǔ)模型,將熱膨脹系數(shù)和降水參數(shù)融入其中,即分別建立受熱膨脹和外界環(huán)境降水影響的相位分量模型。用上述模型取代原線性速率模型,解算時不但可以獲取每個PS點(diǎn)的變形值,還可同時獲取模型中相應(yīng)的環(huán)境物理參數(shù),實(shí)現(xiàn)洞庭湖區(qū)域基礎(chǔ)設(shè)施總形變序列與環(huán)境物理參數(shù)的聯(lián)合反演。獲得的形變序列與環(huán)境物理參數(shù)可為基礎(chǔ)設(shè)施的安全穩(wěn)定分析和災(zāi)害防治提供可靠的數(shù)據(jù)支撐。

1 時序InSAR模型構(gòu)建

1.1 PSInSAR傳統(tǒng)線性速率模型

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

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

(1)

(2)

1.2 熱膨脹雙曲線模型

(3)

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

(4)

(5)

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

(6)

(7)

(8)

(9)

式中,Tempm為研究區(qū)域在SAR影像過境時刻的環(huán)境溫度;Precm為研究區(qū)域當(dāng)月降水量;Th和Pr分別為熱膨脹系數(shù)和降水參數(shù),為待求模型未知參數(shù)。

(10)

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

(11)

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

2 模型的未知參數(shù)估計(jì)

2.1 基于LAMDBA-SVD的時間維參數(shù)估計(jì)

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

(12)

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

(13)

2.2 基于雅可比迭代的空間維參數(shù)估計(jì)

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

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

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

X=B0X+f

(14)

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

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

(15)

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

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

3 試驗(yàn)研究與分析

3.1 測區(qū)概況

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

圖2 研究區(qū)概況Fig.2 Study area overview

3.2 試驗(yàn)處理流程

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

表1 所用SAR影像的基本參數(shù)

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

3.3 試驗(yàn)結(jié)果與分析

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

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

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

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

圖6 模型參數(shù)估計(jì)結(jié)果Fig.6 Results of model parameter estimation

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

(16)

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

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

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

表2 PS點(diǎn)上各形變分量占比

圖8 PS點(diǎn)位置及實(shí)地照片F(xiàn)ig.8 PS points location and field photos

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

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

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

根據(jù)文獻(xiàn)[13],熱膨脹特性參數(shù)能夠反映觀測對象的物理性質(zhì),用簡單的函數(shù)估計(jì)橋梁材料的熱膨脹特性參數(shù),其表達(dá)式為

K=Th/H

(17)

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

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

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

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

圖12 模型參數(shù)敏感性分析結(jié)果Fig.12 Sensitivity analysis results of model parameters

4 精度分析與討論

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

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

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

圖14 模型的平均時間相關(guān)系數(shù)對比Fig.14 Comparison of the average temporal coherence of the models

根據(jù)前述試驗(yàn)結(jié)果及精度分析,現(xiàn)針對本文方法進(jìn)行討論如下。

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

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

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

5 結(jié) 論

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