楊少敏，范雪婷，劉 波，李紅慧

（1.常州市基礎(chǔ)地理勘測(cè)中心，江蘇常州 213001；2.江蘇省基礎(chǔ)地理信息中心，江蘇南京 210013；3.江蘇煤炭地質(zhì)局，江蘇常州 213000）

0 引言

長(zhǎng)江三角洲蘇錫常地區(qū)因地下水開采導(dǎo)致的地面沉降屬于一種緩慢地質(zhì)災(zāi)害現(xiàn)象，可能會(huì)對(duì)當(dāng)?shù)厣鐣?huì)生活、經(jīng)濟(jì)活動(dòng)及生態(tài)環(huán)境帶來(lái)不利影響，嚴(yán)重時(shí)甚至?xí)觿〕鞘泻闈碁?zāi)害，威脅生命安全。常州市地處長(zhǎng)江三角洲經(jīng)濟(jì)帶的核心區(qū)，二十世紀(jì)六七十年代已出現(xiàn)較為明顯沉降現(xiàn)象，八十年代沉降主要集中位于城區(qū)范圍［1］，隨著經(jīng)濟(jì)快速發(fā)展，其范圍也逐漸外擴(kuò)偏移，出現(xiàn)了大面積不均勻沉降中心，因此地面沉降監(jiān)測(cè)得到政府和相關(guān)研究者持續(xù)廣泛的關(guān)注。如于軍等［2］采用Envisat衛(wèi)星ASAR數(shù)據(jù)對(duì)常州地區(qū)2003—2008年的地面沉降進(jìn)行監(jiān)測(cè)；劉波等［3］基于SBAS技術(shù)對(duì)武進(jìn)地區(qū)2007—2011年地面進(jìn)行沉降監(jiān)測(cè)；董少春等［4］采用Sentinel-1A影像對(duì)常州地區(qū)2015—2018年地表沉降進(jìn)行監(jiān)測(cè)。以上研究都能夠監(jiān)測(cè)到該地區(qū)地表沉降時(shí)空分布特征，但仍存在不足：一是監(jiān)測(cè)周期短；二是僅監(jiān)測(cè)獲取沉降信息，對(duì)于監(jiān)測(cè)結(jié)果準(zhǔn)確性及因沉降導(dǎo)致的危險(xiǎn)性分析涉及較少。本文利用2007—2018年影像數(shù)據(jù)，基于多主影像相干目標(biāo)小基線技術(shù)［5］（multiple images coherent targets small baselines InSAR，MCTSB-InSAR）進(jìn)行時(shí)序地表沉降監(jiān)測(cè)，利用傳統(tǒng)水準(zhǔn)測(cè)量數(shù)據(jù)進(jìn)行精度分析評(píng)定，并最終對(duì)常州市地面沉降進(jìn)行危險(xiǎn)性評(píng)價(jià)，為常州地區(qū)地面沉降地質(zhì)災(zāi)害防治與預(yù)警預(yù)報(bào)提供依靠數(shù)據(jù)。

1 數(shù)據(jù)

本次研究分別選取ALOSPALSAR、Radarsat-2、Sentinel-1 A作為常州市2007—2011年、2012—2015年、2016—2018年三階段沉降監(jiān)測(cè)的數(shù)據(jù)源，數(shù)據(jù)皆為單視復(fù)數(shù)據(jù)（SLC），不同SAR數(shù)據(jù)的覆蓋情況如圖1所示。PALSAR影像工作頻率為L(zhǎng)波段，波束模式分為精細(xì)波束雙極化（Fine Beam Dualpolarization，F(xiàn)BD）模式及精細(xì)波束單極化（Fine Beam Single Polarization，F(xiàn)BS）模式，二者具有不同的方位向分辨率，時(shí)序InSAR處理時(shí)將FBS數(shù)據(jù)在距離上進(jìn)行2倍降采樣，使其與FBD數(shù)據(jù)有相同分辨率。單景影像覆蓋范圍約60 km×60 km，共需6個(gè)Frame覆蓋常州全市，每個(gè)Frame獲取16～21期，共114景影像。Radarsat-2影像工作頻率為C波段，寬幅模式數(shù)據(jù)單景覆蓋范圍約150 km×150 km，覆蓋常州需2個(gè)Frame，各獲取25期共計(jì)50景影像。Sentinel-1 A影像工作頻率為C波段，單景影像覆蓋范圍約250 km×180 km，僅1個(gè)Frame即可覆蓋常州全市，獲取47景（期）影像。SAR數(shù)據(jù)詳細(xì)參數(shù)見表1所示。研究中采用分辨率為90 m的SRTMDEM數(shù)據(jù)來(lái)剔除地形相位的影響。

圖1 常州市SAR影像覆蓋范圍

表1 SAR影像參數(shù)

2 方法

2.1 時(shí)間序列InSAR處理

采用MCTSB-InSAR方法對(duì)研究數(shù)據(jù)進(jìn)行時(shí)間序列InSAR處理，具體處理流程如圖2所示，主要處理步驟如下。

圖2 MCTSB-InSAR數(shù)據(jù)處理流程

2.1.1 差分干涉處理

主要涵蓋軌道加密、基線分析、影像配準(zhǔn)、組合干涉像對(duì)、獲取干涉圖、去地形相位、相位濾波等處理，最終生成差分干涉圖。其中，因Sentinel-1A數(shù)據(jù)采用新型的TOPS（Terrain Observation with Progressive Scans SAR）模式，InSAR處理在方位向上有很大的多普勒中心頻率偏移，對(duì)配準(zhǔn)提出了更高的要求［6］，本文采用基于精軌數(shù)據(jù)、強(qiáng)度互相關(guān)、增強(qiáng)譜分集的多級(jí)配準(zhǔn)方法，避免了干涉圖中出現(xiàn)相位跳躍現(xiàn)象，得到優(yōu)于千分之一像素的方位向配準(zhǔn)精度。

2.1.2 線性形變反演

依次進(jìn)行穩(wěn)定點(diǎn)目標(biāo)提取、點(diǎn)目標(biāo)三角網(wǎng)連接、建立差分相位模型、線性形變速率和高程誤差求解等處理。精確的穩(wěn)定點(diǎn)目標(biāo)識(shí)別是高精度時(shí)序InSAR地面沉降信息提取的基礎(chǔ)，穩(wěn)定點(diǎn)目標(biāo)指對(duì)雷達(dá)波的后向散射較強(qiáng)，并且在時(shí)序上較穩(wěn)定的各種地物目標(biāo)，如建筑物與構(gòu)筑物、鐵軌、橋梁、欄桿、路燈等人工目標(biāo)及裸露的巖石等自然目標(biāo)。本文綜合使用基于相干系數(shù)、幅度離差、平均幅度的“三閾值二級(jí)探測(cè)法”［7］提取穩(wěn)定點(diǎn)目標(biāo)，分布示例如圖3所示。然后，利用局部Delaunay三角網(wǎng)連接點(diǎn)目標(biāo)，三角網(wǎng)的參數(shù)為格網(wǎng)中心間隔600 m，格網(wǎng)半徑1 000 m，相鄰點(diǎn)目標(biāo)間的最大距離為2 000 m。借助三角網(wǎng)確定點(diǎn)目標(biāo)之間關(guān)系，對(duì)相鄰點(diǎn)目標(biāo)的差分相位實(shí)施二次差分，構(gòu)建點(diǎn)目標(biāo)鄰域差分相位模型，并根據(jù)文獻(xiàn)［5］進(jìn)行模型參數(shù)反演，最終求取各點(diǎn)目標(biāo)的線性形變速率和高程誤差。

圖3 穩(wěn)定點(diǎn)目標(biāo)分布

2.1.3 非線性形變反演

從原始差分干涉圖相位減去高程誤差和線性形變速率分量，可獲得殘余相位，其主要包含非線性形變相位、大氣相位及噪聲相位。根據(jù)各分量不同的時(shí)空頻譜特征，通過(guò)時(shí)間域低通濾波分離出非線性形變相位；對(duì)低通濾波后的殘余相位進(jìn)行空間域低通濾波，分離得到大氣相位。將提取的線性形變和非線性形變進(jìn)行疊加計(jì)算得出需要的全部形變信息，并基于雷達(dá)波入射角關(guān)系，將雷達(dá)視線向的形變值轉(zhuǎn)換為垂直方向上的形變值。

2.1.4 多景InSAR地面沉降結(jié)果拼接

對(duì)于ALOSPAlSAR和RADARSAT-2數(shù)據(jù)，需多個(gè)Frame影像覆蓋研究區(qū)，在分別獲取每個(gè)Frame監(jiān)測(cè)點(diǎn)結(jié)果后，通過(guò)提取影像間重疊區(qū)域的同名點(diǎn)結(jié)果，采用最小二乘平差方法對(duì)多個(gè)Frame點(diǎn)目標(biāo)結(jié)果進(jìn)行拼接，最終形成區(qū)域協(xié)調(diào)一致的沉降結(jié)果［8］。

2.2 精度評(píng)定方法

對(duì)于InSAR技術(shù)獲取的地表沉降結(jié)果，通常采用高等級(jí)水準(zhǔn)測(cè)量或GPS/CORS測(cè)量數(shù)據(jù)進(jìn)行定點(diǎn)比較的方法評(píng)定其精度。因InSAR監(jiān)測(cè)點(diǎn)與水準(zhǔn)點(diǎn)的空間位置不一定完全一致，為了保證評(píng)定結(jié)果的客觀和準(zhǔn)確，約定鄰近點(diǎn)原則對(duì)初始獲得的水準(zhǔn)點(diǎn)進(jìn)行了篩選：以水準(zhǔn)點(diǎn)為中心，如果在一定距離范圍內(nèi)（經(jīng)驗(yàn)值為80 m）存在至少一個(gè)InSAR點(diǎn)目標(biāo)，則選取該水準(zhǔn)點(diǎn)參與精度評(píng)定，比較水準(zhǔn)點(diǎn)與最鄰近InSAR點(diǎn)目標(biāo)的年平均沉降速率，計(jì)算二者的差值的均方根誤差，作為InSAR地面沉降監(jiān)測(cè)精度評(píng)定的指標(biāo)。

2.3 危險(xiǎn)性評(píng)價(jià)方法

地面沉降危險(xiǎn)性是指在不同外界因素作用下地面發(fā)生沉降的概率及程度［9］。通常研究地面沉降危險(xiǎn)性主要采用累計(jì)沉降量和平均沉降速率兩種評(píng)價(jià)指標(biāo)［10］，前者反映地面總沉降量值，后者反映該區(qū)域的沉降速率，這兩個(gè)指標(biāo)都可呈現(xiàn)研究區(qū)地面沉降的發(fā)展程度，在日常地質(zhì)災(zāi)害危險(xiǎn)性評(píng)估工作中最易被使用。

根據(jù)地面沉降危險(xiǎn)性評(píng)價(jià)指標(biāo)，綜合衡量累計(jì)地面沉降量與地面沉降速率間關(guān)系及其對(duì)地面沉降危險(xiǎn)性的影響后，確定彼此間權(quán)重各占0.5［11］。結(jié)合《江蘇省地質(zhì)災(zāi)害危險(xiǎn)性評(píng)估技術(shù)要求》和《江蘇省地面沉降控制區(qū)劃分方案》以及數(shù)據(jù)的實(shí)際獲取情況，基于監(jiān)測(cè)期平均地面沉降速率和地面累計(jì)沉降量數(shù)據(jù)對(duì)研究區(qū)地面沉降進(jìn)行評(píng)分分級(jí)，分級(jí)結(jié)果如表2所示。

表2 地表沉降危險(xiǎn)性評(píng)價(jià)指標(biāo)參數(shù)及賦值

表3 地表沉降危險(xiǎn)性分區(qū)標(biāo)準(zhǔn)

利用綜合指數(shù)法對(duì)兩個(gè)評(píng)價(jià)指標(biāo)按權(quán)重計(jì)算后疊加得到危險(xiǎn)值，將最終危險(xiǎn)性參照評(píng)分分級(jí)表進(jìn)行分區(qū)（見表3），獲取研究區(qū)地面沉降危險(xiǎn)性評(píng)估分區(qū)結(jié)果。

3 結(jié)果分析與評(píng)價(jià)

3.1 結(jié)果分析

基于上述方法，分別提取了常州市2007—2011年、2012—2015年、2016—2018年的地表沉降速率（見圖4—6）和累計(jì)沉降量等信息。采用克里金（Kriging）插值計(jì)算方法將穩(wěn)定點(diǎn)上的沉降速率值和累計(jì)沉降量值進(jìn)行內(nèi)插，將三期差值后的沉降量進(jìn)行疊加計(jì)算得到2007—2018年的總累計(jì)沉降量圖（見圖7）。

2007—2011年，常州市地表形變相對(duì)較為嚴(yán)重，插值后沉降速率超過(guò)10 mm/a的區(qū)域面積達(dá)84.8 km2。其中鐘樓區(qū)鄒區(qū)鎮(zhèn)存在一處明顯的沉降漏斗，最大沉降速率達(dá)57.7 mm/a；武進(jìn)區(qū)存在一條縱貫?zāi)媳钡某两祹?，涉及牛塘?zhèn)、湖塘鎮(zhèn)、高新區(qū)、禮嘉鎮(zhèn)、洛陽(yáng)鎮(zhèn)、前黃鎮(zhèn)、雪堰鎮(zhèn)、橫山橋鎮(zhèn)，最大沉降速率達(dá)56.5 mm/a，位于高新區(qū)；金壇區(qū)和溧陽(yáng)市個(gè)別地區(qū)出現(xiàn)輕微的不均勻沉降，最大沉降速率分別為22.2和42.2 mm/a，分別位于西城街道和別橋鎮(zhèn)；天寧區(qū)和新北區(qū)未發(fā)現(xiàn)明顯的沉降中心。

2012—2015年，常州市地面沉降速率明顯趨緩，超過(guò)10 mm/a的區(qū)域面積減少為5.6 km2。鐘樓區(qū)，金壇區(qū)與溧陽(yáng)市已無(wú)明顯沉降跡象；武進(jìn)區(qū)的沉降帶依然存在，但范圍減小，主要涉及牛塘鎮(zhèn)、湖塘鎮(zhèn)、高新區(qū)、禮嘉鎮(zhèn)、洛陽(yáng)鎮(zhèn)、雪堰鎮(zhèn)，最大沉降速率降為21.8 mm/a，位于高新區(qū)；天寧區(qū)和新北區(qū)仍較為穩(wěn)定，無(wú)明顯沉降。

圖4 2007—2011年平均沉降速率

圖6 2016—2018年平均沉降速率

圖5 2012—2015年平均沉降速率

圖7 2007—2018年累計(jì)沉降量

2016—2018年，常州市地面沉降范圍擴(kuò)大，沉降速率超過(guò)10 mm/a的區(qū)域面積增加至32.6 km2。鐘樓區(qū)鄒區(qū)鎮(zhèn)出現(xiàn)沉降反彈，又有輕微沉降發(fā)生，最大沉降速率為15 mm/a；武進(jìn)區(qū)禮嘉鎮(zhèn)、高新區(qū)、前黃鎮(zhèn)、洛陽(yáng)鎮(zhèn)沉降已連接成片，最大沉降速率仍位于高新區(qū)，達(dá)27.8 mm/a。金壇區(qū)、溧陽(yáng)市有零星的沉降出現(xiàn)。

依據(jù)2007—2018年累計(jì)沉降量，統(tǒng)計(jì)得出沉降量超過(guò)100 mm的區(qū)域面積達(dá)60.2 km2，超過(guò)200 mm的區(qū)域面積達(dá)12.5 km2，可見對(duì)整個(gè)階段的沉降發(fā)展趨勢(shì)需時(shí)刻保持高度警惕。

3.2 精度評(píng)定

本次共收集到常州市47個(gè)水點(diǎn)準(zhǔn)數(shù)據(jù)，2012—2015年每年觀測(cè)一期，可用于RADARSAT-2監(jiān)測(cè)結(jié)果的驗(yàn)證。根據(jù)上述精度評(píng)定方法，共篩選出33個(gè)水準(zhǔn)點(diǎn)（位置分布如圖8所示），計(jì)算得到的沉降速率的均方差為2.7 mm/a，滿足規(guī)范（CH/T 6006—2018）中精度均優(yōu)于5 mm/a的要求（見表4）。

選取位于武進(jìn)區(qū)禮嘉鎮(zhèn)的水準(zhǔn)點(diǎn)ZY076，利用2012—2015年的時(shí)序?qū)崪y(cè)數(shù)據(jù)與InSAR監(jiān)測(cè)點(diǎn)的時(shí)序沉降量進(jìn)行計(jì)算分析（見圖9），可見二者也具有較高的符合性，進(jìn)一步驗(yàn)證了采用MCTSB-InSAR技術(shù)獲取地表形變的可靠性。

3.3 危險(xiǎn)性評(píng)價(jià)

圖8 水準(zhǔn)點(diǎn)分布位置

本文基于2016—2018年平均地面沉降速率和地面累計(jì)沉降量數(shù)據(jù)結(jié)合評(píng)分分級(jí)表對(duì)常州市地面沉降進(jìn)行評(píng)分分級(jí)；利用綜合指數(shù)法對(duì)2個(gè)評(píng)價(jià)指標(biāo)按權(quán)重計(jì)算后疊加得到危險(xiǎn)值，將最終危險(xiǎn)性參照危險(xiǎn)性分區(qū)標(biāo)準(zhǔn)進(jìn)行分區(qū)，得到常州市地面沉降危險(xiǎn)性評(píng)估分區(qū)圖（見圖10）。

易發(fā)現(xiàn)，常州市無(wú)明顯沉降高危險(xiǎn)區(qū)；中危險(xiǎn)區(qū)面積達(dá)15.85 km2，主要分布于鐘樓區(qū)的鄒區(qū)鎮(zhèn)、武進(jìn)區(qū)的牛塘鎮(zhèn)、高新區(qū)、禮嘉鎮(zhèn)、洛陽(yáng)鎮(zhèn)及雪堰鎮(zhèn)，總體呈西北—東南走向，此外金壇區(qū)的直溪鎮(zhèn)西部、金城鎮(zhèn)與東城街道交界附近也出現(xiàn)中危險(xiǎn)區(qū)；剩余地區(qū)地表基本處于穩(wěn)定狀態(tài)為低危險(xiǎn)區(qū)?；诜謪^(qū)結(jié)果結(jié)合調(diào)查分析得知，常州市地面沉降主要因新增大量建筑工程、新區(qū)開發(fā)利用強(qiáng)度加大導(dǎo)致地面下沉（如鄒區(qū)鎮(zhèn)燈博廣場(chǎng)、高新區(qū)全面開啟特色園區(qū)建設(shè)等）。因常州市全面禁止開采地下水，地下水位得到明顯回升，相較以往地表沉降速率已得到有效減緩甚至出現(xiàn)輕微抬升反彈，但因地下水開采導(dǎo)致含水層中出現(xiàn)的“失衡層”，其恢復(fù)是個(gè)蠕變過(guò)程，在沉降上呈現(xiàn)滯后效應(yīng)。

4 結(jié)語(yǔ)

本文利用114景ALOSPALSAR影像、50景Radarsat-2影像和47景Sentinel-1 A影像，基于MCTSB-InSAR方法對(duì)常州市進(jìn)行沉降監(jiān)測(cè)，獲取了2007—2011、2012—2015、2016—2018三階段年平均沉降速率和2007—2018年累計(jì)沉降量信息。通過(guò)與水準(zhǔn)測(cè)量結(jié)果比較，常州市地面沉降速率的精度評(píng)定達(dá)2.7 mm/a，驗(yàn)證了監(jiān)測(cè)結(jié)果的客觀性和可靠性?；诔两当O(jiān)測(cè)結(jié)果系統(tǒng)分析了常州市沉降漏斗分布特征和時(shí)空變化規(guī)律，地表形變總體呈現(xiàn)先明顯減緩后輕微反彈的趨勢(shì)。為了全面分析地質(zhì)災(zāi)害危險(xiǎn)性，利用綜合指數(shù)法分權(quán)疊加分析，實(shí)現(xiàn)常州市地面沉降災(zāi)害危險(xiǎn)性評(píng)估分區(qū)，結(jié)果表明常州市無(wú)地質(zhì)災(zāi)害高危險(xiǎn)區(qū)，主要在鐘樓和武進(jìn)區(qū)存在西北東南走向中危險(xiǎn)區(qū)，分區(qū)評(píng)價(jià)結(jié)果也可為今后地面沉降災(zāi)害防治工作提供基礎(chǔ)數(shù)據(jù)和支撐依據(jù)。

表4 水準(zhǔn)測(cè)量結(jié)果與InSAR沉降速率對(duì)比 （單位：mm/a）

圖9 ZY076水準(zhǔn)點(diǎn)與InSAR沉降量對(duì)比

圖10 地表沉降危險(xiǎn)性分區(qū)