高 輝，羅孝文*，吳自銀，陽凡林

(1.山東科技大學(xué) 測繪科學(xué)與工程學(xué)院，山東 青島 266590；2.自然資源部 第二海洋研究所，浙江 杭州 310012；3.自然資源部 海底科學(xué)重點實驗室, 浙江 杭州 310012)

0 引言

珠江三角洲位于中國廣東省中南部，毗鄰港澳，區(qū)域內(nèi)城市密集，經(jīng)濟發(fā)達，有“南海明珠”之稱。河口等沿海地區(qū)同時受海洋、陸地及河流的相互作用，生態(tài)環(huán)境非常脆弱[1]。近30 a來珠江口發(fā)展迅速，大量基礎(chǔ)設(shè)施建設(shè)、地下資源開采、軟土固結(jié)壓實等人類活動和自然因素的影響，使該區(qū)域發(fā)生大面積的地面沉降，帶來地面塌陷、基礎(chǔ)設(shè)施破壞等地質(zhì)災(zāi)害，全球海平面上升[2]加劇了咸潮上溯的影響及風(fēng)暴潮災(zāi)害的風(fēng)險[3-4]，愈發(fā)威脅該區(qū)域的生態(tài)環(huán)境及人類的生產(chǎn)生活。因此，及時準確地揭示珠江河口地區(qū)地面沉降的時空分布和形變趨勢具有重要的意義。

合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar，InSAR)技術(shù)具有全天時、全天候、覆蓋范圍廣等特點，能夠獲取高分辨率、高精度的地形及地表形變信息[5-6]。近年來，以永久散射體技術(shù)(Persistent Scatterers，PS)[7]和小基線集技術(shù)(Small BAseline Subset，SBAS)[8]為代表的時序InSAR技術(shù)已廣泛應(yīng)用于河口海岸帶地區(qū)的地面沉降監(jiān)測研究。劉勇[9]、張金芝 等[10-11]及劉一霖 等[12]分別用不同的方法研究了黃河口地區(qū)的地面沉降；JIANG et al[13-14]利用SAR干涉點目標分析法(IPTA)監(jiān)測香港國際機場及澳門填海區(qū)的地面沉降，并揭示了填海與沉降的關(guān)系；WANG et al[15]利用多軌SAR數(shù)據(jù)采用π-rate方法監(jiān)測了珠江三角洲深圳、東莞地區(qū)的地面沉降，得到沿岸500 m內(nèi)年沉降速率為5 mm/a；ERBAN et al[16]探測了湄公河三角洲的地面沉降，并認為地下水抽取是引起大范圍沉降的原因。

這些時序InSAR方法對監(jiān)測點目標的選取十分嚴格，要求在所涉及的干涉紋圖中都保持較高的相干性[17]。而河口海岸帶地區(qū)植被覆蓋密，人工目標稀疏，大氣擾動相對復(fù)雜，難以準確反映地面沉降的時空變化[18-19]，限制了這些方法在河口海岸帶地面沉降檢測中的應(yīng)用。

FERRETTI et al[20]提出SqueeSAR技術(shù)，將同分布散射體(Distributed Scatterer，DS)融入到傳統(tǒng)時序InSAR技術(shù)中，大大提高了有效監(jiān)測點的密度。但該方法求解非線性最優(yōu)相位需要迭代，影響計算效率。針對上述問題，HAN et al[21]提出一種基于特征分解估算DS候選點最優(yōu)相位的方法，該方法大大提高了計算效率，并通過杭州灣沿岸上虞區(qū)的應(yīng)用驗證了在沿海地區(qū)的有效性。鑒于目前關(guān)于珠江口地區(qū)的地面沉降研究都在小區(qū)域范圍內(nèi)[15,22-23]，隨著Sentinel-1A數(shù)據(jù)的免費公開，本文采用HAN et al[21]提出的方法，利用2015—2018年間67景Sentinel-1A SAR影像對珠江口地區(qū)進行大范圍的地面沉降監(jiān)測，獲得了該區(qū)域的地面沉降信息及時空分布，并對幾處沉降嚴重地區(qū)進行了分析。

1 技術(shù)方法

采用基于Kolmogorov-Smirnov(KS)檢驗和特征分解的DS-InSAR時序分析方法，其基本數(shù)據(jù)處理如下[21]：首先將Single Look Complex(SLC)影像配準至主影像，采用KS檢驗算法選擇穩(wěn)定的DS候選點，利用DS候選點目標的相位信息進行空間自適應(yīng)濾波并構(gòu)建協(xié)方差矩陣；然后通過特征分解方法估算最優(yōu)相位并代替DS目標原始相位；最后根據(jù)擬合優(yōu)度檢測確定最終的DS點目標，聯(lián)合PS點進行時序InSAR處理，經(jīng)過地形誤差改正、大氣誤差改正、相位解纏等步驟得到最終的形變速率和形變時間序列。下面重點介紹DS點選取過程中基于特征分解方法估算最優(yōu)相位的關(guān)鍵步驟。

基于特征分解的最優(yōu)相位估計方法最初由FORNARO et al[24]提出，通過從協(xié)方差矩陣或相干矩陣中提取與最大特征值相關(guān)的特征向量的相位估計最優(yōu)相位。對通過KS檢驗算法選取并進行空間自適應(yīng)濾波后構(gòu)建的相干矩陣進行特征分解，得：

(1)

其中：T為候選點的最優(yōu)相位，N為干涉相位，λi為特征值且按降序排序，Ui為特征值對應(yīng)的特征向量。假設(shè)U1為最大特征分量，可表示為最大值約束問題，將最大特征值對應(yīng)的特征向量的相位作為DS候選點的最優(yōu)相位[25]：

(2)

約束條件：

(3)

2 研究區(qū)概況與實驗數(shù)據(jù)

2.1 研究區(qū)概況

珠江三角洲位于珠江入?？谔帲晌鹘?、北江和東江入海時沉積而成，地理坐標為112°30′～114°15′E， 21°17′～23°30′N。珠江三角洲呈倒置三角形，地貌類型以三角洲平原為主體，三角洲前緣及西北部有低山、丘陵分布，中部則綴有零星低丘及臺地，平原上河道汊港縱橫交錯，水塘澤地星羅棋布。地層以第四紀地層分布最廣，大面積展布于三角洲平原，巖性以松軟土層沉積為主要特征。南部為海岸帶，是區(qū)內(nèi)海陸相互作用最強烈的地帶，具有海岸線長、海岸曲折、島嶼眾多等特點，但軟土厚度較大，容易發(fā)生地面沉降地質(zhì)災(zāi)害。珠江三角洲地處亞熱帶季風(fēng)氣候區(qū)，溫暖潮濕，雨量充沛。近年來人類活動的影響和自然條件的變化使該區(qū)域發(fā)生大面積地面沉降，嚴重危害該地區(qū)的經(jīng)濟社會可持續(xù)發(fā)展。本文選擇珠江口地面沉降較為嚴重的地區(qū)作為監(jiān)測區(qū)，面積約3 000 km2，有效覆蓋了珠江口地區(qū)。

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

使用67景Sentinel-1A SAR影像進行數(shù)據(jù)處理，Sentinel-1A衛(wèi)星由歐空局發(fā)射，具有短重復(fù)周期、大幅寬、多極化等特點。數(shù)據(jù)獲取時間跨度為2015年6月—2018年3月，其采用干涉寬幅(IW)數(shù)據(jù)成像模式，極化方式為VV極化，所處波段為C波段，具有250 km的掃描寬度和5 m×20 m的地面分辨率。同時采用美國航空航天局(NASA)提供的SRTM 30 m DEM數(shù)據(jù)，去除地形相位與DEM誤差。

3 監(jiān)測結(jié)果

3.1 InSAR地面沉降結(jié)果

利用基于KS檢驗和特征分解的DS-InSAR時序分析方法，處理SAR影像數(shù)據(jù)，設(shè)定時間基線100 d和空間基線200 m的閾值作為小基線集干涉圖選取限制，生成了178對干涉組合，其組合如圖1所示。采用最小費用流法(minimum cost flow, MCF)進行相位解纏。通過數(shù)據(jù)處理，獲得了珠江口地區(qū)2015年6月—2018年3月地面沉降的年均形變速率(圖2)和累計形變分布(圖3)。其中，負值代表下降，正值代表抬升。從圖中可以看出，監(jiān)測區(qū)時序InSAR獲取的地面形變結(jié)果總體呈現(xiàn)地面沉降趨勢，且沉降在空間上分布不均勻，部分區(qū)域存在明顯沉降現(xiàn)象，抬升區(qū)范圍較小。珠江口西岸沉降相對東岸分布更廣，沉降量更大。珠江三角洲西北部(廣州、佛山)和東南部(東莞南部)為主要沉降區(qū)，其中監(jiān)測區(qū)西北部佛山禪城區(qū)南莊鎮(zhèn)及順德區(qū)西部龍江鎮(zhèn)、東莞市南部長安鎮(zhèn)局部存在嚴重的沉降現(xiàn)象，最大沉降速率可達25 mm/a。監(jiān)測穩(wěn)定區(qū)主要集中在東莞市西部虎門鎮(zhèn)附近，該區(qū)域沒有明顯的下沉現(xiàn)象。

圖1 短基線集干涉組合Fig.1 Combination of small-baseline interferograms

圖2 珠江口地區(qū)InSAR地面年均形變速率Fig.2 Annual mean deformation rate in the Pearl River Estuary from InSAR time series

圖3 珠江口地區(qū)InSAR累積形變量Fig.3 Cumulative deformation in the Pearl River Estuary from InSAR time series

通過DS-InSAR方法在實驗區(qū)獲得的DS目標分布，可以看出在非城市區(qū)域有著較高密度的高相干點，主要集中在瀝青路面、裸土及其他有DS點特性的地物上，有效增強了高質(zhì)量InSAR監(jiān)測點的空間密度，進而提高了時序InSAR形變反演的精度。

3.2 沉降嚴重地區(qū)分析

3.2.1 東莞市長安鎮(zhèn)

長安鎮(zhèn)位于東莞市南端，東鄰深圳市，南臨珠江口，西連虎門港，是地面沉降信號明顯的重點區(qū)域。由圖4可以看出，長安鎮(zhèn)整體都有明顯沉降信號，沉降速率在8～13 mm/a范圍內(nèi)。長安鎮(zhèn)沉降由西北(靠近虎門)向東南表現(xiàn)為沉降加劇，尤其是靠近珠江口區(qū)域，包括南部深圳市寶安區(qū)沙井街道，其最大沉降速率達24.1 mm/a，3 a累積最大沉降量可達45.3 mm。近年來東莞市經(jīng)濟飛速發(fā)展，隨之而進行的大量填?；顒?、基礎(chǔ)設(shè)施建設(shè)、工業(yè)科技園區(qū)建設(shè)是該地區(qū)發(fā)生沉降的主要原因。

3.2.2 廣州市白云區(qū)

廣州市白云區(qū)位于廣州市城區(qū)北部，東鄰增城區(qū)，西界佛山市南海區(qū)，南連荔灣區(qū)、越秀區(qū)、天河區(qū)、黃埔區(qū)4城區(qū)，北接花都區(qū)和從化區(qū)。由圖5可以看出，白云區(qū)沙貝立交附近存在嚴重不均勻分布的地面沉降現(xiàn)象，年均沉降速率為8～15 mm/a；尤其是沙貝立交南部以及東部沙涌村下沙附近，地面沉降中心最大沉降速率可達15.2 mm/a，3 a最大累計沉降量達40 mm；沙貝立交北部為相對穩(wěn)定區(qū)域，年均沉降速率大部分為3～8 mm/a，僅有零星沉降分布。廣州白云區(qū)位于廣花盆地，其地質(zhì)環(huán)境脆弱，相關(guān)調(diào)查結(jié)果表明，基坑工程、地鐵工程、隧道工程等地下空間的開發(fā)引起的巖溶地面塌陷是該地區(qū)地面沉降的主要原因。

3.2.3 佛山市南海區(qū)

佛山市南海區(qū)地處珠江三角洲腹地，東連廣州市芳村區(qū)、番禺區(qū)，南接順德區(qū)、鶴山市、新會區(qū)，西鄰三水區(qū)、高明區(qū)，北鄰花都區(qū)，環(huán)抱佛山禪城區(qū)。由圖6可知對江沙北部整體較為穩(wěn)定，部分地區(qū)存在地面沉降現(xiàn)象，平均沉降速率為3～8 mm/a，最大沉降中心在海南立交南部，沉降速率為10.2 mm/a；對江沙南部存在3個沉降速率大于15 mm/a的沉降中心，從西至東分別位于平勝大橋東部，新五斗大橋南部，三江郊野公園東南部。其中三江郊野公園東南部最為嚴重，沉降范圍最大，最大沉降速率達23.9 mm/a，3 a累積沉降量最大可達65 mm。近年來，南海區(qū)的經(jīng)濟發(fā)展迅速，大規(guī)模城市、地鐵建設(shè)和工業(yè)園開發(fā)以及區(qū)內(nèi)廣佛江珠、廣珠西線高速公路等荷載量大增，為該地區(qū)地面沉降的主要原因。

圖4 東莞市長安鎮(zhèn)地面年均形變速率Fig.4 Annual mean deformation rate in Chang’an Town, Dongguan

圖5 廣州市白云區(qū)地面年均形變速率Fig.5 Annual mean deformation rate in Baiyun District, Guangzhou

圖6 佛山市南海區(qū)地面年均形變速率 Fig. 6 Annual mean deformation rate in Nanhai District, Foshan

4 小結(jié)

利用基于KS檢驗和特征分解的DS-InSAR方法，對67景Sentinel-1A影像數(shù)據(jù)進行時序分析，獲得了2015—2018年間珠江口地區(qū)地面沉降場的時空演變過程，重點分析了3處沉降較明顯區(qū)域的沉降時空變化特征。采用的DSI方法在瀝青路面、裸土等非城市區(qū)域提取了較大密度的高相干點。研究結(jié)果表明，珠江口地區(qū)總體呈地面沉降趨勢，珠江三角洲西北部(廣州、佛山)和東南部(東莞南部)為主要沉降區(qū)，最大沉降速率可達25 mm/a；珠江口西岸沉降相對東岸分布更廣，沉降量更大；東莞市西部虎門鎮(zhèn)附近比較穩(wěn)定，該區(qū)域沒有明顯的下沉現(xiàn)象。

本文監(jiān)測了珠江口地區(qū)的地面沉降，后續(xù)將結(jié)合海岸線、河岸線的變化，基礎(chǔ)設(shè)施建設(shè)，地下水抽取及填海造陸情況，對該地區(qū)地面沉降的因素進行詳細的分析，并研究其和海岸線變化的關(guān)系，進一步為地面沉降的預(yù)防和治理提供參考。