高 輝，羅孝文*，吳自銀，陽(yáng)凡林

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

0 引言

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

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

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

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

1 技術(shù)方法

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

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

(1)

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

(2)

約束條件：

(3)

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

2.1 研究區(qū)概況

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

2.2 實(shí)驗(yàn)數(shù)據(jù)

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

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

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

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

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

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

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

長(zhǎng)安鎮(zhèn)位于東莞市南端，東鄰深圳市，南臨珠江口，西連虎門(mén)港，是地面沉降信號(hào)明顯的重點(diǎn)區(qū)域。由圖4可以看出，長(zhǎng)安鎮(zhèn)整體都有明顯沉降信號(hào)，沉降速率在8～13 mm/a范圍內(nèi)。長(zhǎng)安鎮(zhèn)沉降由西北(靠近虎門(mén))向東南表現(xiàn)為沉降加劇，尤其是靠近珠江口區(qū)域，包括南部深圳市寶安區(qū)沙井街道，其最大沉降速率達(dá)24.1 mm/a，3 a累積最大沉降量可達(dá)45.3 mm。近年來(lái)東莞市經(jīng)濟(jì)飛速發(fā)展，隨之而進(jìn)行的大量填?；顒?dòng)、基礎(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ū)沙貝立交附近存在嚴(yán)重不均勻分布的地面沉降現(xiàn)象，年均沉降速率為8～15 mm/a；尤其是沙貝立交南部以及東部沙涌村下沙附近，地面沉降中心最大沉降速率可達(dá)15.2 mm/a，3 a最大累計(jì)沉降量達(dá)40 mm；沙貝立交北部為相對(duì)穩(wěn)定區(qū)域，年均沉降速率大部分為3～8 mm/a，僅有零星沉降分布。廣州白云區(qū)位于廣花盆地，其地質(zhì)環(huán)境脆弱，相關(guān)調(diào)查結(jié)果表明，基坑工程、地鐵工程、隧道工程等地下空間的開(kāi)發(fā)引起的巖溶地面塌陷是該地區(qū)地面沉降的主要原因。

3.2.3 佛山市南海區(qū)

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

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

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