楊亞彬，李 慧

(1.廣東省地質(zhì)測繪院，廣東 廣州 510800；2．廣東工業(yè)大學(xué)，廣東 廣州 510006)

珠海是廣東省的地級市，是珠江三角洲中心城市之一，地處北緯21°48′～22°27′，東經(jīng)113°03′～114°19′之間。位于廣東省珠江口的西南部，東與香港隔海相望，南與澳門相連，西鄰江門市，北與中山市接壤，是珠三角中海洋面積最大、島嶼最多、海岸線最長的城市，素有“百島之市”之稱[1]。珠海市的城市規(guī)模在過去的30年中得到快速發(fā)展，以此同時在大量而且頻繁的基礎(chǔ)城市建設(shè)中由于土建工程而抽取地下水、基坑開挖等人為因素必然會對土層產(chǎn)生擾動。而珠海地處珠江三角洲，主要城區(qū)的地質(zhì)結(jié)構(gòu)比較不穩(wěn)定[2]，地下水分布廣泛，地質(zhì)災(zāi)害頻繁發(fā)生，尤其是由于地表沉降而發(fā)生的災(zāi)害已經(jīng)給社會的人身以及財產(chǎn)安全帶來巨大的隱患。因此急需一種手段能夠?qū)Υ蠓秶牡乇沓两颠M(jìn)行監(jiān)測，為地質(zhì)災(zāi)害的防治工作提供科學(xué)依據(jù)。

合成孔徑雷達(dá)干涉測量(InSAR)技術(shù)是一種新型的空間對地觀測技術(shù)，具有高時空分辨率、高精度、監(jiān)測范圍廣、穿透性強(qiáng)等優(yōu)勢，是當(dāng)今地學(xué)界較為先進(jìn)的基于面觀測的地表形變監(jiān)測手段[3-5]。目前，已有眾多的國內(nèi)外專家學(xué)者成功應(yīng)用InSAR技術(shù)進(jìn)行地表沉降監(jiān)測工作，并取得讓人矚目的成果[6-11]。本研究基于Sentinel-1A衛(wèi)星數(shù)據(jù)[12]，采用小基線集干涉技術(shù)(SBAS)對珠海市相關(guān)區(qū)域開展地表形變監(jiān)測，分析沉降區(qū)域的時空分布規(guī)律和引起地表沉降的影響因素，為港珠澳大灣區(qū)城市基礎(chǔ)設(shè)施健康診斷提供評價依據(jù)，提高城市災(zāi)害風(fēng)險管控能力。

1 SBAS技術(shù)原理及其數(shù)據(jù)源

短基線集( Small-Baseline Subset，SBAS) 技術(shù)是近十幾年發(fā)展起來的一種高級DInSAR數(shù)據(jù)處理方法，最先是2002年由Berardino和Lanari等學(xué)者提出[13-14]，主要用于研究低分辨率、大尺度上的形變。SBAS方法通過自由組合基線較短的影像對，產(chǎn)生一系列基于不同主影像的時間序列干涉圖子集，再利用矩陣的奇異值分解( SVD)方法，將多個短基線集聯(lián)合起來求解，解決各數(shù)據(jù)集之間空間基線過長造成的時間不連續(xù)問題，提高監(jiān)測的時間分辨率，得到覆蓋整個觀測時間的形變序列和平均沉降速率。

本項目擬對斗門區(qū)、金灣區(qū)、香洲區(qū)以及澳門特別行政區(qū)進(jìn)行監(jiān)測，并對監(jiān)測范圍進(jìn)行適當(dāng)外延，監(jiān)測總面積約8 192 km2，監(jiān)測范圍如圖1中白框所示。

圖1 研究區(qū)范圍

本項目采用32景Sentine-1A SAR影像進(jìn)行數(shù)據(jù)處理。Sentinel-1是專為“哥白尼計劃”設(shè)計的第一顆衛(wèi)星，主要用以監(jiān)測歐洲、加拿大和極地區(qū)域的環(huán)境變化，如海冰、海洋環(huán)境以及地表運(yùn)動，地表森林、水源和土壤動態(tài)變化。Sentinel-1衛(wèi)星是歐洲極地軌道C波段雷達(dá)成像系統(tǒng)，采用預(yù)編程無沖突操作模式，可實現(xiàn)高分辨率全球陸地、海岸帶和航線的成像[13-14]。其相對于ERS-1/1/SAR和ENVISAT/ASAR具有重訪周期短、單景影像覆蓋范圍高的優(yōu)勢，雙星座衛(wèi)星的Sentinel-1重復(fù)周期可縮短至6 d。本研究時間跨度為2016-01-29—2017-11-07，數(shù)據(jù)成像模式為IW(干涉寬幅)，幅寬為250 km，地面分辨率為5 m×20 m，衛(wèi)星影像數(shù)據(jù)信息如表1所示。

數(shù)字高程模型(Digital Elevation Model，DEM) 采用日本宇宙航空開發(fā)機(jī)構(gòu)(JAXA)聯(lián)合美國RESTEC公司利用35 km×35 km的光學(xué)影像生產(chǎn)的30 m分辨率的高精度DEM，該光學(xué)影像采集自日本發(fā)射的ALOS衛(wèi)星所搭載的PRISM(全色立體測繪儀)傳感器，由此獲得的DEM平面及高程精度均可達(dá)到5 m，其數(shù)據(jù)產(chǎn)品能在全球范圍內(nèi)覆蓋。相較于SRTM數(shù)據(jù)，AW3D30的數(shù)字高程產(chǎn)品由于不存在雷達(dá)信號失相干的問題，因此不會出現(xiàn)數(shù)據(jù)空洞的問題。利用該DEM進(jìn)行去除地形處理，能夠精確提供研究區(qū)域的高程數(shù)據(jù)，并且能有效降低由于外部提供的DEM不正確而產(chǎn)生的干涉去平地效應(yīng)出錯等情況，使得InSAR結(jié)果的可靠性得以提高。

表1 SAR衛(wèi)星影像數(shù)據(jù)信息

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

2.1 差分干涉圖的生成

根據(jù)設(shè)定的時空基線閾值生成短基線集，對短基線集進(jìn)行差分干涉處理。與PSI方法相比，SBAS避免了長空間基線的幾何去相干問題和長時間基線的時間去相干問題，提高干涉圖的相干性和時間分辨率，確保地表形變監(jiān)測的準(zhǔn)確性。本項目在SBAS數(shù)據(jù)處理過程中，時間閾值設(shè)定為200 d，垂直基線設(shè)定為45%的臨界基線，共生成239對干涉像對。剔除解纏效果不理想或出現(xiàn)明顯大氣效應(yīng)的干涉像對，使其不參與后續(xù)的處理工作。經(jīng)過上述處理步驟，SAR數(shù)據(jù)集中用于反演的干涉數(shù)據(jù)集共包含92對干涉像對，形成的連接圖如圖2所示。

圖2 SAR干涉對連接圖

2.2 高相干目標(biāo)點的選取

SBAS技術(shù)處理的對象是一些高相干目標(biāo)點，這些點的選擇是通過對一系列雷達(dá)干涉影像進(jìn)行相關(guān)閾值處理獲取。目前相干目標(biāo)點的選取方法主要有：基于像元強(qiáng)度的穩(wěn)定性、基于相位穩(wěn)定性和基于空間相干性等方法。強(qiáng)度穩(wěn)定性方法是根據(jù)像元強(qiáng)度的穩(wěn)定性來判斷相位穩(wěn)定性，該方法需要建立在大量數(shù)據(jù)進(jìn)行統(tǒng)計分析的基礎(chǔ)上，因此要求數(shù)據(jù)量一般不少于30景SAR影像。相位穩(wěn)定性分析的方法根據(jù)相位的時間和空間相干性，對相位殘差(噪聲)進(jìn)行分析的技術(shù)，該方法涉及到復(fù)雜的數(shù)學(xué)方法，而且需采用大量模擬實驗來確定閾值。最后綜合考慮算法的復(fù)雜度并結(jié)合本研究數(shù)據(jù)量小的特點，本研究采用相干性的方法來選取高相干目標(biāo)點。對時間序列的相干圖取均值，利用累積影像的平均相干值來選取相干目標(biāo)。

2.3 高相干點沉降速率和形變場的獲取

在選取的相干目標(biāo)上建立Delaunay三角網(wǎng)，以三角網(wǎng)上的所有邊為單元進(jìn)行建模，建立觀測方程。采用奇異值分解的方法進(jìn)行求解，即可得到各時間區(qū)間的沉降速率。各時段沉降速率在時間域上進(jìn)行積分即可得到各個時間點的累積形變量。

3 結(jié)果分析

本項目通過對32景Sentinel-1A衛(wèi)星IW模式的SAR影像進(jìn)行數(shù)據(jù)處理，得到研究區(qū)的年平均沉降速率圖如圖3所示。從圖3中可看出該監(jiān)測區(qū)沉降速率最大可達(dá)-101.05 mm/a，抬升速率最大達(dá)到75.25 mm/a。監(jiān)測區(qū)南部沿海區(qū)域存在輕微抬升現(xiàn)象，抬升速率處于17.44 mm/a以內(nèi)。珠海市斗門區(qū)大部分地區(qū)處于穩(wěn)定狀態(tài)，年平均沉降速率不超過5 mm/a，但斗門區(qū)白蕉鎮(zhèn)地區(qū)沉降較為嚴(yán)重，最大年平均沉降速率可達(dá)-101.05 mm/a。珠海市金灣區(qū)主要城區(qū)年平均沉降速率基本處于-5～5 mm/a之間，城區(qū)處于穩(wěn)定狀態(tài)，但在高欄港站、平灣四路附近等區(qū)域發(fā)現(xiàn)較為明顯的沉降信號，沉降速率最大達(dá)到-75 mm/a。香洲區(qū)大部分地區(qū)處于年平均形變速率在-17～17 mm/a之間，較為穩(wěn)定，但橫琴等區(qū)域存在明顯的沉降漏斗。澳門大部分區(qū)域的年平均形變速率處于5.88～17.44mm/a之間，少數(shù)地區(qū)存在沉降信號，最大沉降量達(dá)到-40.36 mm/a。

圖3 整體平均沉降速率圖

3.1 斗門區(qū)地表形變結(jié)果

將沉降速率圖疊加在Google Earth上，可以看出珠海市斗門區(qū)大部分地區(qū)年平均沉降速率在5 mm/a以內(nèi)，但斗門區(qū)白蕉鎮(zhèn)地區(qū)出現(xiàn)嚴(yán)重的沉降現(xiàn)象，如圖4所示，最大沉降速率達(dá)到101.05 mm/a，圖中白色方框即為沉降明顯區(qū)域。對此區(qū)域選定采樣點#1、#2、#3、#4進(jìn)行形變量時間序列提取得到圖5，從圖中可看出，采樣點均呈現(xiàn)一致的下沉趨勢，且累計沉降量可以達(dá)到將近-132 mm，沉降趨勢明顯。

圖4 斗門區(qū)白蕉鎮(zhèn)沉降速率分布

圖5 斗門區(qū)白蕉鎮(zhèn)相干目標(biāo)采樣點時間序列

3.2 金灣區(qū)地表形變結(jié)果

由平均沉降速率圖3可以看出，珠海市金灣區(qū)主城區(qū)處于穩(wěn)定狀態(tài)，年平均形變速率不超過5 mm/a，南部沿海區(qū)域出現(xiàn)20 mm以內(nèi)的輕微抬升現(xiàn)象，金灣區(qū)整體較為穩(wěn)定。其中在紅旗鎮(zhèn)、三灶鎮(zhèn)、南水鎮(zhèn)及平沙鎮(zhèn)共發(fā)現(xiàn)8處地方存在較為明顯的沉降信號。在沉降信號明顯區(qū)域抽取有代表性的相干目標(biāo)點，提取其垂直形變的時間序列進(jìn)行分析。

圖6 海洋石油工程有限公司附近區(qū)域沉降速率分布

如圖6所示，在使用SAR衛(wèi)星的IW數(shù)據(jù)進(jìn)行SBAS處理得到的結(jié)果中，可以在海洋石油工程有限公司廠區(qū)觀測到較為明顯的沉降信號，該區(qū)域相干目標(biāo)點的年平均沉降速率可達(dá)到-75.04～-63.48 mm/a。對相干目標(biāo)點進(jìn)行時間序列分析，從圖7可以看到，位于沉降速率較慢地區(qū)的#1和#3點累計沉降量最大值為-44 mm，位于平均沉降速率達(dá)到-75.04 mm/a的#2點累計沉降量達(dá)到-80 mm。由圖8可看出，該地區(qū)為海洋石油工程有限公司廠區(qū)，廠區(qū)內(nèi)存在石油開采等作業(yè)項目，這可能是導(dǎo)致該地區(qū)地面發(fā)生沉降的主要因素。

圖7 海洋石油工程有限公司附近區(qū)域相干目標(biāo)采樣點時間序列

圖8 海洋石油工程有限公司石油開采作業(yè)現(xiàn)場

3.3 香洲區(qū)地表形變結(jié)果

從圖3年平均沉降速率圖可以看出香洲區(qū)從北向南呈現(xiàn)出先沉降后抬升的趨勢，平均形變速率在-17～17 mm/a之間，大部分地區(qū)處于穩(wěn)定狀態(tài)，但在橫琴島等區(qū)域發(fā)現(xiàn)較為明顯的沉降信號。在香洲區(qū)發(fā)現(xiàn)的沉降信號中，最大沉降速率達(dá)到-63.48 mm/a。將有明顯沉降現(xiàn)象的區(qū)域選取采樣點，提取該像元形變值的時間序列。

圖9 港珠澳大橋碼頭沉降速率分布

在香洲區(qū)東北部的港珠澳大橋碼頭發(fā)現(xiàn)小范圍的沉降漏斗，平均下沉速率最大達(dá)到-63.48 mm/a(見圖9)。在港珠澳大橋1號碼頭和2號碼頭選取采樣點，提取其時間序列，從圖10可以看出，采樣點#1、#2、#3點表現(xiàn)出一致的下沉現(xiàn)象和下沉量級，累積沉降量可達(dá)到-80 mm，并且有繼續(xù)下沉的趨勢。經(jīng)過相關(guān)資料收集和分析，港珠澳大橋碼頭正在進(jìn)行人工島及相關(guān)工程施工，而這很可能是導(dǎo)致該地區(qū)出現(xiàn)沉降現(xiàn)象的原因。

圖10 港珠澳大橋碼頭相干目標(biāo)采樣點時間序列

圖11 在珠海北站及附近區(qū)域沉降速率分布

根據(jù)年平均沉降速率圖11可以看出，在珠海北站及附近區(qū)域發(fā)現(xiàn)沉降信號，最大沉降達(dá)到-63.48 mm/a，其中興中路附近區(qū)域沉降較為明顯。對沉降明顯區(qū)域選取采樣點#4、#5、#6、#7點，提取采樣點的形變時間序列得到圖12，可以看出位于沿海區(qū)的#6和#7點平均沉降速率較快，累計沉降量達(dá)到-100 mm,#4和#5點累計沉降量近-60 mm，整個區(qū)域累計沉降量相對其它區(qū)域大，并且沉降速率較快。由TOD項目施工現(xiàn)場(見圖13)可以看出，發(fā)生地表變形的北圍地區(qū)正好被項目施工區(qū)所覆蓋，故該項目的實施可能是導(dǎo)致珠海北站區(qū)域地面沉降的主要原因。

圖12 在珠海北站及附近區(qū)域相干目標(biāo)采樣點時間序列

圖13 珠海北站TOD項目施工現(xiàn)場

4 結(jié) 論

本項目充分發(fā)揮InSAR技術(shù)在面域測量、監(jiān)測精度精度高、監(jiān)測周期短等方面的優(yōu)勢，對研究區(qū)域進(jìn)行整體監(jiān)測。通過數(shù)據(jù)處理分析，確定地表形變發(fā)生、發(fā)展和變化的全過程?；酒詹榱酥楹Ｊ卸烽T區(qū)、金灣區(qū)、香洲區(qū)在2016年1月～2017年11月期間的地表形變演化歷史。在整個監(jiān)測區(qū)域，通過數(shù)據(jù)處理分析發(fā)現(xiàn)大部分地區(qū)平均沉降速率在-17.23～5.88 mm/a之間，主城區(qū)較為穩(wěn)定。將地表形變量較大的區(qū)域標(biāo)定為災(zāi)害隱患點，隱患點多聚集在沿海區(qū)域如橫琴島、白蕉鎮(zhèn)、南水鎮(zhèn)、紅旗鎮(zhèn)等地區(qū)。針對可能發(fā)生災(zāi)害的隱患點進(jìn)行采樣分析，提取相干目標(biāo)點形變的時間序列，通過形變量時間序列分析得到，隱患點的累計下沉量最低為-38 mm，最高達(dá)到-100 mm，并且在中海石油廣東液化天然氣公司天然氣液化廠發(fā)現(xiàn)不均勻沉降現(xiàn)象。同時針對災(zāi)難隱患點收集相關(guān)資料，發(fā)現(xiàn)許多“危險區(qū)”均是受到項目施工的影響而發(fā)生地表垂直方向形變。

通過衛(wèi)星對地綜合探測和動態(tài)觀測，可有效探測研究區(qū)內(nèi)的潛在隱患點。港珠澳地區(qū)的快速城市化需要高度關(guān)注城市基礎(chǔ)設(shè)施的健康診斷，通過InSAR重復(fù)觀測可提高城市災(zāi)害風(fēng)險管控能力，可望成為未來城市基礎(chǔ)設(shè)施建設(shè)和管理標(biāo)準(zhǔn)研究的重要技術(shù)手段。

[1] 祝彥賀, 朱偉林, 徐強(qiáng), 等. 珠江口盆地中部珠海組-珠江組層序結(jié)構(gòu)及沉積特征[J]. 海洋地質(zhì)與第四紀(jì)地質(zhì), 2009 (4): 77-83.

[2] 張舉兼. 珠海市水文地質(zhì)與工程地質(zhì)特征初探[J]. 地質(zhì)災(zāi)害與環(huán)境保護(hù), 1997, 8(3): 6-12.

[3] HU B, WANG H S, SUN Y L, et al. Long-Term Land Subsidence Monitoring of Beijing (China) Using the Small Baseline Subset (SBAS) Technique [J]. Remote Sensing, 2014, 6(5):3648-3661.

[4] HU B. Monitoring of Ground Deformation due to Excessive Withdrawal of Natural Gas Using SBAS, Mathematical Problems in Engineering, vol. 2014, Article ID 674510, 2014. doi:10.1155/2014/674510.

[5] HU B, WANG Hansheng, JIANG Liming. Monitoring of the reclamation-induced ground subsidence in Macao (China) using the PSInSAR technique[J]. Journal of Central South University, 2013, Volume20, Number4, 1039-1046.

[6] 喬書波,李金嶺,孫付平,等. InSAR技術(shù)現(xiàn)狀與應(yīng)用[J]. 天文學(xué)進(jìn)展，2003(1)：11-25.

[7] MASSONNET D. Application of remote sensing data in earthquake monitoring [J]. Advances in Space Research, 1995, 15(11):37-44.

[8] 張景發(fā), 郭慶十, 龔利霞. InSAR測量滄州地下水抽取引起的地面沉降[C]//海峽兩岸空間資訊與防災(zāi)科技研討會， 珠海：2005.

[9] 許文斌, 李志偉, 丁曉利,等. 利用InSAR短基線技術(shù)估計洛杉磯地區(qū)的地表時序形變和含水層參數(shù)[J]. 地球物理學(xué)報, 2012, 55(2):452-461.

[10] Zhang Jingfa, Li Faxiang, 張景發(fā),等. The Processing Procedure of D-InSAR and Influence Factors of Its Application[J]. 地球信息科學(xué), 2000, 2(3):58-64.

[11] 劉曉萌, 常占強(qiáng), 張景發(fā),等. D-InSAR處理中失相干問題的研究[J]. 河北師范大學(xué)學(xué)報(自然科學(xué)版), 2007, 31(2):260-263.

[12] 歐陽倫曦, 李新情, 惠鳳鳴, 等. 哨兵衛(wèi)星 Sentinel-1A 數(shù)據(jù)特性及應(yīng)用潛力分析[J]. 極地研究, 2017, 29(2): 286-295.

[13] BERARDINO P, FORNARO G, LANARI R,et al. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[C]. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2375-2383.

[14] LANARI R, MORA O, MANUNTA M,et al. A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms.IEEE Trans. Geosci. Remote Sens. 2004, 42, 1377-1386.