北京東部平原區(qū)地面沉降時(shí)空演化特征及預(yù)測(cè)

于文， 宮輝力， 陳蓓蓓， 周超凡

(1.首都師范大學(xué)水資源安全北京實(shí)驗(yàn)室，北京 100048； 2.首都師范大學(xué)地面沉降機(jī)理與防控教育部重點(diǎn)實(shí)驗(yàn)室，北京 100048； 3.首都師范大學(xué)城市環(huán)境過(guò)程與數(shù)字模擬國(guó)家重點(diǎn)實(shí)驗(yàn)室培育基地，北京 100048； 4.首都師范大學(xué)京津冀平原地下水與地面沉降國(guó)家野外科學(xué)觀測(cè)研究站，北京 100048)

0 引言

地面沉降是在自然和人為因素作用下地面高程降低的地質(zhì)現(xiàn)象。已成為城市發(fā)展中普遍存在的環(huán)境地質(zhì)問(wèn)題，同時(shí)也是制約社會(huì)經(jīng)濟(jì)可持續(xù)發(fā)展的重要地質(zhì)災(zāi)害之一。我國(guó)的區(qū)域地面沉降不可小覷。據(jù)《全國(guó)地面沉降防治規(guī)劃(2011—2020年)》統(tǒng)計(jì)，已有50多個(gè)城市發(fā)生地面沉降，大于200 mm地面沉降的區(qū)域面積超過(guò)7萬(wàn)km2，其中華北平原、長(zhǎng)江三角洲和汾渭盆地是其中的“重災(zāi)區(qū)”[1]。北京作為首都，城市發(fā)展迅速，近年來(lái)，地面沉降作為北京平原區(qū)的主要地質(zhì)災(zāi)害之一，影響著城市社會(huì)與經(jīng)濟(jì)的健康平穩(wěn)發(fā)展，其潛在的危害和經(jīng)濟(jì)損失已經(jīng)受到社會(huì)和政府的廣泛關(guān)注[2]。然而，傳統(tǒng)監(jiān)測(cè)地面沉降的精密水準(zhǔn)測(cè)量、分層標(biāo)技術(shù)以及全球定位系統(tǒng)等技術(shù)雖然精度較高，但是設(shè)備造價(jià)高、空間分辨率低和難以獲取大范圍信息等問(wèn)題給地面沉降監(jiān)測(cè)工作帶來(lái)困難。20世紀(jì)80年代出現(xiàn)的合成孔徑雷達(dá)干涉測(cè)量技術(shù)(interferometric synthetic aperture Radar，InSAR)彌補(bǔ)了傳統(tǒng)測(cè)量的不足，該技術(shù)具有全天時(shí)、全天候、監(jiān)測(cè)范圍廣和監(jiān)測(cè)精度高的特點(diǎn)[3]，已被廣泛應(yīng)用于各個(gè)領(lǐng)域，如地形測(cè)繪[4]、地震[5]和地面沉降[6]等。

InSAR技術(shù)是空間對(duì)地觀測(cè)技術(shù)的一次革命性飛躍，Gabriel等[7]1989年首次論證合成孔徑雷達(dá)差分干涉測(cè)量(differential interferometric SAR，D-InSAR)技術(shù)可用于探測(cè)亞厘米級(jí)的微小地表形變。但由于D-InSAR技術(shù)存在時(shí)空失相干[8]，相位解纏中參數(shù)值評(píng)估[9]以及大氣相位影響[10]等自身局限性,為了在一定程度上解決常規(guī)D-InSAR在空間、時(shí)間去相關(guān)、大氣誤差、軌道誤差和地形誤差去除中的局限性，多時(shí)相合成孔徑差分干涉測(cè)量技術(shù)(multi-temporal InSAR, MT-InSAR)陸續(xù)應(yīng)運(yùn)而生。其中，永久散射體(persistent scatters，PS)-InSAR 技術(shù)[11]僅有一幅主影像，其空間基線是將所有的從影像與唯一的主影像進(jìn)行配準(zhǔn)，空間基線相對(duì)較長(zhǎng)，在時(shí)間序列上識(shí)別具有穩(wěn)定散射特性的點(diǎn)，對(duì)目標(biāo)點(diǎn)經(jīng)過(guò)處理得到可靠的地表形變估算結(jié)果。

在地面沉降演化特征分析方面，已有相關(guān)研究，有關(guān)學(xué)者利用全局莫蘭指數(shù)和局部莫蘭指數(shù)對(duì)北京典型地面沉降區(qū)的格局特征進(jìn)行了分析，結(jié)果表明，區(qū)域地面沉降全局空間自相關(guān)十分顯著，局部空間格局分為高低聚集區(qū)與非相關(guān)區(qū)3個(gè)亞區(qū)[12]。此外有研究利用標(biāo)準(zhǔn)差橢圓方法對(duì)北京東部地區(qū)演化特征進(jìn)行分析，研究結(jié)合水文地質(zhì)資料，分析地面沉降的空間演化模式。結(jié)果表明，直到2012年，地面沉降主要向西北-東南方向發(fā)展，然后向各個(gè)方向擴(kuò)展[13]。有研究提出用區(qū)域重心法來(lái)分析來(lái)廣營(yíng)的地面沉降時(shí)空演化特征，得到地面沉降遷移方向和距離[14]?，F(xiàn)有研究在地面沉降矩陣場(chǎng)整體時(shí)空特性的研究較少，相關(guān)演化特征分析的方法仍存在不足，需要進(jìn)一步探索新方法。

地面沉降作為一種緩變的城市地質(zhì)災(zāi)害會(huì)對(duì)城市基礎(chǔ)設(shè)施造成損毀，嚴(yán)重威脅著城市安全并造成巨大的經(jīng)濟(jì)損失。建立一個(gè)高效的地面沉降預(yù)測(cè)模型對(duì)于地面沉降的防治和保障城市安全有著重要意義[15]。在地面沉降預(yù)測(cè)研究方面，主要分為物理模型和數(shù)學(xué)模型。物理模型主要通過(guò)水土耦合模型建立地面沉降預(yù)測(cè)模型[16]，這種方法一般需要大量的先驗(yàn)知識(shí)和實(shí)測(cè)數(shù)據(jù)，在大型的工程項(xiàng)目中被廣泛應(yīng)用。由于一些構(gòu)建地面沉降物理模型的土體參數(shù)很難獲取，并且這種模型的計(jì)算效率較低，因此很大程度上限制了這種模型的使用范圍。數(shù)學(xué)模型是利用離散的時(shí)序歷史沉降數(shù)據(jù)的數(shù)學(xué)統(tǒng)計(jì)規(guī)律建立地面沉降預(yù)測(cè)模型[17]。相對(duì)于物理模型，該數(shù)學(xué)模型所需數(shù)據(jù)更容易獲取且運(yùn)算高效。近年來(lái)，地理人工智能快速發(fā)展，一些智能算法被廣泛應(yīng)用于地面沉降。有關(guān)研究利用灰色馬爾科夫模型(grey-Markov model，IGMM)建立了北京平原的地面沉降預(yù)測(cè)模型，發(fā)現(xiàn)該預(yù)測(cè)模型精度較高[18]。隨著研究的深入，深度學(xué)習(xí)越來(lái)越受關(guān)注。有相關(guān)學(xué)者分別利用長(zhǎng)短期記憶(long short term memory，LSTM)模型[19]與循環(huán)神經(jīng)網(wǎng)絡(luò)(recurrent neural network，RNN)模型[20]對(duì)滄州市和撫順市地面沉降進(jìn)行預(yù)測(cè)。其中RNN運(yùn)算中每個(gè)隱含細(xì)胞單元計(jì)算最后需要經(jīng)過(guò)一個(gè)非線性函數(shù)，輸出[0,1]間的結(jié)果，這使多次運(yùn)算后的數(shù)值不斷衰減且無(wú)法記憶較遠(yuǎn)位置數(shù)據(jù)。LSTM是一種特殊的RNN，解決了RNN運(yùn)算的部分局限性，其通過(guò)在隱藏層增加門機(jī)制來(lái)控制信息的流失，再經(jīng)過(guò)反向傳播過(guò)程的動(dòng)態(tài)調(diào)整，使得網(wǎng)絡(luò)可以學(xué)習(xí)長(zhǎng)距離的時(shí)間序列數(shù)據(jù)。盡管該深度學(xué)習(xí)算法取得一定進(jìn)展，但在地面沉降時(shí)序預(yù)測(cè)上仍面臨著挑戰(zhàn)。

本文首先利用PS-InSAR技術(shù)對(duì)北京東部平原區(qū)地面沉降信息進(jìn)行獲取。鑒于當(dāng)前對(duì)地面沉降要素場(chǎng)區(qū)域整體的研究有所不足，運(yùn)用經(jīng)驗(yàn)正交函數(shù)(empirical orthogonal function，EOF)對(duì)地面沉降信息的空間特征與時(shí)序特征進(jìn)行了系統(tǒng)分析與規(guī)律挖掘，并利用一種改進(jìn)的LSTM模型對(duì)區(qū)域地面沉降進(jìn)行預(yù)測(cè)，為城市健康發(fā)展提供依據(jù)。

1 研究區(qū)及其數(shù)據(jù)源

1.1 研究區(qū)概況

北京位于華北平原北部邊緣， 地理位置在E 115° 20′～117° 33′，N39° 23′～41° 05′之間。研究區(qū)(圖1)位于北京平原東部，是地面沉降嚴(yán)重的地區(qū)，是典型的溫帶大陸性季風(fēng)氣候，年均氣溫為11～12 ℃，四季分明。春季多風(fēng)，降雨較少，易發(fā)生干旱； 夏季降雨較多且溫度較高； 秋季秋高氣爽，溫度適中，降雨適中； 冬季盛行西北風(fēng)，寒冷干燥，降雨較少。根據(jù)歷史資料記載，北京西單和東單一帶最早在1935年就發(fā)生了地面沉降，建國(guó)之后，隨著城市的不斷發(fā)展，地面沉降的范圍也在逐步擴(kuò)大。從沉降發(fā)育歷史來(lái)看，北京平原區(qū)地面沉降先后經(jīng)歷了形成階段(1955—1973年)、發(fā)展階段(1973—1983年)、擴(kuò)展階段[21](1983—1999年)、快速發(fā)展階段[22](1999—2014年)、區(qū)域發(fā)展不平衡階段[23](2014—2016年)。其中朝陽(yáng)、通州沉降區(qū)連成一片，成為北京平原區(qū)沉降最發(fā)育的地區(qū)。

圖1 研究區(qū)概況Fig.1 Overview of the study area

1.2 數(shù)據(jù)源

本研究中使用的數(shù)據(jù)集來(lái)自2顆不同的衛(wèi)星，一個(gè)是加拿大的RADARSAT-2衛(wèi)星，另一個(gè)是來(lái)自歐空局的Sentinel-1衛(wèi)星。其中，本研究選取48景存檔RADARSAT-2降軌數(shù)據(jù)，時(shí)間跨度為2010年11月22日—2015年11月20日； 選取61景存檔Sentinel-1升軌數(shù)據(jù)，時(shí)間跨度為2016年1月14日—2018年11月11日。數(shù)據(jù)集參數(shù)如表1所示。研究中利用Sarproz軟件來(lái)處理SAR數(shù)據(jù)集。在PS-InSAR處理過(guò)程中，使用空間分辨率為30 m的航天飛機(jī)雷達(dá)地形任務(wù)(shuttle Radar topography mission，SRTM)數(shù)字高程模型(digital elevation model，DEM)去除地形相位，并對(duì)干涉圖進(jìn)行地理編碼。

表1 S1A雷達(dá)影像信息情況Tab.1 S1A Radar image information

2 研究方法

2.1 PS-InSAR

PS-InSAR彌補(bǔ)了差分干涉測(cè)量的不足，為監(jiān)測(cè)城市微小形變帶來(lái)便利。該方法是由Ferretti等[11]在2000年提出。利用研究區(qū)域內(nèi)的多景SAR影像，通過(guò)分析幅度與相位信息，找出不受時(shí)間、空間基線失相關(guān)等因素影響的穩(wěn)定點(diǎn)目標(biāo)，并將經(jīng)過(guò)長(zhǎng)時(shí)間的推移，仍能保持穩(wěn)定散射特性的點(diǎn)稱為PS點(diǎn)。

該方法主要針對(duì)覆蓋同一個(gè)區(qū)域的N景SAR影像，選擇出最佳主影像，將其他的N-1景SAR影像與主影像進(jìn)行配準(zhǔn)，除去地形相位的影響，得到N-1幅差分干涉結(jié)果圖。則第i幅干涉圖的第j個(gè)點(diǎn)目標(biāo)的相位公式為：

(1)

差分干涉處理后，進(jìn)行PS點(diǎn)的選取，選取散射特性強(qiáng)且穩(wěn)定的像素作為PS點(diǎn)，可以是具有二面角和散射特性強(qiáng)的建筑物，如道路邊緣、橋梁、裸露的巖石等。然后對(duì)提取的PS點(diǎn)相位信息進(jìn)行濾波，去除大氣相位的影響，得到變形相位信息，從而得到地面沉降信息。

為了便于對(duì)從地面獲得的數(shù)據(jù)進(jìn)行比較和分析，需要根據(jù)雷達(dá)成像幾何將視線方向(line of sight,LOS)獲取數(shù)據(jù)轉(zhuǎn)換為垂直數(shù)據(jù)du，公式為：

(2)

式中：dlos為L(zhǎng)OS方向的變形。

2.2 EOF

EOF也稱為特征向量分析，是一種分析矩陣數(shù)據(jù)結(jié)構(gòu)特征并提取原始數(shù)據(jù)特征量的方法，可以用來(lái)分析變量場(chǎng)的結(jié)構(gòu)特征，包括空間模態(tài)和時(shí)間序列。Lorenz[25]在20世紀(jì)50年代首次將其引入氣象和氣候研究，現(xiàn)在被廣泛應(yīng)用于地球科學(xué)、水文學(xué)和其他學(xué)科[26-27]。EOF分析為時(shí)空分解，既反映空間特征，也體現(xiàn)時(shí)間變化。即

X=EOFm×mPCm×n，

(3)

式中：m為PS點(diǎn)數(shù)；n為時(shí)間月數(shù);PCm×n為主成分矩陣。具體算法為：

1)數(shù)據(jù)矩陣標(biāo)準(zhǔn)化預(yù)處理，得到一個(gè)數(shù)據(jù)矩陣Xm×n。

2)計(jì)算協(xié)方差矩陣，即

(4)

3)計(jì)算Am×m的特征根λ和特征向量Vm×m，兩者滿足

Am×mVm×m=Vm×mEm×m，

(5)

式中E是m×m維對(duì)角陣，即

(6)

一般特征根從大到小排列,每一個(gè)非零的特征根λ對(duì)應(yīng)一列特征向量值，也稱EOF。如λ1對(duì)應(yīng)的特征向量稱第一個(gè)EOF模態(tài)。

4)計(jì)算主成分。將EOF投影到原始矩陣，得到每個(gè)空間特征向量對(duì)應(yīng)的時(shí)間系數(shù)，即

(7)

5)計(jì)算方差貢獻(xiàn)率。X的方差大小可以用特征根λ表示，第k個(gè)模態(tài)對(duì)總的方差貢獻(xiàn)率P為：

(8)

2.3 Attention-LSTM模型

LSTM在傳統(tǒng) RNN 的基礎(chǔ)上引入3個(gè)門[28]： 遺忘門、輸入門、輸出門，結(jié)構(gòu)如圖2所示。

圖2 LSTM單元結(jié)構(gòu)示意圖Fig.2 Schematic diagram of LSTM unit structure

具體表達(dá)式為：

ft=σ(Wf·[Ht-1,Xt]+Bf)，

(9)

it=σ(Wi·[Ht-1,Xt]+Bi)，

(10)

(11)

(12)

ot=σ(Wo·[Ht-1,Xt]+Bo)，

(13)

Ht=ot·tanhCt，

(14)

Attention機(jī)制是一種模擬人腦注意力的模型，可以根據(jù)輸入的每項(xiàng)特征對(duì)輸出的影響為神經(jīng)網(wǎng)絡(luò)中的隱層狀態(tài)賦予不同的權(quán)重。本文將該機(jī)制引入地面沉降時(shí)序預(yù)測(cè)模型中，以選擇性地關(guān)注不同時(shí)間步的輸入對(duì)預(yù)測(cè)結(jié)果的影響，從而改善預(yù)測(cè)效果。Attention-LSTM的地面沉降預(yù)測(cè)框架如圖3所示。

圖3 Attention-LSTM的地面沉降預(yù)測(cè)框架Fig.3 Land subsidence prediction framework using Attention-LSTM

3 結(jié)果與分析

3.1 北京東部平原區(qū)地面沉降信息獲取

通過(guò)PS-InSAR方法對(duì)RARDARSAT-2與Sentinel-1數(shù)據(jù)處理，得到時(shí)間序列具有穩(wěn)定后向散射信號(hào)的PS點(diǎn)。利用ArcGIS軟件中克里金插值方法對(duì)PS點(diǎn)在不同時(shí)間的年沉降量屬性信息進(jìn)行插值處理。圖4為北京東部平原2011—2018年地面沉降年累計(jì)圖，可以看出該時(shí)段地面沉降在持續(xù)發(fā)生。截至2018年底，最大累計(jì)沉降量達(dá)1 030 mm。沉降區(qū)域主要分布在北京市朝陽(yáng)—通州交界處。

圖4-1 研究區(qū)累計(jì)沉降量Fig.4-1 Cumulative settlement of the study area

圖4-2 研究區(qū)累計(jì)沉降量Fig.4-2 Cumulative settlement of the study area

圖5為年沉降圖，其中可以看出2011—2018年間年沉降量，最大沉降發(fā)生在朝陽(yáng)金盞，其中2011—2015年平均沉降速率最大約為146 mm/a，2016—2018年平均沉降速率最大約為138 mm/a。通過(guò)時(shí)序年沉降圖，可以明顯看出2015年后，年沉降量有所減緩。為了對(duì)地面沉降結(jié)果進(jìn)行驗(yàn)證，其中以水準(zhǔn)點(diǎn)為中心，利用ArcGIS軟件進(jìn)行緩沖區(qū)分析，獲取所有水準(zhǔn)點(diǎn)200 m緩沖范圍，提取緩沖區(qū)內(nèi)的PS點(diǎn)，分別將各自緩沖范圍內(nèi)的PS點(diǎn)取平均值，然后與相應(yīng)的水準(zhǔn)點(diǎn)進(jìn)行驗(yàn)證，結(jié)果如圖6所示，得到較高的相關(guān)性值，其中2015—2016年RADARSAT-2結(jié)果與相應(yīng)時(shí)期水準(zhǔn)信息的相關(guān)系數(shù)達(dá)0.978； 2017—2018年Sentinel-1結(jié)果與相應(yīng)水準(zhǔn)信息的相關(guān)系數(shù)達(dá)0.954。

圖5 研究區(qū)時(shí)序年沉降量Fig.5 Time series annual settlement map of the study area

(a) 2015—2016年 (b) 2017—2018年圖6 InSAR結(jié)果與水準(zhǔn)結(jié)果驗(yàn)證Fig.6 Verification of InSAR results and leveling results

3.2 基于EOF的北京東部平原區(qū)典型地面沉降時(shí)空演化特征分析

分別對(duì)研究區(qū)2011—2018年時(shí)間段內(nèi)年尺度累計(jì)沉降變化與月尺度沉降變化進(jìn)行EOF分析。由于PS點(diǎn)較密，為了提高運(yùn)行效率，對(duì)點(diǎn)進(jìn)行抽稀預(yù)處理。在整個(gè)年變化時(shí)序處理中，使用9 852個(gè)點(diǎn)，圖7顯示了年沉降變化的空間模態(tài)(EOF-1和EOF-2)以及相關(guān)的主成分分量(PC1和PC2)?？臻g模態(tài)，也稱空間特征向量，反映地面沉降場(chǎng)的空間分布特征； 主成分對(duì)應(yīng)的是時(shí)間變化，也稱時(shí)間系數(shù)。EOF分析的前2種模態(tài)捕捉了長(zhǎng)期趨勢(shì)和季節(jié)變化，分別解釋了總方差的99.5%和0.4%。年沉降變化模態(tài)1的特征值為138 503 654，模態(tài)2的特征值為65 543。從圖7可以看出模態(tài)1的特征向量以朝陽(yáng)金盞、管莊和通州八里橋?yàn)橹行南蛲廨椛?，該中心點(diǎn)為高值中心，且是負(fù)值，該模態(tài)對(duì)應(yīng)的PC1與地面沉降變化的長(zhǎng)期趨勢(shì)和年際變化有關(guān)。結(jié)合振幅時(shí)間序列PC1，負(fù)高值中心與負(fù)時(shí)間系數(shù)變化疊加效應(yīng)為整體增加。結(jié)果表明，東部平原地面沉降年尺度變化在2011—2014年期間不斷增大，在2015—2016年期間沉降變化減緩，高值中心即紅色渲染部分附近的強(qiáng)度較其他區(qū)域變化幅度要大。PC2與主要的年際周期特征有關(guān)，峰值在2015年，谷值分別在2011年和2018年。結(jié)合模式2的空間格局，可以看出該模態(tài)存在2個(gè)大紅渲染區(qū)域正高值中心，位于朝陽(yáng)金盞—樓梓莊區(qū)域與通州臺(tái)湖村，1個(gè)負(fù)高值中心位于朝陽(yáng)化工橋，正高值中心與負(fù)高值中心區(qū)域出現(xiàn)相反的年際周期變化特征。

(a) 2011—2018年模態(tài)1 (b) 2011—2018年模態(tài)2

(c) 模態(tài)1對(duì)應(yīng)時(shí)間系數(shù)PC1 (d) 模態(tài)2對(duì)應(yīng)的時(shí)間系數(shù)PC2圖7 2011—2018年北京東部平原典型沉降區(qū)年尺度沉降的特征向量分布及對(duì)應(yīng)時(shí)間系數(shù)Fig.7 Distribution of subsidence fecture vectors and corresponding time coefficients of typical subsidence areas in the eastern plain of Beijing from 2011 to 2018

由于南水北調(diào)工程實(shí)施影響地面沉降演化[23],為了進(jìn)一步探究該影響，分別對(duì)2012—2014年的10 669個(gè)PS點(diǎn)與2016—2018年的11 914個(gè)PS點(diǎn)這2個(gè)時(shí)間段的月尺度沉降累計(jì)變化進(jìn)行EOF分析，圖8與圖9分別展示了2012—2014年與2016—2018年期間月尺度沉降變化的空間模式(EOF-1和EOF-2)以及相關(guān)的主成分分量(PC1和PC2)。2012—2014年期間，EOF-1與EOF-2分別解釋了總方差的98.8%和0.9%，特征值分別為25 768 800與250 655； 2016—2018年期間，EOF-1與EOF-2分別解釋了總方差的99.4%和0.2%，特征值分別為7 205 189與13 443。

(a) 2012—2014年模態(tài)1 (b) 2012—2014年模態(tài)2

(c) 模態(tài)1對(duì)應(yīng)時(shí)間系數(shù)PC1 (d) 模態(tài)2對(duì)應(yīng)的時(shí)間系數(shù)PC2圖8 2012—2014年北京東部平原典型沉降區(qū)月尺度沉降特征向量分布及對(duì)應(yīng)時(shí)間系數(shù)Fig.8 Distribution of monthly-scale subsidence feature vectors and corresponding time coefficients of typical subsidence area in the eastern plain of Beijing from 2012 to 2014

(a) 2016—2018年模態(tài)1 (b) 2016—2018年模態(tài)2

(c) 模態(tài)1對(duì)應(yīng)時(shí)間系數(shù)PC1 (d) 模態(tài)2對(duì)應(yīng)的時(shí)間系數(shù)PC2圖9 2016—2018年北京東部平原典型沉降區(qū)月尺度沉降的特征向量分布及對(duì)應(yīng)時(shí)間系數(shù)Fig.9 Distribution of monthly-scale subsidence feature vectors and corresponding time coefficients of typical subsidence area in the eastern plain of Beijing from 2016 to 2018

從圖8可以看出模態(tài)1的特征向量以朝陽(yáng)金盞、管莊和通州八里橋?yàn)橹行南蛲廨椛?，該中心點(diǎn)為高值中心，且是正值，該模態(tài)對(duì)應(yīng)的PC1與地面沉降變化的長(zhǎng)期趨勢(shì)保持一致。正高值中心與正值時(shí)間系數(shù)變化疊加效應(yīng)為整體增加。結(jié)合振幅時(shí)間序列PC1，結(jié)果表明，東部平原地面沉降月尺度變化在2012—2013年7月期間不斷增大，在2013年7月—2015年底減小，變化幅度最大也是高值中心處。PC2與主要的地面沉降季節(jié)性特征有關(guān)。PC2的時(shí)間系數(shù)變化振幅近似于正弦波，峰值在每年7—9月，谷值在1—3月。結(jié)合模式2的空間格局，可以看出很少一部分地區(qū)出現(xiàn)明顯的季節(jié)性周期變化，大多數(shù)地區(qū)顯示微弱的季節(jié)性周期變化。

從圖9可以看出月尺度模態(tài)1的特征向量高值中心區(qū)域范圍較2012—2014年有所縮減，時(shí)間系數(shù)值從±8 000減小到±6 000以內(nèi)。該高值中心為正值，該模態(tài)對(duì)應(yīng)的PC1與地面沉降變化的長(zhǎng)期趨勢(shì)有關(guān)。結(jié)合振幅時(shí)間序列PC1，結(jié)果表明，東部平原地面沉降月尺度變化在2016—2017年下旬不斷增大，在2017年下旬—2018年底減小，負(fù)值中心的變化模式相反。PC2與主要的季節(jié)性周期特征有關(guān)。PC2的振幅近似于正弦波，峰值在1—3月，谷值在7—9月。結(jié)合模式2的空間格局，可以看出正值高中心與負(fù)值高中心相間分布，均具有明顯的季節(jié)性周期變化。

3.3 基于Attention-LSTM模型的北京東部平原區(qū)地面沉降預(yù)測(cè)

分別利用傳統(tǒng)LSTM與Attention-LSTM方法對(duì)時(shí)序地面沉降進(jìn)行預(yù)測(cè)。發(fā)現(xiàn)添加Attention網(wǎng)絡(luò)機(jī)制后，預(yù)測(cè)精度得到提高。通過(guò)對(duì)LSTM模型增加Attention機(jī)制，可以充分學(xué)習(xí)不同沉降時(shí)間序列中的非線性關(guān)聯(lián)，進(jìn)而可以捕獲研究區(qū)域中的復(fù)雜沉降機(jī)理。

在模型訓(xùn)練時(shí)，通過(guò)不斷調(diào)節(jié)訓(xùn)練次數(shù)epoch來(lái)提高訓(xùn)練精度，圖10顯示epoch與損失函數(shù)loss的關(guān)系，可以看出，隨著epoch增大，損失函數(shù)均方根誤差不斷減小，并趨于穩(wěn)定?？梢钥闯鯝ttention-LSTM的最終損失函數(shù)均方根誤差優(yōu)于LSTM，其中Attention-LSTM的損失函數(shù)均方根誤差小于0.01，達(dá)到很高的精度。

(a) 2012—2015年LSTM (b) 2012—2015年Attention-LSTM(c) 2016—2018年LSTM(d) 2016—2018 年Attention-LSTM圖10 預(yù)測(cè)模型損失函數(shù)Fig.10 Forecast model loss function

分別對(duì)空間上點(diǎn)(包括PS-InSAR獲取的真實(shí)沉降點(diǎn)、LSTM預(yù)測(cè)點(diǎn)、Attention-LSTM預(yù)測(cè)點(diǎn))進(jìn)行克里金插值處理，得到累計(jì)沉降結(jié)果如圖11所示。從渲染色帶分析，三者結(jié)果的趨勢(shì)大體相同，但局部存在差異，在空間分布上，通過(guò)顏色拉伸結(jié)果，大體可以看出Attention-LSTM預(yù)測(cè)點(diǎn)的分布與PS-InSAR方法獲取地面沉降有一定的吻合度，而LSTM預(yù)測(cè)點(diǎn)與其他兩者局部地區(qū)有明顯的色調(diào)差異。

(a) 2012—2015年真實(shí)值 (b) 2012—2015年LSTM預(yù)測(cè) (c) 2012—2015年Attention-LSTM預(yù)測(cè)圖11-1 LSTM與Attention-LSTM的區(qū)域沉降預(yù)測(cè)結(jié)果Fig.11-1 Regional settlement prediction results of LSTM and Attention-LSTM

(d) 2016—2018年真實(shí)值 (e) 2016—2018年LSTM預(yù)測(cè) (f) 2016—2018年Attention-LSTM預(yù)測(cè)圖11-2 LSTM與Attention-LSTM的區(qū)域沉降預(yù)測(cè)結(jié)果Fig.11-2 Regional settlement prediction results of LSTM and Attention-LSTM

為了進(jìn)一步對(duì)預(yù)測(cè)結(jié)果的值進(jìn)行分析，在研究區(qū)域選擇了2條剖面，一條東西向剖面，一條南北向剖面，如圖11中所示。利用ArcGIS軟件 3D分析工具，對(duì)預(yù)測(cè)值柵格影像與真實(shí)沉降柵格影響剖面所在區(qū)域的值進(jìn)行提取(圖12)，對(duì)剖面所處位置可以看出剖面位置點(diǎn)真實(shí)值與預(yù)測(cè)值的變化趨勢(shì)保持一致性。但同時(shí)看出在剖面上的絕大多數(shù)點(diǎn)，Attention-LSTM預(yù)測(cè)更靠近真實(shí)測(cè)量值。

(a) 2012—2015年?yáng)|西向剖面真實(shí)值與預(yù)測(cè)值 (b) 2012—2015年南北向剖面真實(shí)值與預(yù)測(cè)值

(c) 2016—2018年?yáng)|西向剖面真實(shí)值與預(yù)測(cè)值 (d) 2016—2018年南北向剖面真實(shí)值與預(yù)測(cè)值圖12 研究區(qū)所選剖面真實(shí)值與預(yù)測(cè)值的對(duì)比Fig.12 Comparison between the real value and the predicted value of the selected profile in the study area

4 結(jié)論與展望

地下水的過(guò)度開采導(dǎo)致北京平原區(qū)域內(nèi)產(chǎn)生嚴(yán)重的地面沉降現(xiàn)象。隨著南水北調(diào)工程項(xiàng)目的實(shí)施，地面沉降演化特征將發(fā)生變化。本次研究選取北京平原區(qū)沉降最嚴(yán)重的區(qū)域，東部平原朝陽(yáng)—通州區(qū)域，探討該區(qū)域地面沉降演化情況，利用EOF模型分析地面沉降要素場(chǎng)的空間分布特征與時(shí)間變化。并利用深度學(xué)習(xí)方法對(duì)該區(qū)地面沉降進(jìn)行時(shí)間序列預(yù)測(cè)。研究結(jié)論如下：

1)2011—2018年，北京東部平原區(qū)地面沉降主要發(fā)生在朝陽(yáng)與通州交界處，截至2018年底，最大累計(jì)沉降量達(dá)1 030 mm，2011—2015年平均沉降速率最大約為146 mm/a，2016—2018年平均沉降速率最大約為138 mm/a。將PS-InSAR方法獲取的地面沉降結(jié)果與水準(zhǔn)結(jié)果進(jìn)行驗(yàn)證，得到較高的相關(guān)系數(shù)。

2)利用EOF方法對(duì)2011—2018年年尺度地面沉降場(chǎng)與月尺度地面沉降場(chǎng)的空間特性與時(shí)間變化特性進(jìn)行分析，發(fā)現(xiàn)研究區(qū)域空間模態(tài)1方差貢獻(xiàn)率很大，幾乎代表研究區(qū)域空間的整體演化情況。對(duì)應(yīng)時(shí)間系數(shù)線性趨勢(shì)明顯。模態(tài)2有一定的方差貢獻(xiàn)率，但占比很小，對(duì)應(yīng)的時(shí)間系數(shù)季節(jié)性顯著。

3)利用Attention-LSTM對(duì)地面沉降進(jìn)行預(yù)測(cè)，得到精度較高的預(yù)測(cè)模型，擴(kuò)展了深度學(xué)習(xí)在地面沉降研究方面的應(yīng)用，利用單一變量對(duì)地面沉降進(jìn)行預(yù)測(cè)，避免了由于數(shù)據(jù)源不足的困境。

在以后的研究中，訓(xùn)練模型時(shí)，將盡可能考慮影響地面沉降的機(jī)理要素，并將其考慮到時(shí)間序列預(yù)測(cè)中。

志謝：感謝Sarproz軟件的制造商以及歐空局提供的Sentinel-1數(shù)據(jù)。