基于雷達衛(wèi)星時序分析技術的荊江沿岸堤防形變研究

羅天文

(1.武漢大學 遙感信息工程學院， 武漢 430079；2.長江科學院 空間信息技術應用研究所，武漢 430010;3.貴州省水利水電勘測設計研究院, 貴陽 550002)

1 研究背景

近幾十年來，雷達衛(wèi)星時序分析技術已應用于分析各類與不同現(xiàn)象相關的地面變形[1]。基于雷達衛(wèi)星的時序分析技術可以區(qū)分為2大類：①永久散射體干涉測量(Persistent Scatterer-Interferometric Synthetic Aperture Radar,PS-InSAR)，基于單個主圖像的使用以及時間序列中穩(wěn)定和高度相干散射體的選擇[2-3]。②小基線子集(Small Baseline Subset,SBAS)干涉測量，其組合多個計算的干涉圖，選擇以軌道(基線)之間的小空間間隔為特征對，以限制空間去相關效應[4-5]。Berardino等[6]、Ferretti等[7]提出原始PS-InSAR方法，隨后幾年又提出了幾種改進方法。這些方法已被用于歐洲遙感(ERS)和環(huán)境衛(wèi)星(Environmental Satellite，Envisat)數(shù)據(jù)的分析，證明永久散射體干涉測量(PS-InSAR)是相對低成本高精度的陸地變形監(jiān)測的強大工具。該技術在大范圍區(qū)域高精度形變監(jiān)測領域的巨大優(yōu)勢，使得大量國內(nèi)學者參與到雷達衛(wèi)星時序分析技術的研究中，也取得了飛速的進展與卓越成就。2013年裴媛媛等[8]利用Envisat ASAR影像得到了上海市長江口南岸和杭州灣北側堤壩的沉降速率；2019年李陶等[9]利用COSMO-SkyMed數(shù)據(jù)對公明水庫路堤邊坡的沉降過程進行了監(jiān)測，都取得了很好的效果。得益于Sentinel-1A(歐洲航天局哥白尼計劃中載有C波段合成孔徑雷達的地球觀測衛(wèi)星)的開放存取數(shù)據(jù)，每12 d提供一次10 m分辨率的全球觀測雷達數(shù)據(jù)，可以利用PS-InSAR技術高效快捷地對地表進行大面積形變監(jiān)測。本文展示了如何通過永久散射干涉測量法和Sentinel-1A數(shù)據(jù)監(jiān)測荊江大堤的形變，對雷達圖像進行處理，提取出持久散射體，并對圖像處理得到的數(shù)據(jù)進行分析建模，獲得荊江大堤形變情況。

2 PS-InSAR監(jiān)測形變算法

保持較高時間相干性的相干點目標為永久散射體PS(Permanent Scatterer)，PS-InSAR形變監(jiān)測的基本思路是構建一系列基于相同主影像的累計差分干涉圖，然后分析時序影像的幅度與相位變化，識別選取保持較高時間相干性的永久散射體，進而建立基于永久散射體的時序差分干涉模型，根據(jù)PS點的相位特征與點之間的關系，分離出有效相位，并轉(zhuǎn)換為堤防形變數(shù)據(jù)。

首先在研究區(qū)域獲取的雷達衛(wèi)星時序數(shù)據(jù)中選擇一幅合適的雷達圖像作為主圖像，然后將其他圖像作為輔圖像生成干涉圖。再利用DEM(Digital Elevation Model)去除地形相位，從而得到差分干涉圖。差分干涉相位Δφi可以表示為

Δφi=φtopo+φdef+φatm+φnoi。

(1)

式中：φdef為地表形變相位，φdef由線性形變相位Δφlin(x,y)和非線性形變相位Δφnon-lin(x,y)相加而成；φtopo為DEM誤差相位；φatm為大氣延遲相位；φnoi為噪聲相位(散射體變化，熱噪聲，誤配準造成的誤差)。在雷達時序數(shù)據(jù)中選擇振幅和相位穩(wěn)定的點，作為候選PS點，通過設置一定的閾值，將候選PS點連成空間上相互連通的稀疏網(wǎng)格，其中2相鄰PS點x和y的相位差Δφdiff(x,y)可以表示為

Δφdiff(x,y)=Δφtopo(x,y)+Δφlin(x,y)+

Δφnon-lin(x,y)+Δφatm(x,y)+Δφnoi(x,y)。(2)

運用干涉獲取形變最重要的是獲取形變相位，對于相鄰的PS點來說，Δφnon-lin(x,y)很小，近似為噪聲忽略不計，Δφlin(x,y)為有效形變相位，式(2)可以轉(zhuǎn)換為

Δφatm(x,y)+Δφnoi(x,y)=

Δφdiff(x,y)-Δφtopo(x,y)-Δφlin(x,y)。

(3)

令Δω(x,y)=Δφatm(x,y)+Δφnoi(x,y)，則有

Δω(x,y)=Δφdiff(x,y)-(Δφtopo(x,y)+

Δφlin(x,y)) 。

(4)

求解目標函數(shù)式(4)可獲得有效形變相位和DEM誤差相位。選取某一PS點候選點作為參考點，在所有時序差分干涉圖中迭代搜索，當Δφtopo(x,y)和Δφlin(x,y)滿足所有干涉圖的整體相干性最高時，則可分別作為DEM誤差相位與線性形變相位的最優(yōu)估計。定義PS點的相干性指標為

(5)

式中：M為干涉圖的個數(shù);xi和yi表示第i幅干涉圖中相鄰的2個候選PS點；j為復數(shù)的虛數(shù)單位;γs的取值范圍在[0，1]，候選PS點在所有干涉圖上的相位離散度越低，γs越接近1，因此可通過最大化γs來求得線性形變速率和DEM誤差。

獲得DEM誤差相位最優(yōu)估計后，再用差分干涉相位減去大氣相位，利用時序差分干涉模型，對所有的像素點進行分析，當像素點的γs值滿足設定的閾值時則可最終確定為PS點，然后重新構建稀疏網(wǎng)格，對PS點進行相位解纏，利用測量平差方法和最小二乘平差法估計所有 PS 點的線性形變速率從而獲得堤防的形變量。

3 試驗數(shù)據(jù)與方法

3.1 研究區(qū)域

本文選擇的研究區(qū)域為下荊江石首至監(jiān)利段荊江大堤。荊江位于江漢平原的南麓,是長江在湖北境內(nèi)的部分河段，上起枝城,下迄城陵磯,全長約400 km，它是世界上典型的自由河曲之一[10]，河道蜿蜒曲折,特別是下荊江段。荊江大堤被列為長江防洪重點確保堤，保護了江漢平原大量城鎮(zhèn)，被稱為江漢平原的“生命之堤”，大堤的形變監(jiān)測工作非常必要。

3.2 試驗數(shù)據(jù)

本文利用Sentinel-1A數(shù)據(jù)進行堤防形變研究，Sentinel-1A是一顆C波段雷達衛(wèi)星，于2014年4月3日發(fā)射，具有12 d重復周期，空間分辨率為10 m，中心頻率為5.405 GHz。數(shù)據(jù)由歐洲航天局免費提供(https:∥scihub.copernicus.eu/dhus/)。Sentinel-1A有4種模式:條帶模式(Strip Map，SM)、干涉寬幅模式(Interferometric Wide，IW)、極寬幅模式(Extra-Wide Swath，EW)和波譜模式(Wave)。本研究利用單視斜距復數(shù)影像(Sentinel-1A IW Level 1 (L1) SLC)產(chǎn)品。收集了2018-03-26—2019-04-02觀測到的Sentinel-1A IW模式(包含VV極化方式和VH極化方式)數(shù)據(jù)共32景，數(shù)據(jù)參數(shù)如表1所示。

表1 本研究中Sentinel-1A的數(shù)據(jù)參數(shù)

3.3 試驗方法

PS-InSAR處理技術分為5個階段，具體步驟如下。

(1)配準和生成干涉圖。實現(xiàn)PS-InSAR算法首先在研究區(qū)域獲取的32幅SAR影像中選擇一幅SAR影像作為主圖像，然后將其他圖像作為輔圖像生成31幅干涉圖，于2018年10月4日獲得的Sentinel-1A圖像具有最合適的垂直基線和一致性。因此，該圖像被認為是創(chuàng)建干涉圖堆棧的主圖像，圖1為干涉相對基線分布情況。

圖1 干涉相對基線分布

(2)生成差分干涉圖。在生成干涉圖之前輔影像需要被配準到主影像。首先對圖像進行粗配準，然后對高分辨率圖像和相對于主圖像的亞像素圖像進行精配準。2個配準SAR圖像的組合生成干涉圖。干涉圖是輔影像和主影像各值的乘積，每個像素的輸出結果是新因子的關聯(lián)。干涉圖采用二維相位差圖的形式，在本文中設置相關性閾值為0.75。在生成干涉圖后，本文中利用空間分辨率為30 m的SRTM version-4數(shù)據(jù)從生成的干涉圖中去除地形相位，得到?jīng)]有地形相位影響的差分干涉圖。

(3)選擇永久散射體。得到干涉圖以后需要選擇候選的持久散射點，候選持久散射點的選擇一般分為2個步驟。第1步根據(jù)后向散射振幅選擇初始的持續(xù)散射點;第2步根據(jù)校正后的相位選擇主要的持續(xù)散射點，這些點上的相位在時間上應該是穩(wěn)定的。在PS點的選擇上，由于地形、軌道和地面變形等因素沒有考慮干涉相位，因而不能利用干涉相位的相干信息來選擇持久散射點，所以采用基于振幅的方法來選擇初始點。傳統(tǒng)的基于振幅的PS點選擇方法有基于振幅色散指數(shù)的振幅色散法[11]和信雜波比法。本研究采用振幅色散指數(shù)，理論上來說隨著時間的推移，相位變化越小的像素，振幅色散指數(shù)也越低。通過在振幅色散指數(shù)上定義特定的閾值限制來選擇初始候選像素。

(4)相位解纏。2幅SAR復圖像獲得的干涉相位差值是被周期折疊后位于(-π，π)之間的相位主值，想要得到正確的形變信息，必須進行相位解纏,恢復被模糊掉的相位周期，獲得目標的真實相位差。關于相位解纏算法國內(nèi)外進行了深入的研究，如區(qū)域增長法[12]、最小費用流法[13]、最小二乘法[14]等。在本研究中，采用最小費用流法進行相位解纏。該方法適用于由于存在大面積低相干性區(qū)域或因生長限制因素而導致解纏困難的情況，在這種情況下相對于其他方法最小費用流法的性能更好。

(5)地理編碼。最后將展開相位轉(zhuǎn)換為相應的形變并進行地理編碼。地理編碼的結果導入ArcGIS?10.5中，對PS點進行適當?shù)牟逯?，以便更好地展示PS-InSAR結果。通過對研究區(qū)域的地質(zhì)統(tǒng)計，并對基于點空間分布的數(shù)據(jù)進行正態(tài)性檢驗，選擇逆距離加權(IDW)算法，利用已知PS點單元計算未知單元，生成柵格圖，作為荊江堤防的位移模型。在插值圖制作完成后，利用Sentinel-1A圖像處理生成的位移值，生成時間序列形變圖，以便更好地比較位移隨時間的變化。

4 試驗結果

利用PS-InSAR技術計算的形變速率結果如圖2所示。 對試驗結果進行分析， 從圖2的結果中可以看出， 沿荊江堤防主干區(qū)域出現(xiàn)不均勻的緩慢形變現(xiàn)象(17～30 mm/a)， 形變導致堤頂高程相對變化， 堤防的形變速率在不同區(qū)域表現(xiàn)較大差異， 最大形變可以達到30 mm/a。 形變在軟土基礎上的堤防段表現(xiàn)較為嚴重， 而混凝土基礎的堤防形變程度則較輕。 荊江堤防有分布廣、 岸線長的特征， 不同區(qū)域的土質(zhì)不同， 形變情況也不同， 而不均勻形變會對堤防的安全性將構成危險。 基于時序SAR影像監(jiān)測形變的技術可以高精度、 大面積地監(jiān)測堤防形變， 為堤壩監(jiān)管和防汛抗風暴潮提供了有效信息。

圖2 最終形變速率結果

為了更精細地監(jiān)測形變結果，分析在1 a時間里堤防的形變變化，本文在目標沿線選擇采樣點進行時間序列的形變分析，采樣點選擇的依據(jù)為：采樣點必須在荊江堤防沿線，并且周圍的PS點數(shù)量足夠，具體表現(xiàn)為該點100 m范圍內(nèi)至少存在5個PS點，通過鄰域分析可減小形變速率誤差。這里選取了8個采樣點進行分析，分別包含軟土為基礎的堤防和混凝土為基礎的堤防，各類采樣點的分布情況見圖3。雷達數(shù)據(jù)時間序列跨度為2018-03-26—2019-04-02，平均每12 d有一次時間序列的結果，由此得到的結果如圖4所示。

圖3 采樣點分布(A-H共8個點)

圖4 采樣點的時間序列形變

從圖4可以發(fā)現(xiàn)，所有采樣點都呈現(xiàn)出了比較明顯的共同規(guī)律：①2018年9月前后，形變速率變化規(guī)律明顯不同，從明顯的地表沉降轉(zhuǎn)變?yōu)榈乇淼奶?；②形變速率?018年9月左右達到最大，在2018年12月左右達到最大抬升。分析原因是江漢平原地區(qū)上半年雨水豐富，特別在7月份的梅雨季節(jié)，荊州地區(qū)經(jīng)歷了持續(xù)降雨，根據(jù)國家氣象中心的數(shù)據(jù)，降水量達到了100～180 mm，持續(xù)的降水造成了堤防的沉降，而到了秋冬季雨水減少、氣溫降低、土壤凝結，從而出現(xiàn)了抬升的現(xiàn)象，到了1月份開始出現(xiàn)降雪，雨水增多，又出現(xiàn)沉降。

5 結 語

由于越來越多地使用高分辨率雷達圖像，近年來在基礎設施空間監(jiān)測方面取得了重大進展。利用雷達衛(wèi)星時序分析技術對堤防的形變進行監(jiān)測，其消耗成本遠比地面測量手段的消耗成本低。本文利用Sentinel-1A數(shù)據(jù)得到了石首至監(jiān)利段在2018年4月至2019年4月間的平均沉降速率，結果表明，堤防在不同的時間，其形變速率不同，特別在雨季和旱季形變情況尤為不同。試驗證明利用雷達衛(wèi)星時序分析技術可實現(xiàn)荊江堤防的沉降情況高效大范圍監(jiān)測，表明該技術在長江堤防形變監(jiān)測上具有很大的應用潛力。