楊亞彬，李 慧

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

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

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

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

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

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

圖1 研究區(qū)范圍

本項(xiàng)目采用32景Sentine-1A SAR影像進(jìn)行數(shù)據(jù)處理。Sentinel-1是專為“哥白尼計(jì)劃”設(shè)計(jì)的第一顆衛(wèi)星，主要用以監(jiān)測(cè)歐洲、加拿大和極地區(qū)域的環(huán)境變化，如海冰、海洋環(huán)境以及地表運(yùn)動(dòng)，地表森林、水源和土壤動(dòng)態(tài)變化。Sentinel-1衛(wèi)星是歐洲極地軌道C波段雷達(dá)成像系統(tǒng)，采用預(yù)編程無(wú)沖突操作模式，可實(shí)現(xiàn)高分辨率全球陸地、海岸帶和航線的成像[13-14]。其相對(duì)于ERS-1/1/SAR和ENVISAT/ASAR具有重訪周期短、單景影像覆蓋范圍高的優(yōu)勢(shì)，雙星座衛(wèi)星的Sentinel-1重復(fù)周期可縮短至6 d。本研究時(shí)間跨度為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) 采用日本宇宙航空開(kāi)發(fā)機(jī)構(gòu)(JAXA)聯(lián)合美國(guó)RESTEC公司利用35 km×35 km的光學(xué)影像生產(chǎn)的30 m分辨率的高精度DEM，該光學(xué)影像采集自日本發(fā)射的ALOS衛(wèi)星所搭載的PRISM(全色立體測(cè)繪儀)傳感器，由此獲得的DEM平面及高程精度均可達(dá)到5 m，其數(shù)據(jù)產(chǎn)品能在全球范圍內(nèi)覆蓋。相較于SRTM數(shù)據(jù)，AW3D30的數(shù)字高程產(chǎn)品由于不存在雷達(dá)信號(hào)失相干的問(wèn)題，因此不會(huì)出現(xiàn)數(shù)據(jù)空洞的問(wèn)題。利用該DEM進(jìn)行去除地形處理，能夠精確提供研究區(qū)域的高程數(shù)據(jù)，并且能有效降低由于外部提供的DEM不正確而產(chǎn)生的干涉去平地效應(yīng)出錯(cuò)等情況，使得InSAR結(jié)果的可靠性得以提高。

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

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

2.1 差分干涉圖的生成

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

圖2 SAR干涉對(duì)連接圖

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

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

2.3 高相干點(diǎn)沉降速率和形變場(chǎng)的獲取

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

3 結(jié)果分析

本項(xiàng)目通過(guò)對(duì)32景Sentinel-1A衛(wèi)星IW模式的SAR影像進(jìn)行數(shù)據(jù)處理，得到研究區(qū)的年平均沉降速率圖如圖3所示。從圖3中可看出該監(jiān)測(cè)區(qū)沉降速率最大可達(dá)-101.05 mm/a，抬升速率最大達(dá)到75.25 mm/a。監(jiān)測(cè)區(qū)南部沿海區(qū)域存在輕微抬升現(xiàn)象，抬升速率處于17.44 mm/a以內(nèi)。珠海市斗門區(qū)大部分地區(qū)處于穩(wěn)定狀態(tài)，年平均沉降速率不超過(guò)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)較為明顯的沉降信號(hào)，沉降速率最大達(dá)到-75 mm/a。香洲區(qū)大部分地區(qū)處于年平均形變速率在-17～17 mm/a之間，較為穩(wěn)定，但橫琴等區(qū)域存在明顯的沉降漏斗。澳門大部分區(qū)域的年平均形變速率處于5.88～17.44mm/a之間，少數(shù)地區(qū)存在沉降信號(hào)，最大沉降量達(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ū)域。對(duì)此區(qū)域選定采樣點(diǎn)#1、#2、#3、#4進(jìn)行形變量時(shí)間序列提取得到圖5，從圖中可看出，采樣點(diǎn)均呈現(xiàn)一致的下沉趨勢(shì)，且累計(jì)沉降量可以達(dá)到將近-132 mm，沉降趨勢(shì)明顯。

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

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

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

由平均沉降速率圖3可以看出，珠海市金灣區(qū)主城區(qū)處于穩(wěn)定狀態(tài)，年平均形變速率不超過(guò)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處地方存在較為明顯的沉降信號(hào)。在沉降信號(hào)明顯區(qū)域抽取有代表性的相干目標(biāo)點(diǎn)，提取其垂直形變的時(shí)間序列進(jìn)行分析。

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

通過(guò)衛(wèi)星對(duì)地綜合探測(cè)和動(dòng)態(tài)觀測(cè)，可有效探測(cè)研究區(qū)內(nèi)的潛在隱患點(diǎn)。港珠澳地區(qū)的快速城市化需要高度關(guān)注城市基礎(chǔ)設(shè)施的健康診斷，通過(guò)InSAR重復(fù)觀測(cè)可提高城市災(zāi)害風(fēng)險(xiǎn)管控能力，可望成為未來(lái)城市基礎(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è)量滄州地下水抽取引起的地面沉降[C]//海峽兩岸空間資訊與防災(zāi)科技研討會(huì)， 珠海：2005.

[9] 許文斌, 李志偉, 丁曉利,等. 利用InSAR短基線技術(shù)估計(jì)洛杉磯地區(qū)的地表時(shí)序形變和含水層參數(shù)[J]. 地球物理學(xué)報(bào), 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處理中失相干問(wèn)題的研究[J]. 河北師范大學(xué)學(xué)報(bào)(自然科學(xué)版), 2007, 31(2):260-263.

[12] 歐陽(yáng)倫曦, 李新情, 惠鳳鳴, 等. 哨兵衛(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.