小基線技術在玉溪市的地面沉降研究

張航，左小清，李勇發，熊鵬，王志紅，管慶丹，游洪

(1.昆明理工大學 國土資源工程學院，云南 昆明 650093； 2.貴州工程應用技術學院，貴州 畢節 551700)

1 引 言

地面沉降是自然和人為因素引發的，松散地層壓縮導致的地面高程降低的地質現象[1]。近年來，由于我國經濟迅速發展，城鎮化進程加快，城市地面沉降速度明顯加快，地面沉降給工農業、交通物流、城市建設造成嚴重危害，并使居民生活質量下降。目前，中國在21個省份中超過96個城市發生了不同程度的地面沉降，累計沉降量超過 200 mm的總面積超過7.9萬平方公里，地面沉降問題愈發嚴峻，對城市進行地面沉降監測刻不容緩[2，3]。然而傳統的城市地面監測方法如水準測量、三角高程測量、GPS(global positioning system)測量等，不僅存在工作成本高、觀測周期長、觀測過程中測點缺失等缺點，更不能進行大范圍沉降監測[4，5]。合成孔徑雷達干涉測量(interferometric synthetic aperture radar，InSAR)技術的出現，以其全天候、低成本、范圍廣、高精度等優勢在城市地表形變監測中得到廣泛應用，使得對于城市的大范圍變形監測更加快捷、方便[6]。2002年SBAS-InSAR技術的出現，對D-InSAR技術中的大氣影響、失相干等問題，可以通過多景雷達影像和小基線技術，有效降低時空失相干和大氣因素的影響，以得到更加精確的形變結果[7，8]。

近年來，國內外學者運用SBAS-InSAR技術進行城市地面沉降監測，已取得一系列研究成果。2020年莫瑩等利用小基線集技術對紹興市上虞區進行了監測，獲取了該地區的沉降信息，并利用同時期的水準數據驗證了SBAS-InSAR監測結果的可靠性[9]。2018年謝文斌等利用SBAS-InSAR對昆明市進行了沉降監測，得到相關形變結果，昆明最大沉降速率為 -37.6 mm/a[10]。2020年熊鵬等對昆明市進行SBAS-InSAR雙極化實驗，其沉降區域與謝文斌所得到的區域沉降速率基本吻合[11]。2020年Shi Wei等利用了Sentinel-1A與ALOS兩種數據對西安市進行監測，結果表明，西安市在2015年～2019年由于地下水位的上升，出現顯著的局部回彈變形，最大抬升速率為 22 mm/a[12]。2021年冉培廉等采用SBAS-InSAR技術對西安市進行了沉降監測，結果顯示，西安市的地面沉降趨勢正在減緩，但其總體沉降依舊十分嚴重[13]。2017年JimingGuo等將SBAS-InSAR得到的天津市總體沉降結果與2015年地面沉降年報進行對比分析，結果表明，兩者的沉降趨勢和分布基本一致，并運用46個參考站數據驗證了沉降結果的可靠性[14]。2020年張文哲等研究天津市SBAS時序沉降結果，發現與GNSS點觀測結果相吻合，并依據《地面沉降干涉雷達數據測量技術規程》計算出兩者中誤差為 12 mm，滿足技術規程要求[15]。因此，運用SBAS-InSAR技術進行地面沉降監測，可以得到較好的城市地面沉降效果。

目前，還未有學者利用SBAS-InSAR技術對玉溪市進行地面沉降監測，本文利用SBAS-InSAR技術適宜大范圍地面沉降監測的優勢，對玉溪市進行大面積沉降監測，可以克服傳統監測方法的缺陷。本文基于SBAS-InSAR技術選取覆蓋玉溪市中心區域的Sentinel-1A數據影像進行實驗，以得到玉溪市地面沉降結果。并且結合水文地質因素分析玉溪市地面沉降產生的原因，為玉溪市的地面沉降防治提供相關依據。

2 研究區概況與實驗數據

2.1 研究區概況

玉溪市是云南省第五大城市，地處云南省中部，玉溪市東南靠紅河哈尼族彝族自治州，西南連普洱市，東北接省會昆明市，西北鄰楚雄彝族自治州，總面積 15 285 km2，玉溪市地質構造屬于揚子準地臺，位于其西南邊緣。若按二級構造劃分，玉溪市東部為滇東臺褶帶，是元古代至三疊紀的蓋層沉積，主要巖石為震旦系的砂頁巖、碳酸鹽巖，厚度達到了 800 m～2 600 m。并且由于斷層陷落的影響東部存在撫仙湖、星云湖、杞麓湖三個大面積湖泊，其中撫仙湖是我國最大的深水型淡水湖，平均水深達到 95 m。玉溪市西部為川滇臺背斜，地質年代為中晚元古界晉寧亞構造層，地面主要是昆陽群砂泥質板巖以及碳酸鹽巖，巖層厚度在 10 km左右[17]。本文研究區范圍如圖1所示，主要選取紅塔區、江川區全區，以及通海縣、華寧縣、峨山彝族自治縣部分地區，其總面積 2 700余平方公里。

圖1 研究區范圍示意圖

2.2 實驗數據

本文實驗數據，包括星載SAR數據，精密軌道數據，數字高程模型(digital elevation model，DEM)。歐洲航天局(European Space Agency)于2014年發射了Sentinel-1A衛星，該衛星包含一個用于數據采集的C波段合成孔徑雷達，本文實驗主要使用該衛星提供的SLC數據。Sentinel-1A數據為2018年4月～2020年11月范圍內的69景降軌數據，影像主入射角為39.3°，距離向分辨率為 2.3 m，方位向分辨率為 13.9 m，極化方式采用垂直極化(VV)，詳細信息如表1所示。精密軌道數據為哨兵數據對應軌道數據。DEM數據為美國航空航天局(NASA)的SRTM-1 DEM數據，平面分辨率為 30 m。

實驗數據基本信息 表1

3 SBAS-InSAR數據處理

3.1 SBAS-InSAR原理

小基線干涉測量(synthetic baseline sub-set interferometry，SBAS-InSAR)的基本原理是選取多主影像，將多主影像和超級主影像組合出短時間基線、空間基線的干涉圖，然后用奇異值分解(Singular Value Decomposition，SVD)法進行反演，以得到最終的沉降速率結果[16]，SBAS-InSAR技術路線如圖2所示。大致原理如下：

假定某一研究區的影像有N+1幅，選取其中一幅作為主影像，將其余的影像配準到主影像的成像空間上，根據空間基線和時間基線的閾值組合干涉對，得到M幅差分干涉圖，M數量與影像數量有以下關系(圖2)：

圖2 SBAS-InSAR技術路線圖

(1)

對經過兩個不同時刻tA和tB(tA>tB)獲取的影像生成第j景差分干涉圖，該差分干涉圖中的干涉相位可以表示為：

(2)

上式(1)中的j代表任意一景差分干涉圖，dtA和dtB分別為tA和tB時刻相對于t0的累積形變量。當忽略高程和大氣等因素的影響時，可以簡化為：

(3)

為了簡化公式(3)，假定形變速率在相鄰時間間隔以內是線性變化的，在整個的時間段內都是分段線性的，上式(3)的形變相位值可以寫成：

(4)

將處理后的相位寫成矩陣形式為：

Bv=δφ

(5)

式中：B為一個M×N的矩陣。在理想情況下同一小基線合內包含了所有影像，可以用最小二乘法直接進行求解以得到形變相位。但實際上這種可能性很小，因為在進行差分干涉圖組合時，由于設置不同的時間及空間基線閾值，會出現在單個集合內采樣不夠，此時會使得矩陣B出現秩虧。為了得到最小范數意義上的最小二乘解，奇異值分解(SVD)方法可以解決秩虧問題，以求得最終的形變量信息。

3.2 數據處理

本文使用基于Envi5.3平臺的SARscape 5.2.1軟件進行處理。在數據處理之前，首先對哨兵數據進行預處理，將哨兵數據進行數據導入，之后根據本文研究區范圍，進行裁剪。哨兵數據在進行導入、裁剪等預處理操作后，將進行以下五步處理操作，以得到研究區的沉降數據。詳細處理步驟如下：

(1)生成基線連接圖對輸入的69景數據進行配對，生成的連接圖會自動選擇最優的組合方式進行配對，結果會出現時間、空間基線連接圖。生成的像對會進行干涉工作流處理，用于SBAS反演。本文處理數據時自動選擇2018年11月5日的影像作為超級主影像，所有的像對與超級主影像進行配準。由于本文時間跨度大，故時間基線閾值設置為100天，生成的空間連接如圖3所示。

圖3 基線連接圖

(2)干涉處理軟件生成干涉工作流，主要是對通過配對的干涉對進行干涉處理，將所有的數據對都配準到超級主影像上。在多視視數設置時，通過查詢影像的多視視數，設置為5和1。濾波采用Goldstein方法減少噪聲，Minimum Cost Flow方法進行解纏。這一步完成后，查看去平和干涉后的結果，發現干涉圖總體干涉效果達到了預期。

(3)軌道精煉和重去平這一步主要是為了估算和去除殘余的恒定相位以及解纏后還存在的相位坡道。在衛星不準確或者DEM定位不準確的情況下，用地面控制點(GCP)來修正SAR數據。為達到以上目的，主要通過選擇一個去平干涉效果好的干涉圖選取GCP，利用GCP進行重去平。選擇GCP是在干涉圖中無條紋，且平滑之處。

(4)SBAS反演分兩次進行。第一次反演是SBAS反演的核心步驟，第一次估算形變速率地形，并且進行第二次空間解纏。處理完成后加載結果，符合預期結果。第二次反演是在第一次反演結果的基礎上，通過時間域上進行高通濾波，空間域上進行低通濾波，還利用GCP去除殘余的相位，得到更加準確的相位形變信息。

(5)地理編碼該步驟是將處理得到的相位信息，轉化為形變信息，并且將其轉換到地理坐標系上，得到相關的矢量文件和柵格文件。

4 結果與分析

4.1 沉降結果

本文基于SBAS時序分析方法，通過對Sentinel-1SAR影像進行處理，得到玉溪市2018年4月～2020年11月的沉降結果，如圖4所示，有七處明顯的沉降區域(A～G)。其中紅塔區沉降最為嚴重，出現四個明顯的沉降區域(ABCD)，江川區出現兩個沉降漏斗(EF)，通海縣出現一處沉降漏斗(G)。七處沉降漏斗中A、C、E沉降漏斗呈片狀分布，B、D、F、G沉降漏斗面積小，并且較為分散。本文統計沉降點的個數如表2所示，發現研究區內共 3 208 507個沉降點，其中超過85%的沉降點速率在 -13.58 mm/a～6.04 mm/a之間，表明研究區內總體保持穩定狀態，沉降速率在 -13.59 mm/a以下的占比11.91%，抬升速率超過 6.04 mm/a占比僅2.07%，抬升和沉降的區域只占一小部分。

A沉降區主要以北城街道為中心，形成了一個巨大的沉降漏斗，并且其面積達到 30.1 km2，該區域總體的沉降速率在 -20 mm/a左右，累積沉降量最大為 -139.8 mm，最大沉降速率達到 -59.23 mm/a。該區域不僅包含北城街道的中心，還有其周圍的梅園村、蓮池村、劉張營等村落，不僅有建筑用地，還包括一大部分農業用地。分析其沉降原因，發現該研究區不僅城市擴張迅速，北城西側為蓮花池洪積扇，東北部是洪積裙及沖積臺地，向內擠壓造成了大范圍的沉降。并有一定數量的深井開采地下水，導致北城地區地下水水位下降，引起該地區的沉降漏斗。

B沉降區官村-黃泥田一帶，沉降點分散，累積沉降量最大達 -108.59 mm，最大沉降速率 -49.72 mm/a，沉降面積達 3.92 km2。有兩處明顯的形變區，沉降漏斗出現在官村-小犁鏵營以及橄欖灣-黃泥田之間。不難看出，該區域西側為蓮花池洪積扇，受地質因素影響較大，判斷主要因素為二元結構的以粗礫沉積為主的沖積物的影響。

C沉降區大營街-高倉呈現兩片相連的趨勢，尤其是杯湖村到上小屯之間，沉降現象非常嚴重，累積沉降量為 -117 mm，最大沉降速率達到 -72.83 mm/a。沉降區在中間的沖洪積層構成沖積臺的西側，而臺面向西傾斜，造成表層黏土層發生沉降變化。隨著該地區洗浴休閑為主的地熱開發的蓬勃發展，已經造成了地熱水資源的嚴重超采，地熱水屬于一種深層地下水，導致該地區地下水下降區域廣，間接導致了地面沉降的發生。

D沉降區研和街道整體沉降區域分為兩部分即潘井-和樂村與向家莊-南廠村，造成的沉降區域范圍較小，沉降面積為 3.41 km2，但累積沉降量達到了 -146.04 mm，最大值沉降速率達 -67.27 mm/a。分析其沉降原因，研和街道地表巖層由昆陽群淺變質砂板巖組成，即第四系松散巖其厚度在 10 m以上，易造成沉降現象的發生。

E沉降區為江城鎮中心城區，是一個較大的沉降漏斗，該地區位于地處星云湖北岸，撫仙湖西南，江城鎮鎮中心沉降尤為突出，沉降面積達到 7.08 km2，累積沉降量為 -142.38 mm，最大值沉降速率達 -68.67 mm/a。沉降漏斗主要為鎮中心的城市建筑范圍內，該處居民生活用水為撫仙湖水，地下水因素影響較小，該區域城鎮化進程明顯，分析為該地處地面載荷較大引起的地面沉降。

F沉降區位于星云湖南岸，該沉降區沉降范圍分散，存在兩個主要的沉降漏斗，即大寨河入湖區與煙草小區—魚文化廣場之間區域，總沉降面積達到 6.55 km2，累積沉降量為 -138.62 mm，最大值沉降速率達 -80.20 mm/a。分析其原因，大寨河入湖區初步判斷為該區域常年被星云湖湖水浸泡，河流豐水期沖刷河道沖卓導致岸邊土質松軟，易導致地面沉降的發生。煙草小區—魚文化廣場之間，則是一處大型工地，建設過程中造成該地區的地面發生沉降。

G沉降區位于杞麓湖南岸。兩處明顯的沉降區域，即通海縣杞麓湖一號別墅區、獨房子村以及其周邊區域。該區域沉降面積達到 3.83km2，累積沉降量為 -138.25 mm，最大值沉降速率達 -72.70 mm/a。杞麓湖一號別墅區，北面毗鄰杞麓湖，東面有人工湖，地基受湖水影響，造成該區域的沉降較大。獨房子村以及其周邊地區，距離杞麓湖有一定區域，但由于湖水的浸泡，湖浪和沿岸流的沖刷和搬運作用形成各種侵蝕地形和沉積砂體，在農業生產時主要利用地下水進行灌溉，造成該處沉降的發生。

時序沉降點個數統計 表2

地面沉降是一個長期緩慢的過程，不易被人察覺，為了分析玉溪市地面沉降隨時間變化的規律，提取各研究區沉降中心的沉降點(圖4中1-7點)，形成如圖5的時序累積沉降圖。從圖中可以明顯看出，七處沉降區域沉降隨時間的增加而逐漸加重。1(A)在2019年3月～2019年11月沉降加劇，沉降量達到 -40 mm，在2019年11月后，沉降趨勢減緩。2(B)在2018年9月～2019年3月和2019年9月～2020年3月，整體保持穩定，分別在 -20 mm和 -60 mm上下波動，在2019年4月～2019年8月，急劇沉降 -30 mm。3(C)整體呈現線性分布，只在2020年5月～2020年7月沉降加速。4(D)在2019年2月～2019年7月，其沉降量達到了 -50 mm，其總體沉降量達到了 -132.55 mm。5(E)和7(G)兩者沉降趨勢大概一致，在2018年6月～2019年6月，兩者差距在 20 mm，之后趨于重合，7(G)在2020年9月，急劇沉降，最后達到了 -129.92 mm。6(F)在2018年12月～2019年7月沉降速度加快，之后趨于緩慢，但在2019年7～10月與2020年5月數據出現異常，分析為兩期數據去噪效果差。

圖5 玉溪市典型沉降區累積形變序列

4.2 形變典型區域分析

為了進一步揭示沉降嚴重地區的垂直空間特征，本文選取A、C、E三個典型沉降漏斗區域，提取剖面線a-a′、b-b′、c-c′，位置如圖6所示，得到沉降剖面線如圖7所示。

圖6 玉溪市三大沉降區域示意圖

圖7 沉降剖面線圖

A區中北城橫斷面以大米飯館為起點至春和居委會北部為終點，跨度約為 6 km。該區域出現了兩個典型的“V”字形沉降漏斗(1)與(2)，分別位于陳大場村南部和高橋東南部，最大沉降速率分別為 -44.83 mm/a、-40.73 mm/a，該區域總體平均沉降速率在 -25 mm/a左右，并且兩個沉降漏斗有進一步發展的趨勢。C區中為大營街、高倉一帶橫斷面從蔡官營開始至棠梨樹東部結束，該區域有兩個“U”字形沉降漏斗即(3)和(4)，(3)沉降漏斗的寬度約為 0.6 km，沉降漏斗沉降速率在 -35 mm/a以下，(4)沉降漏斗的寬度約為 1.4 km，該沉降區域沉降速率在 -30 mm/a以下。從兩個沉降漏斗的發展來看，兩個沉降漏斗的下沉趨勢正在減緩。E區江城鎮橫截面起點自伊旗村東北與張官營交接處，終點在云巖村東北部，沉降區的跨度約為 4.5 km，該區域整體沉降速率在 -30 mm/a左右，在約 1.5 km處有一個明顯的沉降漏斗(5)，江城鎮財政所最大沉降速率達到 -55.99 mm/a，該區域存在兩處沉降脊，沉降脊兩側沉降漏斗有進一步發育的趨勢。

4.3 形變原因分析

現代社會每年用水量巨大，大致用途可分為農業、工業、服務業用水。而玉溪市內有珠江、紅河等水系以及星云湖、杞麓湖、撫仙湖等湖泊，主要居民點的服務業用水以地表水為主，少部分以地下水供水。玉溪市的工農業用地占土地總面積的15%左右，工農業用水量超過總用水量的80%，地表用水已不足以供應其用水，絕大部分抽取地下水用以生產，造成地下水水位一直處于下降狀態。2000年以來，而伴隨著休閑服務行業的興起，玉溪市內出現了多家以地熱水為主的洗浴休閑中心，僅紅塔區內處于開采中的地熱井估計在10眼以上，年采熱水超過122.6×104m3，將紅塔區的沉降區域與張苗紅[19]研究得到的地下水沉降區域進行對比發現，兩者有較高的吻合度，表明地下水水位的下降是導致地面沉降的一個重要因素。

5 結 論

本文基于2018年7月～2020年11月Sentinel-1A的69景SAR數據，提取出玉溪市的沉降速率圖。結果表明，有7處明顯的形變區，最大累計沉降量為 -175.27 mm，最大年沉降速率達 -85.25 mm/a。北城沉降區沉降面積最大，達到了 30.1 km2。分析7處典型區域中心點的累積形變序列，發現在2019年3月～2019年9月沉降速度加快，之后沉降速度減緩。并對3個典型區域進行垂直空間特征分析，大部分漏斗的沉降速率在 -30 mm/a～-40 mm/a范圍內，江城鎮沉降漏斗最大沉降速率達到 -55.99 mm/a。通過進一步分析玉溪市沉降原因，發現玉溪市的經濟快速發展導致的地下水位下降是主要原因。