李蓉蓉，楊維芳*，李得宴

(1. 蘭州交通大學(xué) 測繪與地理信息學(xué)院，蘭州 730070；2. 地理國情監(jiān)測技術(shù)應(yīng)用 國家地方聯(lián)合工程研究中心，蘭州 730070；3. 甘肅省地理國情監(jiān)測工程實驗室，蘭州 730070)

地面沉降是目前各大城市的一個主要工程地質(zhì)問題，它給城市的經(jīng)濟(jì)和公共安全帶來極大的危害.蘭州市位于青藏高原東北側(cè)的黃河河谷盆地內(nèi)，地質(zhì)災(zāi)害風(fēng)險指數(shù)和風(fēng)險等級在黃河上游地區(qū)處于首位[1].目前該區(qū)域地面沉降具體情況未知，監(jiān)測措施極少，采用一定的監(jiān)測手段對其沉降狀況進(jìn)行定量研究并分析其時空演化特征及其主要驅(qū)動力因素對該市的健康發(fā)展和地質(zhì)災(zāi)害的早期識別和預(yù)防具有重要意義.

傳統(tǒng)的水準(zhǔn)和GNSS監(jiān)測手段，監(jiān)測精度高，但只能獲得布設(shè)點上的沉降狀況，并不能獲取連續(xù)的沉降信息.相較于傳統(tǒng)的沉降監(jiān)測手段，合成孔徑雷達(dá)干涉測量技術(shù)(interferometry synthetic aperture radar,InSAR)可以進(jìn)行高精度、高空間分辨率、大范圍的沉降監(jiān)測.目前常用于地面形變監(jiān)測的InSAR技術(shù)主要有合成孔徑雷達(dá)差分干涉測量(differential InSAR,D-InSAR)技術(shù)、永久散射體干涉測量(persistent scatterer InSAR,PS-InSAR)技術(shù)以及短基線集干涉測量(small baseline subsets InSAR,SBAS-InSAR)技術(shù)等[2-5].但是D-InSAR技術(shù)易受大氣擾動、時空基線去相干的影響，PS-InSAR技術(shù)對研究區(qū)及影像數(shù)量有較高要求，而SBAS-InSAR技術(shù)降低了對影像數(shù)量的要求并且在山區(qū)也能監(jiān)測.目前SBAS-InSAR已被廣泛用于地面沉降監(jiān)測并且精度和可靠性也得到了驗證[6-9].標(biāo)準(zhǔn)差橢圓(standard deviation ellipse,SDE)法是一種基于空間統(tǒng)計來揭示地理要素空間變化特征的有效方法.該方法通過標(biāo)準(zhǔn)差橢圓各參數(shù)的變化來定量解釋要素空間分布的時空變化特征，其直觀性與有效性已得到廣泛應(yīng)用[10-12].因此可利用標(biāo)準(zhǔn)差橢圓法來定量分析地面沉降的時空變化特征.

蘭州市城區(qū)地面沉降頻發(fā)，地質(zhì)條件特，地形變化復(fù)雜[13]，很難找到足夠的PS點，PS-InSAR技術(shù)并不適合該地沉降監(jiān)測，故本文利用SBAS-InSAR技術(shù)監(jiān)測該市地面沉降狀況，利用SDE法、時序分析法對地表沉降的時空間變化特征進(jìn)行分析并結(jié)合各種資料深入探討引起地面沉降的原因，為蘭州市主城區(qū)規(guī)劃發(fā)展和地質(zhì)災(zāi)害早期識別、防治提供參考.

1 研究區(qū)與實驗數(shù)據(jù)

1.1 研究區(qū)

蘭州市地處中國西北地區(qū)，地勢東北較低，西部、南部較高，黃河自西向東穿城而過，城區(qū)東西狹長，南北由兩山夾峙.蘭州市地質(zhì)條件復(fù)雜，是國內(nèi)少有的受地質(zhì)災(zāi)害困擾的省會城市，其土質(zhì)類型為黃土，黃土結(jié)構(gòu)疏松，在水力、風(fēng)力、重力和人為因素影響下發(fā)生地質(zhì)災(zāi)害的可能性極大[14].本文以蘭州市主城區(qū)為研究區(qū)，研究區(qū)地理位置如圖1所示，經(jīng)度范圍為103°34′-103°58′，緯度范圍為35°59′-37°10′，主要包括城關(guān)區(qū)、安寧區(qū)、七里河區(qū)、西固區(qū)、榆中縣以及皋蘭縣的部分區(qū)域.

1.2 實驗數(shù)據(jù)

本文選取覆蓋研究區(qū)的29景Sentinel-1A，升軌影像數(shù)據(jù)進(jìn)行實驗，數(shù)據(jù)類型為干涉寬幅模式的單視復(fù)數(shù)數(shù)據(jù).時間跨度為2017-07-13至2020-04-28，極化方式為VV同向極化.此外實驗還用到30 m分辨率的SRTM1高程數(shù)據(jù)以及各影像對應(yīng)的精密定軌星歷數(shù)據(jù).

2 研究方法與數(shù)據(jù)處理

2.1 SBAS-InSAR方法

SBAS-InSAR是一種基于分布式目標(biāo)算法的時序InSAR技術(shù)，能夠在一定程度上克服時間失相關(guān)、空間失相關(guān)的不利影響.該方法首先基于多項式形變模型解算低頻形變，同時對高程誤差進(jìn)行建模和解算，然后從原始差分干涉相位中扣除低頻形變和高程誤差后，再重新對殘差相位進(jìn)行相位解纏，再將低頻形變分量加回，此時再利用上述解算過程求解每個時間段內(nèi)的相位變化速率并恢復(fù)相位時間序列.最后，從相位時間序列中扣除之前求得的低頻形變，對殘差相位時間序列進(jìn)行時空濾波得到大氣延遲相位時間序列，從原始相位時間序列中扣除大氣延遲相位即可得到形變相位序列，通過相位到形變的轉(zhuǎn)換可得到形變時間序列.

圖1 研究區(qū)地理位置Fig.1 Geographical location of the study area

2.2 SDE法

SDE是一種基于空間統(tǒng)計來揭示地理要素空間分布特征的有效方法.本文采用標(biāo)準(zhǔn)差橢圓的圓心、面積、長短軸之比以及方位角的變化來衡量該區(qū)域地面沉降中心、范圍、發(fā)展主方向及發(fā)展主方向與正北方向的偏離角度的變化.由于SBAS-InSAR技術(shù)獲得的地面沉降量是WGS-84坐標(biāo)系下的，故本文以SBAS監(jiān)測結(jié)果矢量化得到的m個點的地面沉降量為權(quán)重建立WGS-84坐標(biāo)系下的加權(quán)標(biāo)準(zhǔn)差橢圓.

本文在計算加權(quán)標(biāo)準(zhǔn)差橢圓時主要計算了圓心、方位角、以及橢圓長短軸的標(biāo)準(zhǔn)差.

加權(quán)標(biāo)準(zhǔn)差橢圓的圓心(SDEx,SDEy)計算公式為

(1)

(2)

其中：

(3)

加權(quán)標(biāo)準(zhǔn)差橢圓的方位角θ的tanθ計算公式為

(4)

加權(quán)標(biāo)準(zhǔn)差橢圓在長軸和短軸方向上的標(biāo)準(zhǔn)差(δx,δy)計算公式為

(5)

(6)

其中：m表示SBAS柵格矢量化得到的點的個數(shù)；wk表示第k個點的地面沉降量；xk和yk表示第k個點的坐標(biāo)；(xc,yc)表示m個點的坐標(biāo)的加權(quán)平均中心.

2.3 數(shù)據(jù)處理

SBAS-InSAR數(shù)據(jù)處理采用ENVI中的SARscape進(jìn)行，首先對數(shù)據(jù)進(jìn)行預(yù)處理，設(shè)置最大時空基線閾值分別為365 d、45%，生成影像連接圖.選擇2018-02-11影像作為超級主影像，進(jìn)行主輔影像配準(zhǔn).相位解纏方法選擇最小費用流法(minimum cost flow,MCF)，解纏相干系數(shù)閾值為0.35，濾波方法選擇Goldstein濾波，進(jìn)行干涉工作流處理.選擇35個高相干的控制點估算并去除殘余的恒定相位和解纏后仍存在的相位坡道，基于線性模型計算出所有像對的形變和高程，并采用MCF進(jìn)行再次解纏.最后設(shè)置相干系數(shù)閾值為0.4，空間濾波窗口為1 200 m×1 200 m，時間濾波窗口為365 d，進(jìn)行定制的大氣濾波，估算和去除大氣相位，得到時間序列上的最終形變結(jié)果.

SDE及其要素的計算，本文采用ArcGIS空間統(tǒng)計工具中的標(biāo)準(zhǔn)差橢圓工具來計算.分別輸入SBAS各時段的累計沉降量作為權(quán)重，設(shè)置橢圓大小為包含點數(shù)量為68%，運行即可得到標(biāo)準(zhǔn)差橢圓及其參數(shù).

3 沉降結(jié)果及其時空變化特征分析

3.1 SBAS-InSAR監(jiān)測結(jié)果

本文采用SBAS-InSAR技術(shù)對蘭州市主城區(qū)2017-07-13至2020-04-28間的地面沉降進(jìn)行了監(jiān)測，結(jié)果如圖2～3所示.圖中負(fù)值表示地面沉降，正值表示地面抬升.由圖2可以看出該技術(shù)得到的年均沉降速率范圍為[-35.5,22.3] mm/a，研究區(qū)中沉降速率較大的區(qū)域主要分布在蘭州市文創(chuàng)城(A)、城關(guān)區(qū)青白石片區(qū)(B)、城關(guān)區(qū)恒大附近的物流園(C)、蘭州資源環(huán)境職業(yè)技術(shù)學(xué)院(D)、伏龍坪附近(E)、西固區(qū)孔家營和杏胡臺村附近(F).年均沉降速率最大的區(qū)域位于城關(guān)區(qū)青白石片區(qū)附近，它的最大沉降速率可達(dá)35.5 mm/a.從圖3可以看出在監(jiān)測期間研究區(qū)的累計沉降量范圍為[-117.6,42.7] mm，城關(guān)區(qū)青白石片區(qū)的累計沉降量最大.監(jiān)測期間在研究區(qū)東南部的榆中縣及城關(guān)區(qū)部分區(qū)域出現(xiàn)了大范圍的沉降并且累計沉降量較大，隴海鐵路及蘭新高鐵等鐵路沿線出現(xiàn)了明顯的帶狀沉降區(qū)域，南北兩山山麓地帶的一些村莊及土地開發(fā)區(qū)內(nèi)也發(fā)生了較為明顯的沉降.地質(zhì)災(zāi)害防治部門應(yīng)加大對這些沉降嚴(yán)重區(qū)域的監(jiān)管.

圖2 蘭州市2017-2020年年均沉降速率Fig.2 The average annual subsidence rate of Lanzhou from 2017 to 2020

圖3 蘭州市2017-2020年累計沉降量Fig.3 The cumulative settlement of Lanzhou from 2017 to 2020

SBAS-InSAR處理結(jié)果柵格矢量化得到1 215 462個點，將這些點的年均沉降速率和累計沉降量按不同的量值區(qū)間進(jìn)行分類并統(tǒng)計年均沉降速率和累沉降量在各區(qū)間內(nèi)點數(shù)如圖4所示.

圖4 年均沉降速率、累計沉降量統(tǒng)計Fig.4 Statistics of annual average settlement rate and cumulative settlement

由圖4(a)可知，對于年均沉降速率，分布在[-5,5] mm/a區(qū)間內(nèi)點數(shù)最多，分布在[10,25] mm/a區(qū)間點數(shù)最少.統(tǒng)計約有50%的點的年均沉降速率處在0 mm/a之下.由圖4(b)可知對于累計沉降量，分布在[-10,0] mm區(qū)間內(nèi)點數(shù)最多，分布在[0,10] mm區(qū)間的百分比僅次之，分布在[-120,80] mm區(qū)間內(nèi)的點數(shù)最低，約有67%的點的累計形變量處在0 mm之下.經(jīng)以上分析發(fā)現(xiàn)監(jiān)測期間，蘭州市城區(qū)整體上發(fā)生了沉降，但沉降速率及累計沉降量均較小.

3.2 地面形變時空變化特征分析

本文采用SDE法和時序分析法來定量分析監(jiān)測期間蘭州市地面沉降時空變化特征.

首先采用SDE法定量分析蘭州市沉降的空間變化特征.本文以SBAS-InSAR技術(shù)柵格轉(zhuǎn)矢量生成的數(shù)據(jù)中2017-08-18、2018-07-08、2019-06-09、2020-04-28的累計沉降量分別作為權(quán)重，計算得到各時間的加權(quán)標(biāo)準(zhǔn)差橢圓，對各時間上標(biāo)準(zhǔn)差橢圓的參數(shù)進(jìn)行統(tǒng)計如表1所列，繪制各橢圓參數(shù)隨時間變化如圖5所示.由圖5(a)可以看出，整體上標(biāo)準(zhǔn)差橢圓的中心點向東北方向移動即地面沉降中心向東北方向移動；由圖5(b)可以看出標(biāo)準(zhǔn)差橢圓面積先減小后增大又減少，故監(jiān)測期間蘭州市地面沉降范圍呈現(xiàn)出減少增加減少的趨勢；由圖5(c)中可以看出短軸/長軸的比值均小于0.35，地面沉降具有明顯的方向性.從2017-08-18至2020-04-28標(biāo)準(zhǔn)差橢圓短軸明顯縮短，沉降在短軸方向上即西南-東北方向發(fā)展相對減緩；由5(d)可知，從2017-08-18至2020-04-28標(biāo)準(zhǔn)差橢圓的方位角由94.178 °增加至98.950 °，在此期間地面沉降的標(biāo)準(zhǔn)差橢圓在空間演化上表現(xiàn)為順時針旋轉(zhuǎn).綜上所述，從整體上看監(jiān)測期間蘭州市地面沉降的中心逐漸向東北方向移動，形變范圍逐漸減小，在西北-東南方向上具有明顯的方向性，西南-東北方向上沉降相對減緩，地面沉降不斷朝順時針方向發(fā)展，并且沉降重心向黃河北岸移動.

表1 標(biāo)準(zhǔn)差橢圓參數(shù)

圖5 標(biāo)準(zhǔn)差橢圓各參數(shù)變化Fig.5 Changes of various parameters of the standard deviation ellipse

其次采用時序分析法分析累計沉降量隨時間變化特征.以B、C、D、F區(qū)域為例，對各區(qū)域內(nèi)SBAS點(SBAS-InSAR柵格轉(zhuǎn)矢量生成的點)在各個時相上的累計沉降量求平均值，繪制各區(qū)域內(nèi)平均累計沉降量隨時間變化曲線，分析各區(qū)域累計沉降量隨時間變化特征，結(jié)果如圖6所示.從圖6可以看出D區(qū)域的平均累計沉降量最大，約為32.1 mm，C區(qū)域僅次之約31.2 mm，B區(qū)域的平均累計沉降量最小，約為16.5 mm.C、D區(qū)域的平均累計沉降量達(dá)到30 mm以上，地質(zhì)災(zāi)害防治部門應(yīng)加大對這兩個區(qū)域的監(jiān)測.從整個監(jiān)測期間來看，各區(qū)域的平均累計沉降量隨時間波動性較大并且呈現(xiàn)非線性下降趨勢，在2019-07-15之后沉降速率均明顯增大，沉降速率和累計沉降量均有進(jìn)一步加大的趨勢.

圖6 各區(qū)域平均累計沉降量隨時間變化折線Fig.6 Broken line of the average cumulative settlement of each region over time

4 沉降原因分析

監(jiān)測期間在蘭州市南北兩山坡腳地帶的一些村莊、鐵路沿線、土地開發(fā)區(qū)、物流園附近等地均出現(xiàn)明顯的沉降.這些區(qū)域發(fā)生沉降主要與蘭州市的土體類型、建筑施工建設(shè)運營及降水量有關(guān).

4.1 土體類型對地面沉降的影響

蘭州市城區(qū)土體類型主要包括碎礫土、砂性土、粘性土、黃土、黃土狀土以及少量淤泥質(zhì)土和人工填土.砂性土分布零星，黃土、黃土狀土分布最為廣泛.IV級階地及IV級以上丘陵均為黃土分布區(qū)，河谷地帶多為黃土狀土.將蘭州市城區(qū)地表沉降空間分布特征與南北兩山、黃河高階地地區(qū)地質(zhì)類型分布特征對比，發(fā)現(xiàn)兩地表下沉最為顯著的幾個區(qū)域均分布在黃土、黃土狀土上.I、II、III、IV級階地的黃土狀土及黃土都具有濕陷性[15].這些區(qū)域土質(zhì)結(jié)構(gòu)疏松、孔隙率較高，當(dāng)土層處于欠壓密狀態(tài)時受水浸濕，土層結(jié)構(gòu)會迅速破壞，強度降低，在自重應(yīng)力和附加應(yīng)力共同作用下發(fā)生顯著附加下沉[16].

4.2 建筑施工等工程效應(yīng)對地面沉降的影響

蘭州市為了緩解用地緊張和人口激增所面臨的居住問題，修建了大量的高層建筑.建筑在開挖大基坑時，會造成支護(hù)結(jié)構(gòu)失穩(wěn)，導(dǎo)致基坑周邊地區(qū)地面沉降[17].深基坑和地下建筑物開挖時，為防止邊坡基坑失穩(wěn)和水下作業(yè)，會進(jìn)行基坑降排水，降水方案不當(dāng)會引起周圍地面沉降.施工完成時，在建筑物荷載的作用下，土體產(chǎn)生附加壓力，土體的超靜水壓力迫使土中水產(chǎn)生流動的滲流，形成水頭差，土地孔隙比改變，飽和土體在外荷載作用下發(fā)生明顯的沉降.從圖3可以明顯看出監(jiān)測期間在隴海鐵路、包蘭鐵路、蘭新鐵路沿線出現(xiàn)了明顯的沉降.這主要與列車的運營和建設(shè)有關(guān).列車在運營期間自身的重力以及長時間往復(fù)施加的循環(huán)荷載會對土體的動力穩(wěn)定性構(gòu)成影響，在循環(huán)荷載振動下土性以及孔隙水壓力會發(fā)生變化[18]，進(jìn)而導(dǎo)致地面沉降.在施工建設(shè)期間對地面沉降的影響主要是由于進(jìn)行地下工程施工會引起地層原始應(yīng)力狀態(tài)的改變、土體的固結(jié)、土體蠕變效應(yīng)以及地層損失等.

4.3 降水量對地面沉降的影響

本文以A、D區(qū)域為例，分別獲取A、D區(qū)域的平均時序累計沉降量及月平均降水量，分析A、D區(qū)域的平均時序累計沉降量與月平均降水量之間的關(guān)系，結(jié)果如圖7所示.由圖7可知,在A、D區(qū)域內(nèi)每年1、2、3、12月份降水量相對較少，而對應(yīng)時段的沉降速率在增大，沉降量在增加.這主要是由于在這幾個月大氣降水較少，地面較為干旱，地下水的抽取力度加大，含水層的水位降低，孔隙應(yīng)力減小，隔水層的有效應(yīng)力增大，引起了隔水層粒間骨架的壓縮，從而產(chǎn)生地面沉降[19].A、D區(qū)域在每年7-9月份降水量是全年最大的時期，對應(yīng)時段內(nèi)的沉降速率較小，沉降量變化量較少.這主要是由于在這幾個月里，降水量較大，降水通過地下深水井或者土壤滲透到地下水位，使對應(yīng)的地下水位呈現(xiàn)上升狀態(tài)，孔隙水壓力減弱，有效的緩解了地面沉降速率，從而導(dǎo)致沉降量減少.通過以上分析可知，降水量的多少影響著沉降速率及累計沉降量，當(dāng)降水量增多時地表沉降有所減緩，當(dāng)降水量減少時地表沉降加劇.

圖7 平均累計沉降量與月平均降水量關(guān)系Fig.7 Relationship between average cumulative subsidence and monthly average precipitation

5 結(jié)論

本文采用SBAS-InSAR技術(shù)對蘭州市2017-07-13至2020-04-28間的地面沉降狀況進(jìn)行監(jiān)測并且詳細(xì)分析了沉降的時空變化特征和影響因素，得到以下結(jié)論：

1) 監(jiān)測期間蘭州市城區(qū)整體上發(fā)生了沉降，但沉降速率及累計沉降量均較小.榆中縣沉降區(qū)域最多，年均沉降速率及累計形變量最大的區(qū)域分布在城關(guān)區(qū)青白石片區(qū).蘭州資源環(huán)境職業(yè)技術(shù)學(xué)院附近的平均累計沉降量最大，恒大山水城附近的物流園僅次之.地質(zhì)災(zāi)害防治部門應(yīng)加大對恒大附近的物流園以及蘭州資源環(huán)境職業(yè)技術(shù)學(xué)院附近沉降嚴(yán)重區(qū)域的監(jiān)管.

2) 在監(jiān)測期間，蘭州市地面沉降中心先向東北方向移動再向西北方向移動，整體上向東北方向移動；沉降范圍呈現(xiàn)出先增加后減少的趨勢；地面沉降具有明顯的方向性，在東南-西北方向上沉降明顯，在西南-東北方向發(fā)展相對減緩并且沉降不斷朝順時針方向發(fā)展.累計沉降量隨時間均呈現(xiàn)非線性下降趨勢，并且沉降速率和累計形變量有進(jìn)一步加大的趨勢.

3) 監(jiān)測期間，蘭州市在南北兩山坡腳地帶、鐵路沿線、以及部分物流園附近出現(xiàn)了明顯的沉降，這主要與蘭州市土體類型、建筑和施工等工程效應(yīng)以及降水量有關(guān)，并且它們對地面沉降的影響具有明顯的復(fù)合效應(yīng).

毋庸置疑，城市需要開發(fā)建設(shè)，但是在開發(fā)建設(shè)的過程中，一定要科學(xué)規(guī)劃，合理的規(guī)避這些沉降嚴(yán)重的區(qū)域并且地質(zhì)災(zāi)害防治部門也要加強對沉降嚴(yán)重區(qū)域的監(jiān)管，及時了解沉降狀況，這對城市的健康發(fā)展和地質(zhì)災(zāi)害的預(yù)防具有重要意義.