李克沖,董張玉,楊學(xué)志
1(合肥工業(yè)大學(xué) 計(jì)算機(jī)與信息學(xué)院,合肥 230601)
2(工業(yè)安全與應(yīng)急技術(shù)安徽省重點(diǎn)實(shí)驗(yàn)室,合肥 230601)
3(合肥工業(yè)大學(xué) 軟件學(xué)院,合肥 230601)
地面沉降是在自然變化和人類(lèi)活動(dòng)綜合影響下,引起地殼表面標(biāo)高緩慢變化的一種局部工程地質(zhì)現(xiàn)象.隨著全球城市化進(jìn)展不斷加快,城市擴(kuò)張和人口快速增長(zhǎng),城市基礎(chǔ)設(shè)施的大量修建和地下資源的過(guò)度開(kāi)采,導(dǎo)致城市地面沉降日益嚴(yán)重.地面沉降可對(duì)建筑物和生產(chǎn)設(shè)施造成嚴(yán)重破壞,不利于城市建設(shè)和資源勘測(cè)開(kāi)發(fā),甚至可能在沿海地區(qū)造成海水倒灌,對(duì)人民生命財(cái)產(chǎn)安全和基礎(chǔ)設(shè)施安全造成了重大安全隱患[1,2].傳統(tǒng)的調(diào)查和偵查方法如全球定位系統(tǒng)(global positioning system,GPS)和水準(zhǔn)測(cè)量等監(jiān)測(cè)方法雖然精度較高,但是只能在小范圍的點(diǎn)和線進(jìn)行測(cè)量,難以對(duì)大范圍區(qū)域進(jìn)行監(jiān)測(cè),且對(duì)監(jiān)測(cè)點(diǎn)的長(zhǎng)期維護(hù)也比較困難[3–6].因此,合成孔徑雷達(dá)干涉測(cè)量(interferometric synthetic aperture radar,InSAR)憑借著不受時(shí)間氣候環(huán)境影響和高精度、高分辨率等優(yōu)點(diǎn)在地面沉降監(jiān)測(cè)中顯示出極大的應(yīng)用潛力[7–9].
InSAR 技術(shù)最早被應(yīng)用于地形測(cè)繪領(lǐng)域,文獻(xiàn)[10]首次提出合成孔徑雷達(dá)差分干涉測(cè)量(differential interferometric synthetic aperture radar,D-InSAR)技術(shù),研究表明D-InSAR 監(jiān)測(cè)地表形變的精度可以達(dá)到厘米級(jí),自此D-InSAR 技術(shù)開(kāi)始被應(yīng)用于地表形變監(jiān)測(cè)中;文獻(xiàn)[11]利用D-InSAR 技術(shù)監(jiān)測(cè)伊朗達(dá)馬內(nèi)市的地面沉降;文獻(xiàn)[12] 利用永久散射體技術(shù)(persistent scatterer interferometric synthetic aperture radar,PS-InSAR)對(duì)意大利的托斯卡納滑坡地區(qū)進(jìn)行形變監(jiān)測(cè)研究;文獻(xiàn)[13]利用短基線集技術(shù)對(duì)克羅地亞的薩格勒布地區(qū)進(jìn)行地面形變分析.
我國(guó)將InSAR 技術(shù)應(yīng)用于地面形變監(jiān)測(cè)較晚,但也取得了許多成果.文獻(xiàn)[14]利用PS-InSAR 技術(shù)對(duì)上海高架路進(jìn)行沉降監(jiān)測(cè),分析提取的高架路沉降的空間分布以及時(shí)間變化情況;2019年文獻(xiàn)[15]采用時(shí)序InSAR 技術(shù)監(jiān)測(cè)三峽庫(kù)區(qū)藕塘滑坡的穩(wěn)定性;文獻(xiàn)[16]利用Sentinel-1A 影像分別采用PS-InSAR 技術(shù)和SBASInSAR 技術(shù)監(jiān)測(cè)南京市地面沉降.D-InSAR 技術(shù)可以達(dá)到厘米級(jí)的測(cè)量精度,但是難以對(duì)同一區(qū)域進(jìn)行時(shí)間序列形變監(jiān)測(cè);PS-InSAR和SBAS-InSAR 都可利用大量合成孔徑雷達(dá)(synthetic aperture radar,SAR)影像(一般多于20 幅)進(jìn)行時(shí)間序列分析,但是在研究SAR影像較少的區(qū)域時(shí),上述兩種方法的應(yīng)用均受到限制.
傳統(tǒng)的干涉圖疊加方法(interferogram stacking)[17–19]是將地面形變量視為一種線性變化,將相互獨(dú)立的干涉圖中所包含大氣擾動(dòng)的相位誤差在時(shí)間上近似為隨機(jī)量,將研究區(qū)域的地表形變信息視為線性變化,通過(guò)線性疊加,提高了疊加相位圖中形變信息與大氣干擾噪聲之間的信噪比.該方法在僅有少量的SAR 數(shù)據(jù)的支撐下也可實(shí)現(xiàn)較長(zhǎng)時(shí)間的形變監(jiān)測(cè).但是在沒(méi)有考慮干涉圖質(zhì)量權(quán)重的情況下,使相干性較低的干涉圖參與疊加,勢(shì)必會(huì)影響監(jiān)測(cè)結(jié)果的精度.
綜上所述,本文提出一種面向地面沉降的基于多幅主影像干涉圖加權(quán)疊加的方法.該方法通過(guò)選取多幅主影像,并在傳統(tǒng)Stacking 法上在加入圖像質(zhì)量權(quán)重的方法對(duì)地表進(jìn)行形變監(jiān)測(cè),相對(duì)傳統(tǒng)的監(jiān)測(cè)方法提高監(jiān)測(cè)精度.
本文的其余章節(jié)將簡(jiǎn)要安排如下:第1 節(jié)介紹研究區(qū)域以及數(shù)據(jù).第2 節(jié)介紹多幅主影像干涉圖加權(quán)疊加方法的基本原理,并就所提方法的主要思想、流程及計(jì)算公式進(jìn)行詳細(xì)介紹.第3 節(jié)就監(jiān)測(cè)結(jié)果從主觀和客觀兩個(gè)方面進(jìn)行評(píng)價(jià)并給出分析,總結(jié)部分將在第4 節(jié)給出.
美國(guó)的加利福尼亞州一直以來(lái)存在著不同程度的沉降.在2014年加利福尼亞州出現(xiàn)了近50年來(lái)最嚴(yán)重的地面沉降,后續(xù)預(yù)計(jì)地面沉降將繼續(xù)惡化.
選取圣迭戈縣作為實(shí)驗(yàn)區(qū)域,圣迭戈縣的地形較為復(fù)雜,其西側(cè)是蜿蜒的海岸線,東北側(cè)主要是大雪覆蓋的雪峰,而東部則為大片的金黃色沙漠.因此這次實(shí)驗(yàn)主要選取城區(qū)監(jiān)測(cè)沉降,避開(kāi)山脈,沙漠等相干性低的區(qū)域.
選用的實(shí)驗(yàn)數(shù)據(jù)包括哨兵1A (Sentinel-1A) 數(shù)據(jù)、數(shù)字高程模型(digital elevation model,DEM)、以及地下水監(jiān)測(cè)數(shù)據(jù),這一數(shù)據(jù)來(lái)自美國(guó)地質(zhì)調(diào)查局(United States geological survey,USGS).
Sentinel-1A 衛(wèi)星是歐洲航天局于2014年發(fā)射的首顆地球觀測(cè)衛(wèi)星,載有C 波段合成孔徑雷達(dá),可以在各種環(huán)境條件下提供影像.對(duì)同一區(qū)域的重訪周期為12 天.Sentinel-1A 數(shù)據(jù)基本參數(shù)見(jiàn)表1所列.本文使用的數(shù)字高程數(shù)據(jù)為SRTM-1 DEM 數(shù)據(jù),平面分辨率為30 m,每個(gè)DEM 子單元的數(shù)據(jù)幅寬為經(jīng)緯度1°.
表1 Sentinel-1A 衛(wèi)星4 種工作模式基本參數(shù)
考慮到多幅主影像干涉圖加權(quán)疊加方法監(jiān)測(cè)形變的優(yōu)越性,選擇圣迭戈縣城區(qū)作為試驗(yàn)區(qū),避免相干性低的區(qū)域,以提高監(jiān)測(cè)精度.采用2019年7月至2019年12月Sentinel-1A的9 景IW 模式的SLC (single looking complex)數(shù)據(jù)進(jìn)行試驗(yàn).
首先從這9 幅SAR 影像中擇優(yōu)選取3 幅作為主影像,使其余影像與其配準(zhǔn)并生成干涉對(duì),可利用綜合函數(shù)模型來(lái)進(jìn)行確定.
其中,dl是以第l幅影像為主影像時(shí)正則化距離的平均值,選取計(jì)算出的最小值所對(duì)應(yīng)的影像作為公共主影像[20];l代表SAR 影像序號(hào)(按時(shí)間前后排序);Bc、Tc和fc分別為垂直空間基線臨界值、時(shí)間基線臨界值和多普勒質(zhì)心頻率差臨界值;Bi,l、Ti,l和fi,l分別為第l幅影像作為主影像與第i幅輔影像之間的垂直空間基線、時(shí)間基線和多普勒質(zhì)心頻率差.根據(jù)上式,計(jì)算所有影像的平均正則化距離d,選取平均值最小的3 幅影像作為主影像,最終確定軌道號(hào)為027961、028836、029711的3 幅影像作為主影像.
基于D-InSAR的具體工作流程,首先選取027961、028836、029711 共3 幅影像作為干涉主影像,選取027961 作為參考影像.并將所有的圖像配準(zhǔn)到參考影像上;接著使用精密衛(wèi)星軌道數(shù)據(jù)消除衛(wèi)星軌道殘余相位誤差;選擇所有空間垂直基線中小于50 m的11 個(gè)干涉對(duì)進(jìn)行干涉處理,以抑制空間去相關(guān)引起的相位誤差,生成的干涉像對(duì)見(jiàn)表2所列.主要過(guò)程包括相干性計(jì)算、去平、濾波和相位解纏.其中濾波過(guò)程采用Goldstein 方法對(duì)干涉圖進(jìn)行濾波,這種濾波方法可以提高干涉條紋的清晰度、減少了由空間基線或時(shí)間基線引起的失相干噪聲.再將配準(zhǔn)的干涉圖與STRM1 DEM 數(shù)據(jù)進(jìn)行差分干涉,生成差分干涉圖集,選取部分如圖1所示.隨后的相位解纏過(guò)程使用最小費(fèi)用流(minimum cost flow,MCF)方法[21]進(jìn)行,弱化低相干區(qū)域?qū)?shí)驗(yàn)結(jié)果的影響,最終得到各解纏相位.
圖1 差分干涉圖集
表2 Sentinel-1A 影像和其構(gòu)成的干涉像對(duì)
多幅主影像干涉圖加權(quán)疊加方法是選擇較優(yōu)的多幅主影像,剩余的影像與主影像進(jìn)行配準(zhǔn),對(duì)配準(zhǔn)到同一主影像生成的干涉像對(duì)篩選出時(shí)空基線較短的干涉對(duì)生成不同的集合,各個(gè)集合之間的基線較長(zhǎng),通過(guò)最小二乘法得到每個(gè)集合的地面形變信息,再使用奇異分解法將各個(gè)集合聯(lián)合后求解.這種方法可以在數(shù)據(jù)集較少的基礎(chǔ)上獲得較高的監(jiān)測(cè)精度.
多幅主影像干涉圖加權(quán)疊加方法數(shù)據(jù)的處理流程如圖2所示,主要包括:數(shù)據(jù)預(yù)處理、選取多幅主影像;基線估算生成連接圖,參考DEM和軌道數(shù)據(jù)生成干涉圖;干涉圖濾波,并計(jì)算相干性;相位解纏;剔除相干性差的像對(duì),干涉圖加權(quán)疊加;相位轉(zhuǎn)形變,地理編碼;生成形變圖.
圖2 多幅主影像干涉圖加權(quán)疊加方法流程圖
在地表形變速率計(jì)算過(guò)程中,對(duì)于任一幅經(jīng)過(guò)相位解纏的干涉圖,其干涉相位可表示為:
其中,φtopo為地形相位,φdefo為兩次 獲取SAR 影像時(shí)地表在視線方向上形變引起的相位變化,φobj是由目標(biāo)散射特性變化導(dǎo)致的相位變化,φorb為軌道誤差引起的相位變化,αatmos為大氣擾動(dòng)引起的相位變化,nnoise是由熱噪聲引起的相位變化.
在D-InSAR 處理過(guò)程中,nnoise一般忽略;多普勒質(zhì)心頻率在使用同種數(shù)據(jù)源時(shí)的變化很小,φobj可以忽略不計(jì);根據(jù)DEM 數(shù)據(jù)和各原始影像對(duì)應(yīng)的精密軌道數(shù)據(jù)計(jì)算得到 φorb和φtopo.
在干涉圖疊加過(guò)程中,會(huì)將所有經(jīng)過(guò)解纏后的干涉圖進(jìn)行相位疊加,但這樣也會(huì)讓相干性差的干涉圖參與相位疊加過(guò)程,導(dǎo)致解算出的形變量產(chǎn)生較大誤差,降低監(jiān)測(cè)精度.因此,需要對(duì)干涉圖疊加時(shí)各圖像的權(quán)重進(jìn)行合理的設(shè)定,以剔除相干性差的干涉圖.首先需要基于干涉處理過(guò)程中生成的相干圖來(lái)統(tǒng)計(jì)相干點(diǎn)的數(shù)量,然后計(jì)算相干圖中各相干點(diǎn)的相干值,根據(jù)平均相干值來(lái)設(shè)置權(quán)重.將平均相干值最大的相干圖權(quán)重wi,jk設(shè)置為1,當(dāng)各相干圖中平均相干值的大小明顯低于平均大小時(shí)(經(jīng)試驗(yàn)證明,當(dāng)相干圖中平均相干值的大小低于最高相干圖平均相干值大小的3/4 時(shí)),需要將此相干圖疊加權(quán)重wi,jk設(shè)置為0,即不參與干涉圖疊加過(guò)程,否則將此相干圖疊加權(quán)重wi,jk設(shè)置為1.
干涉圖相位疊加的數(shù)學(xué)模型可表述為:
經(jīng)過(guò)干涉疊加后的累積相位值 φcum可表示為:
時(shí)間累積值tcum可表示為:
根據(jù)平均相干值給干涉圖設(shè)定權(quán)重wi,jk,相位累積值 φcum可表示為:
時(shí)間累積值tcum可表示為:
多幅主影像的干涉圖加權(quán)疊加方法得到的平均形變相位變化速率可表示為:
其中,M為主影像的個(gè)數(shù),N為第i幅主影像參與生成的干涉圖的個(gè)數(shù);φki,j與wki,j分別為在第i幅主影像參與生成的第j幅差分干涉圖中第k個(gè)像元的解纏形變相位與其疊加權(quán)重;ΔTi,j為第i幅主影像參與生成的第j幅差分干涉圖中主影像與輔影像之間的時(shí)間基線;λ為雷達(dá)波的波長(zhǎng).
由于SAR 影像中的大氣延遲相位是隨機(jī)分布,干涉像對(duì)經(jīng)差分處理生成的干涉圖中大氣延遲相位也隨機(jī)分布.根據(jù)誤差傳播定律,在對(duì)所有干涉圖進(jìn)行疊加后,高質(zhì)量的相干點(diǎn)受到的大氣延遲影響可表示為:
將所有干涉圖加權(quán)疊加后,對(duì)于高質(zhì)量相干點(diǎn),大氣延遲對(duì)線性形變速率的影響可表示為:
由式(8)和式(10)可得,像元k經(jīng)過(guò)干涉圖加權(quán)疊加后的平均速率可表示為:
從式(11)可以看出,疊加后平均形變速率中,形變相位信息與大氣延遲間的信噪比得到了提高.
為了驗(yàn)證方法的有效性,收集相干圖制成相干圖集,如圖3所示.根據(jù)前文所述設(shè)置權(quán)重的理論,計(jì)算相干圖集中的平均相干值,將相干性低的權(quán)重設(shè)置為0.因?yàn)檠芯繀^(qū)域左下角是低相干性的大片水域,在相干圖中顯示為黑色,因此只需研究主城區(qū)的相干性.對(duì)比圖3中的11 幅圖可知,相干圖028836-029711 在主城區(qū)的黑色區(qū)域最少,相干圖027961-030236 在主城區(qū)的黑色區(qū)域最多,統(tǒng)計(jì)其像素點(diǎn)相干值的分布,如圖4所示,并計(jì)算出像對(duì)028836-029711的像素點(diǎn)平均相干值為0.50,該像對(duì)相干性最好;像對(duì)027961-030236的像素點(diǎn)平均相干值為0.37,該像對(duì)相干性最差,且占相干性最好像對(duì)的平均相干值的比重小于3/4,因此需將該相干圖的權(quán)重設(shè)置為0,避免其相干圖參與疊加解算,降低最終的監(jiān)測(cè)精度.
圖3 相干圖集
圖4 像素點(diǎn)相干值分布圖
將以上符合要求的相位解纏圖疊加解算,解算得到平均沉降速率圖并疊加在Google Earth 上,如圖5所示.從圖5可以看出,在監(jiān)測(cè)過(guò)程中,該地區(qū)的地表整體上處于穩(wěn)定狀態(tài),這與選取的SAR 影像時(shí)間間隔較短有一定關(guān)系.結(jié)果表明在圣迭戈縣中心區(qū)域的地表有較大抬升,年均地面抬升速率在10 mm/a 左右,局部區(qū)域的最大抬升速率能達(dá)到30 mm/a.在圣迭戈縣的北部以及周邊區(qū)域都呈現(xiàn)出了不同程度的沉降,大部分區(qū)域的年均地面沉降速率在-16 mm/a,引起沉降的主要原因是該地區(qū)持續(xù)的農(nóng)業(yè)灌溉和城市建設(shè)用水來(lái)源于含水層儲(chǔ)水,導(dǎo)致對(duì)地下水過(guò)度開(kāi)采,造成了地下水位的下降,含水層透支且難以恢復(fù),從而進(jìn)一步惡化地面沉降.后來(lái)由于當(dāng)?shù)貙?duì)于地下水的一系列政策保護(hù),總體上城區(qū)的形變基本保持穩(wěn)定,與歷史監(jiān)測(cè)結(jié)果相比,地表形變速率具有明顯減緩的特征,說(shuō)明該地區(qū)的地面沉降已經(jīng)得到緩解.
圖5 年平均形變速率圖
為了驗(yàn)證結(jié)果的準(zhǔn)確性,采用傳統(tǒng)的D-InSAR 方法對(duì)2019年7月4日和2019年12月31日兩景SAR數(shù)據(jù)進(jìn)行形變監(jiān)測(cè),結(jié)果如圖6所示.從圖6中可以看出,傳統(tǒng)的D-InSAR 方法得到的監(jiān)測(cè)結(jié)果并不準(zhǔn)確.以左下角區(qū)域?yàn)槔M(jìn)行分析,左下角實(shí)驗(yàn)區(qū)域?yàn)楹F矫?是低相干區(qū)域,卻產(chǎn)生了不符實(shí)際的異常形變,監(jiān)測(cè)誤差較大.而本文提出的多幅主影像加權(quán)疊加方法結(jié)果相較而言更加準(zhǔn)確.此外,傳統(tǒng)的D-InSAR 方法是監(jiān)測(cè)起始時(shí)間之間的累積形變量作為最終的形變結(jié)果,用這種結(jié)果代表監(jiān)測(cè)期間的平均形變速率有失結(jié)果準(zhǔn)確性和代表性.本文通過(guò)選取多幅主影像,在采用D-InSAR 處理流程的基礎(chǔ)上,對(duì)得到的干涉圖集進(jìn)行加權(quán)疊加.這不僅可以代表監(jiān)測(cè)時(shí)間序列的沉降速率,又可同時(shí)獲得時(shí)間序列內(nèi)部包含各個(gè)時(shí)間段的形變速率,結(jié)果更具代表性.
圖6 傳統(tǒng)D-InSAR 監(jiān)測(cè)結(jié)果圖
地下水監(jiān)測(cè)結(jié)合了不同類(lèi)型的水位觀測(cè)測(cè)量(周期的、連續(xù)的和實(shí)時(shí)的),并結(jié)合這些記錄來(lái)提供井底隨時(shí)間變化的響應(yīng)的評(píng)估,所以使用USGS 網(wǎng)站地下水網(wǎng)絡(luò)數(shù)據(jù)進(jìn)行驗(yàn)證.圖7(a)是多幅主影像加權(quán)疊加方法在地下水監(jiān)測(cè)點(diǎn)的時(shí)間序列變化圖,可以看出該點(diǎn)的形變量先減小再增加,這意味著該區(qū)域地表先沉降又抬升,本文采用SPSS 軟件中的Analyse 模塊對(duì)地下水位與地面沉降的關(guān)系進(jìn)行分析,分析得到地下水水位與地面沉降量相關(guān)系數(shù)為0.941,兩者具有非常高的相關(guān)程度,如圖7(b).以地下水監(jiān)測(cè)點(diǎn)為基礎(chǔ)建立深層地下水位與地面沉降量關(guān)系的數(shù)據(jù)見(jiàn)表3,可以看出本文方法得到的地表形變趨勢(shì)與地下水水位的變化基本保持一致,兩者最大誤差為1.13 mm,由此證明本文方法用于監(jiān)測(cè)地表形變是可靠的,監(jiān)測(cè)精度可達(dá)毫米級(jí)別.
表3 地下水位監(jiān)測(cè)點(diǎn)水位與地面沉降數(shù)據(jù)表 (mm)
圖7 地下水?dāng)?shù)據(jù)精度驗(yàn)證
因城市地下水的抽取而引起地下水位的變化不是產(chǎn)生地表形變的唯一因素,同時(shí)還有人類(lèi)活動(dòng)和其他因素也會(huì)對(duì)地表形變?cè)斐捎绊?使得實(shí)驗(yàn)結(jié)果在數(shù)值上同地下水觀測(cè)數(shù)據(jù)具有一定的偏差.
本文在傳統(tǒng)D-InSAR 監(jiān)測(cè)方法的基礎(chǔ)上,提出一種基于較少影像實(shí)現(xiàn)時(shí)間序列監(jiān)測(cè)的多幅主影像干涉圖加權(quán)疊加方法.
該方法首先綜合考慮了垂直空間基線、時(shí)間基線和多普勒質(zhì)心頻率差異的影響,利用綜合函數(shù)模型選取三幅主影像,削弱主影像大氣延遲影響,為后續(xù)干涉圖加權(quán)疊加過(guò)程提供更多數(shù)量的相干圖,然后對(duì)相干圖質(zhì)量(依據(jù)平均相干值)進(jìn)行定權(quán),解決低相干影像圖對(duì)監(jiān)測(cè)結(jié)果的負(fù)面影響,通過(guò)加權(quán)疊加后,提高了相位圖的形變信息與大氣噪聲之間的信噪比,最后提取形變信息,實(shí)現(xiàn)對(duì)目標(biāo)區(qū)域的時(shí)序形變監(jiān)測(cè).
綜上所述,本文提出的方法通過(guò)采用圣迭戈縣城區(qū)的9 景Sentinel-1A 影像與傳統(tǒng)監(jiān)測(cè)方法進(jìn)行對(duì)比分析,實(shí)驗(yàn)結(jié)果表明,本文所提方法所得地表形變監(jiān)測(cè)結(jié)果與地下水監(jiān)測(cè)數(shù)據(jù)的線性擬合程度為0.941,監(jiān)測(cè)結(jié)果明顯優(yōu)于傳統(tǒng)D-InSAR 方法得到的結(jié)果.這種方法彌補(bǔ)了D-InSAR 技術(shù)不能進(jìn)行時(shí)間序列監(jiān)測(cè)的缺點(diǎn),解決了時(shí)序InSAR 技術(shù)對(duì)數(shù)據(jù)量要求較多的問(wèn)題,降低了傳統(tǒng)Stacking 方法中低相干影像參與疊加帶來(lái)的不利影響,提高了監(jiān)測(cè)精度,方法有效可行.同時(shí),本文方法在長(zhǎng)時(shí)間序列形變監(jiān)測(cè)中具有重要研究意義,后續(xù)將會(huì)收集更多資料進(jìn)行進(jìn)一步研究.