基于小基線集干涉測量的天津地區地表沉降時空分析

朱子林, 任 超,2, 周 呂,3,4*, 施顯健, 李現廣

(1.桂林理工大學測繪地理信息學院，桂林 541004；2.廣西空間信息與測繪重點實驗室，桂林 541004；3.武漢大學地球空間環境與大地測量教育部重點實驗室，武漢 430079；4.城市空間信息工程北京市重點實驗室，北京 100038)

城市地表沉降是當今世界各國城市化發展進程中廣泛存在的一種緩慢的區域性地質災害，該現象是由自然地質災害(如地震、火山等)或人類活動(如地下水開采、建筑荷載等)引起的，城市建設、農業發展、防洪排澇等方面帶來嚴重影響[1-3]，嚴重影響著人民的生命財產安全以及社會的可持續發展。天津市大部分地區由于大量的地下水開采以及工業和人口的迅速增長，從而導致地面發生嚴重沉降，監測整個天津地區沉降過程，保證人民群眾的人身和財產安全是亟待完成的工作[4]。目前，對于地表形變監測和分析，傳統的方法主要包括水準測量、GNSS(global navigation satellite system)測量等可以得到研究區域內精度較高以及高時間分辨率的形變量[5],但是這些方法的監測成本較高，監測點在整個研究范圍內密度有限，而且其空間分辨率較低，很難有效地分析整個研究區域范圍內的形變。因此，對城市地表形變進行大規模的觀測，分析更為復雜的地表沉降現象是當前研究的熱點問題。常規的合成孔徑雷達差分干涉測量技術(differential interferometric synthetic aperture radar，D-InSAR)因其在地表沉降監測方面具有較高的空間分辨率，逐漸成為國際上公認的開展地表沉降監測的有效方法[6-7]。但是因為空間失相干和大氣延遲的影響，傳統的D-InSAR技術在應用過程中精度和可靠性會降低，為了減小上述影響，小基線集合成孔徑雷達干涉測量(small baseline subset interferometric synthetic aperture radar,SBAS-InSAR)技術就是在常規的D-InSAR基礎上新興的一種技術，該技術不僅限制了時空失相干，提高影像的利用率，同時也豐富了影像的觀測信息[8]。姜德才等[9]已經利用MCTSB-InSAR(multiple images coherent targets small baselines interferometric synthetic aperture radar)及其相關技術監測天津市地鐵沉降的時序變化，并驗證天津市地鐵沿線沉降比較嚴重的區域；白澤朝等[10]基于Sentinel-1A數據利用PS-InSAR(permanent scatters interferometric synthetic aperture radar)技術監測天津市地表形變并確定天津市地表形變是地下水過度開采造成的；雷坤超等[11]利用PS-InSAR技術研究天津市地表沉降并確定其沉降區域有向周邊延伸的趨勢，但同時利用SBAS-InSAR技術和Sentinel-1A數據研究天津地區地表沉降的工作尚不多見。因此，利用2017—2019年的30景Sentinel-1A 數據，采用SBAS-InSAR技術監測天津市地表沉降，獲取了天津地區的沉降分布情況以及沉降形變速率，通過與相同時間段內天津市二等水準測量數據對比分析，驗證了SBAS-InSAR技術監測城市地表沉降的可行性，主要討論分析了天津地區的地表沉降時空變化與地下水和城市化發展的關系，以期為掌握天津地區地表形變提供理論支撐。

1 SBAS-InSAR技術原理

SBAS-InSAR是由Bernardino等[12]和Lanari等[13]于2002年正式提出的一種時序InSAR的分析方法。此方法主要是將多景SAR影像進行一系列的排列組合，得到相關空間基線差分干涉圖，同時可以連接被較大空間基線分開的SAR數據集，確保生成的像對中沒有孤立的集群，以增加觀測數據的時間采樣率，可以更好地克服空間失相干現象。SBAS-InSAR技術主要使用奇異值分解法(singular value decomposition,SVD)來求解研究范圍內的形變速率，進而獲得高相干點的沉降形變速率和時間序列。其主要步驟如下[14-15]。

(1)在(t0,t1,…,tn)時間內，對獲取的N+1景同一研究區域的Sentinel-1A影像數據，選擇其中任意一幅影像作為超主影像，并將其他輔影像配準到超主影像上。假設每景影像都至少有一景影像與之干涉，則N+1景影像生成M景差分干涉圖，M滿足條件：

(1)

(2)對于任意差分干涉圖j在輔影像tA時刻和主影像tB(tB>tA)時間內獲得的影像進行干涉處理，方位向坐標x和距離向坐標r像素的干涉相位(忽略大氣延遲相位、殘余地形相位和噪聲相位的影響)可以表示為

δφj=φB(x,r)-φA(x,r)≈

(2)

式(2)中：φ為干涉相位；j為差分干涉圖的景號，j∈(1,2,…,M)；λ為雷達信號波長；φA(x,r)和φB(x,r)分別為tA和tB時刻SAR影像的相位值；d(tA,x,r)和d(tB,x,r)分別為tA時刻和tB時刻相對于d(t0,x,r)=0的視線方向上的累積形變量。

(3)為了更加清晰地表示地表沉降的時間序列，可以將干涉相位值用兩個采集時間內的平均相位速率vj和該時間段之間的乘積來表示，即

(3)

則第j幅干涉圖的相位值可以表示為

(4)

即在主影像時間間隔和輔影像時間間隔上每個時間周期上速度的積分，表達成矩陣的形式,即

Bv=δφ

(5)

在計算系數矩陣B(m×n)的過程中，因為SBAS-InSAR技術是采用多主影像時空基線的方法獲取差分干涉圖，所以矩陣B可能會出現秩虧的情況。矩陣B的廣義逆矩陣可以用奇異值分解法(SVD)求解，從而得到速度矢量的最小范數解，最終通過對每個時間段的速度進行積分得到每個時間段內的形變量。

2 研究區概況及數據源

天津市是中國北方的一個沿海城市，部分地區被渤海環繞，該城市中心坐標約為117°20′E和39°12′N。由于大量的地下水開采、工業和人口的快速增長，天津大部分城市地區遭受了嚴重的沉降。水準歷史資料顯示，自20世紀80年代以來最大累積沉降已經超過1 m，自天津市開始實施限水政策控制后，沉降速率逐漸變慢，但是天津市中心城區以及周邊城區仍然缺乏控制和沉降監測[4,16]。近年來，天津地區城市與軌道交通建設飛速發展，現在已基本形成完整的交通網絡，新的沉降中心也在不斷產生。為了進一步研究天津地區地表沉降情況，選擇天津地區升軌模式下Sentinel-1A數據，由于原始影像數據量較大，后期工作量較大，處理起來比較困難，因此提前對原始數據進行裁剪，裁剪后的研究區范圍如圖1所示。

圖1 研究區地理范圍Fig.1 Geographical scope of the study area

研究選擇了覆蓋天津地區的30景升軌模式下Sentinel-1A影像數據[均為垂直發射水平接收(VH)極化方式]，時間跨度為2017年1月15日—2019年7月16日由歐洲航天局(European Space Agency，ESA)提供。精密軌道數據使用的是歐洲航天局提供的衛星精密軌道數據POD(precise orbit ephemerides)，用來提高配準和基線估算的精度[17]。此外，數據處理過程中所采用的DEM數據來自美國國家航空局(national aeronautics and space administration,NASA)提供的SRTM V4 DEM，分辨率為90 m。

3 數據處理與結果分析

3.1 數據處理

使用SARscape數據處理軟件和30景Sentinel-1A影像數據監測天津市2017—2019年的地表沉降情況，主要處理流程有：①數據預處理：由于原始數據量較大，處理起來比較困難，因此將30景原始Sentinel-1A影像數據進行預處理，并通過裁剪得到本文研究區域；②生成連接像對：選取2018年5月10日獲取的影像為超主影像，將其他影像配準到超主影像上，選取時間基線閾值為400 d和空間基線閾值為400 m，共生成139對小基線差分干涉圖集，其相對空間基線連接圖如圖2所示；③干涉工作及相位解纏：將139對干涉像對進行干涉處理，主要包括生成干涉圖，干涉圖的去平、濾波和相位解纏以及相干性計算，將所有的數據對都配準到超主影像上；④軌道精煉和重去平：該步驟主要用于估算和去除不符合要求的殘余恒定相位以及解纏后還存在的坡道相位；⑤生成時間序列結果：利用奇異值分解法對解纏后的相位進行解算，得到每個時段內的形變速率，最后通過地理編碼得到每個時間段在視線方向(LOS)的形變量(負值表示沉降量，正值表示抬升量)。

綠點代表SAR影像；線段代表干涉圖；黃點代表超主影像圖2 SAR干涉圖時空基線分布Fig.2 Temporal and spatial baseline distribution of SAR interferogram

3.2 結果分析

3.2.1 研究區域沉降速率分析

圖3為天津市2017—2019年地表沉降形變速率。從圖3可以看出：研究時間段內80%的研究區域沉降形變速率在-15.8～4.1 mm/a(正值表示地表上升，負值表示地表下降)。天津中心城區的沉降范圍為-5.9～4.1 mm/a，表明天津市中心城區地表相對穩定，這與天津市在20世紀80年代以來施行市內六區嚴格控制地下水的開采政策有關，因此天津地區市內六區的地表小沉降相對較小且比較穩定；西青區以及勝芳鎮(河北省霸州市所屬)出現明顯的沉降漏斗，其中王慶坨鎮沉降中心較為嚴重，沉降中心的工業區沉降最為明顯，年平均沉降速率達到-134.7 mm/a，最大累積沉降達到 324.6 mm，工業區以及居民小區都集中在該地區，工業用水以及生活用水都需要地下水，地下水開采嚴重以及地面負荷的增加是王慶坨鎮沉降漏斗形成的主要原因；勝芳鎮沉降漏斗位于河北省霸州市與西青區臨近，該區域的郝家堡村沉降較為嚴重，年平均沉降速率達到-114.9 mm/a，最大累積沉降量達到230.8 mm，工業化的發展以及城市居民增加導致該地區地下水的抽取較為嚴重，過度的地下水開采是導致該地區出現沉降漏斗的主要原因；研究區域內靜海縣沉降區中團泊新城沉降相對比較明顯，該區域年平均沉降速率達到-85.2 mm/a，最大累積沉降量達到167.4 mm，團泊鎮被列為第三批國家新型城鎮化綜合試點地區，近年來城鎮化發展較為迅速，新城突出“水和生態”逐漸成為多功能于一體的中國北方水上旅游城市，城鎮化的迅速發展，城市擴張以及城鎮化建設是該地區地表沉降的主要原因。

圖3 天津市2017—2019年地表沉降形變速率Fig.3 Deformation rate of surface subsidence in Tianjin from 2017 to 2019

3.2.2 SBAS-InSAR結果與水準數據對比分析

研究區域內主要沉降中心位于王慶坨鎮、勝芳鎮以及靜海縣，與Li等[16]利用SBAS-InSAR技術獲取的結果時空分布以及最大沉降區域(王慶坨鎮和勝芳鎮)基本吻合，驗證了本文方法獲取天津地區地表沉降形變的可靠性。為了進一步驗證基于Sentinel-1A數據采用SBAS-InSAR準確性，選取研究區域內相同時間段內15個二等水準點數據(水準點 1～15分布如圖4所示)，將二等水準測量數據與SBAS-InSAR相同時間段內計算結果進行比對，結果如圖5所示。從圖5可以看出，相同時間段內SBAS-InSAR計算的結果與二等水準數據結果吻合較好。對比分析兩種數據結果可知，兩種方法的最大相對誤差為8.1 mm，最小相對誤差為-1.54 mm，中誤差為2.9 6 mm，結果表明兩種方法得到的結果具有較好的一致性。

1～15為水準點位置圖4 天津地區二等水準分布情況Fig.4 Second-class level distribution in Tianjin

圖5 水準測量與SBAS-InSAR處理結果對比Fig.5 Comparison between leveling and SBAS-InSAR processing results

3.2.3 研究區地表沉降時空變化分析

利用時序特征點(特征點分布如圖6所示)分析天津地區地表沉降能夠較為直觀地看出各特征點以及重點區域的沉降時空變化過程，選擇四處不同沉降區域進行時序特征分析，結果如下：總體來看天津地區地表沉降不均勻較為明顯(圖7、圖8)，天津市中心城區沉降相對穩定，王慶坨鎮和勝芳鎮出現明顯的沉降漏斗，兩地的最大累積沉降量分別達到324.6和230.8 mm，與這兩個明顯的沉降漏斗相比，靜海縣的沉降相對平緩，最大累積沉降為167.4 mm。由圖7、圖8(a)可以看出，王慶坨鎮工業區沉降較為明顯，出現明顯的沉降漏斗，特征點沉降量成逐年增加趨勢，研究時段內該區域的最大累積沉降達到324.6 mm，工業化發展以及地下水過度開采是導致沉降的主要原因；研究時段內勝芳鎮沉降區的沉降隨時間變化逐漸增加，其中郝家堡村沉降最為明顯，最大累積沉降達到230.8 mm，該地區工業化發展較為迅速，工業和生活用水基本來自地下水的開采，過度開采地下水和工業的快速發展是該地區沉降的主要原因；靜海縣的團泊鎮相較于王慶坨鎮和郝家堡沉降相對較為平緩，最大累積沉降量為167.4 mm，團泊新城的發展迅速，城市化擴張是該地區沉降的主要原因；相較于三個沉降區明顯的沉降漏斗，天津市中心城區沉降相對穩定，由圖8(d)可以看出研究時段內沉降量隨時間的變化有明顯的起伏現象，中心城區發展相對較為穩定，同時中心城區施行限水政策嚴格控制地下水的開采，因此中心城區的沉降相對較為穩定。

1～15為時序特征點位置圖6 時序特征點分布情況Fig.6 Time series feature point distribution

圖7 沉降漏斗研究區Fig.7 Subsidence funnel study area

圖8 各特征點累積沉降Fig.8 Accumulated settlement of characteristic each points

3.2.4 工業區沉降分析

王慶坨鎮和勝芳鎮這一帶區域沉降較為嚴重，出現明顯的沉降漏斗，并且沉降區域有向西延伸的趨勢，針對兩地沉降的情況，單獨分析兩地沉降形成的主要原因。圖9(a)為王慶坨鎮沉降區域形變速率情況，圖9(b)王慶坨鎮沉降區地表利用情況，該地區工業化發展較快，素有“自行車之鄉”之稱的王慶坨鎮是中國北方最大的自行車生產基地，有260多家自行車廠家坐落于此，在工業化發展的同時該地區的城鎮化建設達到了新水平，該地區的工業和生活用水均來自地下水，地下水的過度開采以及地面負荷的增加等因素是造成該地區沉降的主要原因。由圖9可以發現，沉降最為嚴重的區域主要集中在工廠園區，年平均沉降速率最大超過-134 mm/a，最大累積沉降達到了324.6 mm，工業的快速發展帶來了地下水的過度開采以及較大的地表空間利用率對地表沉降起著至關重要的影響。

圖9 王慶坨鎮沉降形變速率及地表空間利用Fig.9 Settlement deformation rate and utilization of surface space in Wangqingtuo town

近年來，勝芳鎮經濟持續發展，隨著中國科學院中試基地和中國金屬玻璃家具研發中心等國家級科研基地落戶于此，勝芳鎮被稱為“中國金屬玻璃家具產業基地”，勝芳鎮工業發展相對較為迅速，主要有金屬玻璃家具產業、鋼鐵產業、板材產業以及鋼鐵產業。由圖10可以看出，勝芳鎮沉降最為嚴重的區域位于郝家堡村，年平均沉降速率達到-134 mm/a，最大累積沉降量達到230 mm，鋼鐵產業和板材產業聚集于此，工業快速發展的同時帶來了地下水的過度開采和地面荷載的增加，城鎮化建設使地表空間利用率增大，地下水的過度開采和地面荷載的增加等因素是勝芳鎮沉降的主要原因。

圖10 勝芳鎮沉降形變速率地表空間利用Fig.10 Settlement deformation rate and utilization of surface space in Shengfang town

2017—2019年，西青區和勝芳鎮兩個地區地下水過度開采導致地下水位基本處于下降狀態。同時，由于工業化和人口增長的影響，兩地的城鎮化區域逐漸增加，因此天津地區地表沉降日益嚴重，且沉降區域有向西延伸的趨勢，需要指出的是，勝芳鎮和西青區之間的一大片區域受耕地和農作物的影響，沉降的連貫性缺失。

4 結論

研究采用SBAS-InSAR技術對天津地區30幅Sentinel-1A影像數據進行處理，揭示了天津地區2017—2019年的地表沉降時空變化特征。得到如下結論。

2017—2019年研究區域內出現明顯的不均勻沉降，天津地區中心城區地表沉降相對平穩，平均沉降速率-5.9 mm/a。王慶坨鎮和勝芳鎮是天津地區兩個最大沉降漏斗，年平均沉降速率分別達到了-134、-114 mm/a。此外研究區內還發現三個沉降較為明顯的區域，即王慶坨鎮和勝芳鎮周邊區域(-95.1～-65.3 mm/a)、這兩個城鎮之間的大片區域(-65.3～-55.4 mm/a)和靜海縣團泊鎮(-104.9～-75.2 mm/a)；將SBAS-InSAR處理結果與同時期天津市二等水準數據進行對比分析，得到同一沉降趨勢和分布。最大相對誤差為8.1 mm，最小相對誤差為-1.54 mm，中誤差為2.96 mm，兩種方法得到的數據結果具有較好的一致性，驗證了利用SBAS-InSAR技術和Sentinel-1A數據用于城市地表沉降監測的可靠性；地下水的過度開采、工業化的迅速發展以及城鎮化建設是天津地區地表沉降的主要因素，同時也是導致天津市地表沉降范圍和幅度不斷增大的主要原因。