戴宜韋，吳希文，閔新穎，王 華，朱煥廉，彭林才

（廣東工業(yè)大學(xué) 土木與交通學(xué)院，廣東 廣州 510006）

地面沉降是全球性問題，由于人類活動的影響，加劇了這種地質(zhì)災(zāi)害發(fā)生的程度與范圍. 現(xiàn)如今，全球性氣候變暖明顯，導(dǎo)致海平面上升，地面沉降對沿海地區(qū)的影響尤為嚴重. 據(jù)國土資源部地質(zhì)災(zāi)害應(yīng)急辦不完全統(tǒng)計，僅僅2017年上半年我國由于地質(zhì)災(zāi)害造成的直接經(jīng)濟損失就達到2.6億元.

近年來，隨著合成孔徑雷達技術(shù)(InSAR，Synthetic Aperture Radar Interferometry)的發(fā)展，InSAR是一種利用微波波段特性進行相干成像系統(tǒng)技術(shù)，在監(jiān)測地面沉降方面已經(jīng)成為一種可靠的新手段. 1999年，F(xiàn)erretti等提出永久散射體干涉測量PSInSAR(Permanent ScattererInSAR)技術(shù)，利用多景SAR圖像從中找出散射特性穩(wěn)定的高相干點作研究對象，并成功獲取長時間序列的大范圍地面形變緩慢的演化過程[1]. 在先前研究中，2014年王華等使用美國噴氣推進實驗室(JPL)和加州理工學(xué)院聯(lián)合開發(fā)的ROI-PAC軟件對覆蓋廣佛地區(qū)的Envisat衛(wèi)星影像進行D-InSAR干涉圖生成，并利用π-rate軟件對干涉圖進行后續(xù)分析，獲得了廣佛地區(qū)的地面沉降速率與時間序列，驗證了利用InSAR技術(shù)在廣佛地區(qū)的沉降監(jiān)測可以達到毫米級精度，廣佛地區(qū)趨于平穩(wěn)狀態(tài)，局部地區(qū)發(fā)現(xiàn)有沉降[2]. 2018年Ng，A.H.M等利用時間序列InSAR原理對廣佛地區(qū)2011年至2017年COSMO-SkyMed衛(wèi)星數(shù)據(jù)進一步進行地面沉降分析，整體區(qū)域的沉降速率范圍在-35～10 mm/a之間[3].由于原先的雷達數(shù)據(jù)很大部分都是需要購買的或需要向衛(wèi)星所屬機構(gòu)申請獲得，有一定的局限性. 隨著2014年Sentinel-1A衛(wèi)星的發(fā)射成功，2015年4月5日開始穩(wěn)定運行，數(shù)據(jù)免費向公眾開放，利用Sentinel-1A數(shù)據(jù)監(jiān)測地面沉降可以進一步減少所需費用. 在過去的研究中并沒有人利用過Sentinel-1A來研究珠江三角洲廣佛及東莞地面變形情況，本文研究證明了Sentinel-1A可以應(yīng)用于珠江三角洲地面變形監(jiān)測，為以后Sentinel-1A應(yīng)用于監(jiān)測氣候濕潤的三角洲地區(qū)的地面變形提供了參考依據(jù).

本文利用時間序列InSAR技術(shù)對珠江三角洲主要城市的主城區(qū)進行沉降觀測，獲取了廣州、佛山、東莞在內(nèi)的沉降演變特征，并分析各個城市的地面沉降速率圖，預(yù)測出潛在的危險沉降點.

1 時序InSAR原理

Persistent Scatterer Interferometry(PSI)技術(shù)的基本原理是利用多景同一地區(qū)的SAR影像，通過統(tǒng)計分析時間序列上的幅度與相位信息的穩(wěn)定性，探測不受時間、空間基線去相關(guān)影響的穩(wěn)定點目標(biāo)[4-5]. 對探測出的穩(wěn)定目標(biāo)點進行時間序列分析，從而反演出研究區(qū)地物的形變量. 其主要處理在于干涉圖生成、利用DEM進行差分干涉圖計算、候選點(PSC)選取以及對候選點進一步進行篩選得到最終的PS點.相位解纏所使用的函數(shù)模型如式(1)所示.

其中， φtopo表示由參考DEM誤差引起的地形相位；φdefo表示點的形變相位；φatmo表示大氣延遲相位；φnoise表示不相干的噪聲相位. 在觀測相位中含有φatmo大氣延遲相位，為了去除減弱大氣相位對最終結(jié)果的影響，通常通過相鄰PS點的大氣相位空間相關(guān)性，采用基于相鄰PS點相位差分模型提取地表沉降信息[6]. 一般對PSI做的準備工作包括:基于振幅離差指數(shù)方法對PS候選點(PSC)選取，基于整體相干性最大法估計出線性形變速度增量及高程誤差以及對大氣相位的濾波去除[4]. 對于大氣相位的消除，是基于大氣相位在空間區(qū)域相關(guān)度很高，呈現(xiàn)出低頻狀態(tài)，在時間域上大氣相位相關(guān)度不高，呈現(xiàn)出高頻狀態(tài)的特征[2，4，7]. 從而采用在時間上高通濾波，在空間上采用低通濾波來分離出非線性沉降形變與大氣延遲相位;即整體相位減去模型相位得殘余相位，殘余相位中包含大氣相位(APS)、非線性形變與其余隨機誤差(noise)[2，4，7]. 隨機誤差對最后的結(jié)果影響較少，可以忽略不計. 最終得到研究區(qū)的形變相位，進行時序分析.

2 研究區(qū)概況

珠江三角洲的植被覆蓋率高，空氣濕熱，年均溫度21～23 ℃，年降水量大約1 700 mm，以南亞熱帶季風(fēng)氣候為主. 多雨季節(jié)與高溫季節(jié)同步，土壤肥沃，河道縱橫. 珠江三角洲陸上部分由西江斷裂、瘦狗嶺斷裂、南崗—太平斷裂所圍攏的斷陷盆地組成. 根據(jù)廣東地質(zhì)調(diào)查局的報告，有三種地質(zhì)構(gòu)造，即該地區(qū)有軟、疏松的沉積物、埋藏的巖溶和其他類型的基巖. 前兩個地質(zhì)單元均有沉降，對周圍建筑物的安全構(gòu)成威脅，基巖由分布在增城丘陵地區(qū)和廣花盆地北部的不同類型的花崗巖組成[3，8]. 三角洲及濱海平原軟土區(qū)則有軟土地基變形、砂土液化、河港淤積，脹縮土等分布區(qū)，常見的地質(zhì)災(zāi)害有地裂縫，部分地區(qū)出現(xiàn)地面沉降[9]. 本文研究的主要城市主城區(qū)如圖1所示

圖1 研究區(qū)域(廣州、佛山和東莞)Fig.1 Research area (Guangzhou, Foshan and Dongguan)

3 Sentinel-1A雷達數(shù)據(jù)及處理方法

本文數(shù)據(jù)處理的范圍覆蓋廣州、佛山和東莞等主城區(qū). 本文處理數(shù)據(jù)采用2015年12月～2018年5月的Sentinel-1A衛(wèi)星數(shù)據(jù)，總共66景IW模式SAR影像.具體衛(wèi)星與數(shù)據(jù)參數(shù)見表1.

表1 Sentinel-1A衛(wèi)星數(shù)據(jù)參數(shù)Tab.1 Sentinel-1A satellite data parameters

本次研究共處理了66景Sentinel-1A衛(wèi)星影像. 衛(wèi)星軌道數(shù)據(jù)采用歐空局(ESA)發(fā)布的精密軌道數(shù)據(jù)，本文采用ESA發(fā)布的精密軌道數(shù)據(jù)，可以認為軌道殘余誤差對最終結(jié)果的影響很小，進一步忽略軌道誤差造成的影響. 采用美國JPL發(fā)布的SRTM的1’DEM數(shù)據(jù)去除地形相位. 本文采用的66景影像及生成65幅干涉圖對應(yīng)的時間基線與垂直基線如圖2所示.

圖2 Sentinel-1A空間基線與時間基線分布圖Fig.2 Sentinel-1A spatial baseline and time baseline distribution map

首先對于時間序列InSAR，選取覆蓋珠江三角洲的66幅單視復(fù)數(shù)圖像(SLC)進行輻射矯正，根據(jù)時空基線，選取獲取日期2017年3月12日作為主影像，其余影像作為輔影像，生成65幅干涉圖. 利用通過設(shè)定幅度離差閾值方法進行PS候選點選取，選取幅度離差指數(shù)低于0.25的PS點構(gòu)建PS網(wǎng)絡(luò). 通過導(dǎo)入外部DEM生成差分干涉圖，進行模型的反演，估算出殘余地形與形變速率. 基于振幅離差指數(shù)[1，6]自動選出一個或者多個像素(參考點)，參考點設(shè)置為25 km2一個參考點，計算出相位偏移[6]. 在干涉圖生成后，就從所有的干涉圖中去除相位偏移. 第一次反演后用線性模型從所有差分干涉圖中估算的形變速率和殘余高程信息，殘差是大氣信息(噪聲和非線性目標(biāo)). 對相干性高的像元進行計算，進行第二次反演，得到估算大氣延遲相位，去除大氣延遲相位得到最終的形變速率[6]；最后對沉降速率、高程殘差以及時間沉降量進行地理編碼，投影到當(dāng)前所用的WGS84地理坐標(biāo)系.

4 研究區(qū)沉降監(jiān)測結(jié)果

由于本文所用衛(wèi)星數(shù)據(jù)從2015年12月至2018年5月底，所得出的研究成果通過跟2015年之前的研究及本文研究期間他人利用其他衛(wèi)星平臺所發(fā)表的研究成果進行比較，從而驗證本文研究成果的準確性.

4.1 廣佛地區(qū)

本文最終計算珠江三角洲廣佛主城區(qū)地面沉降速率精度誤差及精度誤差統(tǒng)計如圖3所示.

圖3 Sentinel-1A廣佛雷達衛(wèi)星影像沉降速率精度圖及精度分布頻率圖Fig.3 Sentinel-1A radar satellite image subsidence rate accuracy map and accuracy distribution frequency chart in Guangzhou and Foshan

從廣佛地區(qū)沉降速率精度圖可以看出，研究區(qū)內(nèi)處理的最終結(jié)果精度大部分分布在0.3～0.6 mm/a范圍內(nèi)，標(biāo)準差(std)0.15 mm/a，平均精度(mean)0.51 mm/a，監(jiān)測結(jié)果穩(wěn)定可靠. 圖4為處理最終的沉降速率分布圖，紅色代表抬升，藍色代表沉降，從圖中廣佛整體的區(qū)域分布上研究，城區(qū)沉降量不大.

本次研究廣佛地區(qū)共生成6 833 697個PS點，生成的PS點沉降速率的標(biāo)準差(std)為2.22 mm/a，平均值(mean)為-0.06 mm/a，可以認為研究區(qū)域整體平穩(wěn)，局部發(fā)生沉降，與原先王華等對廣州佛山沉降研究結(jié)論整體相符[2].

4.2 東莞地區(qū)

東莞主城區(qū)地面沉降速率精度誤差及精度誤差統(tǒng)計如圖5所示

圖5 Sentinel-1A雷達衛(wèi)星影像沉降速率精度圖及精度分布頻率圖Fig.5 Sentinel-1A radar satellite image subsidence rate accuracy map and accuracy distribution frequency chart in Dongguan

從東莞地區(qū)沉降速率精度圖可以看出，研究區(qū)內(nèi)處理的最終結(jié)果精度大部分分布在0.3～0.7 mm/a范圍內(nèi)，標(biāo)準差(std)為0.14 mm/a，精度平均值(mean)為0.50 mm/a，可以認為此次處理的結(jié)果精度穩(wěn)定可靠. 此次實驗研究東莞地區(qū)共生成4 213 293個PS點，沉降速率圖如圖6所示，標(biāo)準(std)為1.56 mm/a，平均值(mean)為-0.05 mm/a，可以認為東莞研究區(qū)域整體平穩(wěn).

圖6 東莞地區(qū)沉降速率分布圖及速率頻率統(tǒng)計圖Fig.6 Subsidence rate distribution map and rate frequency chart of Dongguan

5 變形成因分析

本次研究重點分析了圖4廣佛地區(qū)沉降速率分布圖與圖6東莞地區(qū)沉降速率分布圖中所標(biāo)注的A、B、F、G、M、N、U及V點.

5.1 廣佛地區(qū)變形成因分析

據(jù)廣州日報報道，佛山五區(qū)范圍內(nèi)約有6 000～8 000口井抽取地下水，尤其在農(nóng)村和城鄉(xiāng)接合部地區(qū)地下水破壞尤為嚴重[10]. 此前有研究表明特別在工業(yè)密集區(qū)域及隧道施工區(qū)域，由于非法不合理的對地下水抽取，沉降現(xiàn)象更加明顯. 例如:2007年因武廣高鐵施工而導(dǎo)致金沙洲發(fā)生數(shù)起地面塌陷事故[11-12]，再如2008年7月21日，佛山市順德區(qū)北滘職業(yè)技術(shù)學(xué)校全校幾乎所有的建筑都出現(xiàn)了不同程度的沉降現(xiàn)象，并導(dǎo)致許多墻體開裂[2，12]. 廣州地下水抽取與北方地區(qū)例如北京、天津等城市因地下水開采造成大范圍沉降有所不同[13-14]. 北方地區(qū)特別是華北地區(qū)地下水開采幾乎都在深水層，而珠江三角洲地區(qū)地下水開采在潛水層. 潛水層可以通過降水或地表水下滲補給，由于地下水開采在廣佛主城區(qū)地區(qū)大部分都是在潛水層，距離地面5 m左右[3]. 所以廣佛地區(qū)包括周邊相鄰地區(qū)并未出現(xiàn)由于地下水開采而發(fā)生的大范圍地面沉降，只是出現(xiàn)的局部地面沉降.

圖4中A點標(biāo)注的季華西路沿線一帶在研究期間出現(xiàn)重大安全事故，2018年2月7日，該路段下方在建的佛山地鐵二號線綠湖島至湖涌盾構(gòu)期間發(fā)生地面塌陷事故，造成11人死亡、1人失蹤、8人受傷，直接經(jīng)濟損失達到5 323.8萬元[15]. 此處地陷事故是由于該區(qū)間地質(zhì)條件復(fù)雜，存在深厚富水粉砂層且臨近強透水的中粗砂層，且地下水具有承壓性. 據(jù)中央電視臺央廣網(wǎng)報道由于地鐵右線施工突發(fā)透水，導(dǎo)致軌道管片變形受損，外側(cè)土壤與管片應(yīng)力失衡，引發(fā)地面季華西路30多米路段坍塌[16]. A點在2015年12月至2018年5月沉降的時間序列圖，如圖7所示可以明顯看出季華西路一帶在研究期間出現(xiàn)沉降，最高累計沉降量達到52 mm，2015年12月至事故發(fā)生時已出現(xiàn)明顯沉降.

圖7 A點時間序列圖Fig.7 Point A time series chart

圖4中B點位于佛山樂從鎮(zhèn)大墩村，在建的佛山3號線經(jīng)過此地，并且附近有明顯關(guān)于地鐵的施工建筑. 主要沉降區(qū)域發(fā)生在嶺南大道南路與荷月路交口東約600 m處及岳步工業(yè)區(qū)附近一帶，其中嶺南大道南路與荷月路交口東約600 m處最大沉降速率為-23.17 mm/a沉降點附近交通繁雜，人口居住多，且沉降點西100 m處有工廠，此處安全風(fēng)險較高. 岳步工業(yè)區(qū)東北500 m處有明顯地鐵開挖口，附近施工活動較多. 據(jù)南方都市報報道，大墩站于2016年8月31日先行施工. 從圖8中可以看出在2016年9月(黑色五角星標(biāo)注)左右出現(xiàn)沉降速率加快，此后的沉降趨勢基本上為與施工周期保持一致. 2017年3月16日(黑色正方形標(biāo)注)至2019年12月31日大墩村建設(shè)工程實施封閉施工，從圖8可以看出從2017年3月中旬開始，沉降速率進一步加快. 可以認定B點沉降原因主要在于佛山3號線施工引起，由于地鐵施工破壞土壤穩(wěn)定，造成該處發(fā)生沉降. 圖8為B點附近一沉降點的累計沉降量達到48 mm.

圖8 B點時間序列圖Fig.8 Point B time series chart

圖4中F點位于佛山鋼鐵世界至黎湖工業(yè)區(qū)一帶，沉降區(qū)域內(nèi)鋼鐵廠密集分布，物流公司分布密集，還分布著眾多家具廠，監(jiān)測獲得的最高沉降速率為-24.8 mm/a. F點附近兩沉降點累計沉降量如圖9所示，達到了46.36 mm. G點與F點類似，位于密集的工業(yè)區(qū)內(nèi)，周圍皮革廠眾多，人口密集；G點最高沉降速率達到-25 mm/a，沉降量達到48.86 mm. 從G點與F點位置上分析，此處的沉降點與工業(yè)分布十分吻合，初步認定為工業(yè)不合理開采地下水及當(dāng)?shù)厝藛T抽取地下水引起的地面沉降. NG等利用廣佛地區(qū)2011年至2017年COSMO-SkyMed衛(wèi)星數(shù)據(jù)做時序InSAR分析，也顯示此處周圍有沉降，但未具體分析[3].F點與G點周圍部分沉降點累計沉降量如圖9所示.

5.2 東莞地區(qū)變形成因分析

圖9 F、G點時間序列圖及工廠分布圖Fig.9 Point F, G time series chart and factory distribution map

東莞主城區(qū)的沉降速率圖如圖6所示，整體平穩(wěn)，局部發(fā)現(xiàn)地面沉降，范圍較小，只是呈現(xiàn)出點狀分布. 近年來東莞市經(jīng)濟飛速發(fā)展，有世界工廠之稱.如圖11所示，研究區(qū)內(nèi)的點狀分布的沉降點也主要集中在工廠密集區(qū)域內(nèi). M點、N點分別位于東莞市與深圳市、廣州市接壤處，工業(yè)密集，有各種印刷廠、涂料廠、造紙廠、五金廠等等. 此前也有報道證明目前深圳關(guān)外地區(qū)(東莞)的工廠以前開挖了大量的非法地下水水井，目前很多水井還在運行，無序開采將使地下水環(huán)境進一步惡化，可能會導(dǎo)致地質(zhì)災(zāi)害和海水倒灌[17]. 這兩塊區(qū)域中發(fā)現(xiàn)的沉降點初步歸因于工廠及人口對地下水不合理開采以及人類活動所致. M點與N點區(qū)域的明顯沉降點時間序列圖如圖10所示，點M與點N的累積沉降量分別為達到39 mm、38 mm.

圖10 M、N時間序列圖Fig.10 Point M and N time series chart

M、N、U與V點周圍的工業(yè)分布如圖11所示，圖中紅色數(shù)字代表的是工業(yè)聚集點，紅色圓圈點代表大工業(yè)廠房分布地. 從圖11中與圖6東莞沉降圖可以看出，工業(yè)聚集點M、N、U與V都有明顯沉降，這同時也印證了文獻[3， 7]中工廠及當(dāng)?shù)厝藛T不合理抽取地下水導(dǎo)致地面沉降的觀點.

圖11 M、N、U、V點工廠分布圖Fig.11 Industrial distribution around the M, N, U, and V points

U點沉降明顯區(qū)域位于光正實驗學(xué)校周圍至茶山塘工業(yè)區(qū)一帶，學(xué)校周圍有明顯的施工建筑，施工活動頻繁. 工業(yè)區(qū)內(nèi)密集分布著制衣作坊、紙廠及印刷廠等. V點位于東莞市田坑工業(yè)區(qū)，周圍有村尾第一工業(yè)區(qū)、四海工業(yè)園等，工業(yè)區(qū)內(nèi)密集分布著模具廠、塑膠廠、五金廠等各種耗水量巨大的工廠. 此外，此處人口密集，大量人口抽取地下水量也是巨大的，與M點、N點區(qū)域極為相似. 從圖12可以看出，V點工業(yè)區(qū)沉降速率更快，超過-18 mm/a，累計沉降量也更為明顯，U點與V點累積沉降量分別達到45 mm、52 mm.

圖12 U、V時間序列圖Fig.12 Point U and V time series chart

6 結(jié)論

本文利用Sentinel-1A數(shù)據(jù)研究了廣州、佛山與東莞地區(qū)的地面變形情況，結(jié)果表明2015年12月至2018年6月期間，整體區(qū)域平穩(wěn)，局部發(fā)現(xiàn)明顯沉降.將此次研究的廣佛地區(qū)與其他關(guān)于廣佛地區(qū)的研究進行對比發(fā)現(xiàn)，此次研究結(jié)果可靠. 局部地面變形區(qū)域主要集中在工程施工區(qū)與工廠密集區(qū). 本次研究的主要的幾個危險性較高的沉降區(qū)域中，A點(季華西路)發(fā)生了重大事故，B點(大墩村附近)暫未發(fā)生事故，但不排除未來有事故發(fā)生的可能. 希望并建議有關(guān)部門做好風(fēng)險防控措施，預(yù)防進一步的災(zāi)害發(fā)生.其中此次試驗研究發(fā)現(xiàn)F點(佛山市六涌工業(yè)區(qū)與鋼鐵世界一帶)及G點(佛山市亞洲國際家具材料交易中心周圍一帶)沉降明顯，建議有關(guān)部門采取相關(guān)措施進一步查明沉降原因. 東莞市U點(光正實驗學(xué)校至茶山塘邊工業(yè)區(qū)一帶)、V點(田坑工業(yè)區(qū))、M點及N點周圍沉降明顯，希望有關(guān)部門進一步查明原因，給予防控措施.

總體來說，本次實驗結(jié)果可靠，證明了利用Sentinel-1A也可以較好地監(jiān)測出珠江三角洲所研究城市的地面變形情況，因本次研究缺少珠三角地區(qū)的水文資料、水準數(shù)據(jù)與后續(xù)地鐵沿線施工情況數(shù)據(jù)，后續(xù)研究中將繼續(xù)結(jié)合不同資料進一步分析珠江三角洲地區(qū)的沉降情況.

致謝：感謝歐洲空間局Sentinel-1A衛(wèi)星數(shù)據(jù)的免費使用.