利用DS-InSAR技術(shù)監(jiān)測(cè)黃河三角洲地表形變

曹建濤 鄭翔元 范洪冬 李國(guó)華 黃 晨

1 中國(guó)礦業(yè)大學(xué)自然資源部國(guó)土環(huán)境與災(zāi)害監(jiān)測(cè)重點(diǎn)實(shí)驗(yàn)室，江蘇省徐州市大學(xué)路1號(hào)，221116

2 山東省國(guó)土測(cè)繪院，濟(jì)南市經(jīng)十路2301號(hào)，250013

黃河流域是我國(guó)重要的生態(tài)屏障和經(jīng)濟(jì)地帶，在推動(dòng)我國(guó)社會(huì)經(jīng)濟(jì)發(fā)展與生態(tài)建設(shè)方面都具有非常重要的作用。黃河三角洲作為我國(guó)最大的三角洲，蘊(yùn)含著豐富的鹵水、石油、天然氣等自然資源。隨著人們對(duì)地下鹵水的抽取及石油的不斷開(kāi)采，黃河三角洲存在多處不同程度的地表沉降區(qū)，嚴(yán)重影響當(dāng)?shù)鼐用竦纳a(chǎn)生活。因此，準(zhǔn)確監(jiān)測(cè)黃河三角洲地區(qū)地表沉降對(duì)于減少該地區(qū)的自然災(zāi)害和經(jīng)濟(jì)損失具有重要意義[1]。

目前已有研究人員利用水準(zhǔn)測(cè)量和全球定位系統(tǒng)對(duì)黃河三角洲地區(qū)地表沉降進(jìn)行監(jiān)測(cè)，這些測(cè)量方法雖有較高的精度，但存在費(fèi)時(shí)費(fèi)力、工作效率較低、地面測(cè)量點(diǎn)位易被破壞、測(cè)量點(diǎn)密度低、測(cè)站點(diǎn)不一定與沉降中心重合、難以進(jìn)行大范圍監(jiān)測(cè)等諸多不足。2000年以來(lái)，永久散射體干涉測(cè)量(persistent scatterers InSAR, PS-InSAR)、短基線集干涉測(cè)量(small baseline subset, SBAS)等時(shí)序InSAR技術(shù)已得到較為廣泛的應(yīng)用。有學(xué)者利用PS-InSAR或SBAS方法對(duì)黃河三角洲區(qū)域進(jìn)行地表沉降監(jiān)測(cè)[2]，雖然研究人員將InSAR監(jiān)測(cè)結(jié)果與CORS站結(jié)果進(jìn)行了精度驗(yàn)證[2-3]，但該類(lèi)時(shí)序InSAR方法多依賴于具有強(qiáng)散射特性的地物目標(biāo)，而在裸地和多植被覆蓋區(qū)域，其選點(diǎn)密度不高，難以反映沉降區(qū)域詳細(xì)的沉降特征，降低了監(jiān)測(cè)結(jié)果的準(zhǔn)確性。此外，現(xiàn)有研究監(jiān)測(cè)時(shí)段較早，無(wú)法反映黃河三角洲沉降的現(xiàn)狀。

為解決常規(guī)時(shí)序InSAR的局限性，F(xiàn)erretti等[4]提出SqueeSAR方法，將相干性高的DS點(diǎn)與PS點(diǎn)聯(lián)合進(jìn)行形變解算，自此揭開(kāi)了分布式目標(biāo)InSAR(distributed scatterers InSAR，DS-InSAR)方法研究的序幕。在此基礎(chǔ)上DS-InSAR不斷發(fā)展，2015年Cao等[5]詳細(xì)論述各種相位優(yōu)化方法的原理，并將相干矩陣采用特征值分解方法實(shí)現(xiàn)相位優(yōu)化過(guò)程；2016年蔣彌等[6]提出建立在快速同質(zhì)點(diǎn)選取下的干涉處理框架，利用置信區(qū)間估計(jì)來(lái)判斷像元的同質(zhì)點(diǎn)，得到美國(guó)加州Lost Hills油田2002～2004年間的地表形變監(jiān)測(cè)數(shù)據(jù)，表明該方法在運(yùn)算效率和增加空間觀測(cè)密度上具有明顯優(yōu)勢(shì)；2019年Zhao等[7]利用非線性相位優(yōu)化，并基于StaMPS算法將PS點(diǎn)與優(yōu)化后DS點(diǎn)聯(lián)合處理，對(duì)美國(guó)圖森地區(qū)進(jìn)行地表形變監(jiān)測(cè)，有效提高了相干點(diǎn)密度，但計(jì)算效率方面優(yōu)勢(shì)不大。目前DS-InSAR方法在山體滑坡、煤礦開(kāi)采區(qū)[8]、鐵路沿線、水利設(shè)施[9]等領(lǐng)域均得到成功應(yīng)用，這些成功應(yīng)用案例表明，DS-InSAR能夠提高地表形變監(jiān)測(cè)點(diǎn)的密度及低相干性區(qū)域形變監(jiān)測(cè)的精度。

考慮黃河三角洲大范圍區(qū)域存在大量非城市區(qū)域的狀況，本文利用DS-InSAR方法獲取2019-12～2020-12期間黃河三角洲區(qū)域的地表形變監(jiān)測(cè)數(shù)據(jù)，并與PS-InSAR結(jié)果進(jìn)行對(duì)比分析和交叉驗(yàn)證，證明結(jié)果的可靠性；最后，根據(jù)該地區(qū)人類(lèi)活動(dòng)特征及沉降中心分布情況分析黃河三角洲地表沉降的影響因素。

1 研究區(qū)域及數(shù)據(jù)概況

1.1 區(qū)域概況

黃河三角洲位于我國(guó)山東省東營(yíng)市，北倚渤海灣，東鄰萊州灣?，F(xiàn)代黃河三角洲是1855年黃河改道攜泥沙匯入渤海后，泥沙淤積形成的扇形沖積平原，面積約為5 450 km2，是我國(guó)最大的三角洲。該地區(qū)地勢(shì)平坦，海拔高度多在6 m以下?，F(xiàn)代黃河三角洲植被覆蓋范圍廣，濕地覆蓋面積約4 500 km2，且具有豐富的石油、天然氣等自然資源。在特殊的地質(zhì)構(gòu)造下，現(xiàn)代黃河三角洲存在約60%面積的軟土層，且分布不均勻。軟土分布區(qū)含水量高、承載力低，因此易發(fā)生地表沉降。

本文研究區(qū)覆蓋河口區(qū)、利津縣及墾利縣的大部分區(qū)域，包含勝利油田(中國(guó)第二大油田)、孤東油田、東營(yíng)港口、東營(yíng)機(jī)場(chǎng)及大量鹽田區(qū)。由于黃河三角洲地區(qū)特殊的地質(zhì)構(gòu)造特點(diǎn)及大量開(kāi)采油田、地下鹵水等人類(lèi)活動(dòng)導(dǎo)致地層壓力發(fā)生變化，該地區(qū)存在多個(gè)沉降區(qū)，且沉降量級(jí)大，嚴(yán)重影響了當(dāng)?shù)厝嗣竦纳a(chǎn)生活。

1.2 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)選取覆蓋研究區(qū)的26景Sentinel-1A SAR影像，時(shí)間跨度為2019-12-19～2020-12-13，成像模式為IW模式，極化方式為VV極化，距離向和方位向的像元尺寸分別為2.33 m和13.93 m。為去除地形相位，本文采用90 m分辨率的SRTM數(shù)字高程模型(DEM)作為外部DEM數(shù)據(jù)。實(shí)驗(yàn)以2020-04-29影像作為公共主影像，分別與其他影像生成25個(gè)干涉對(duì)。

2 DS-InSAR方法原理

2.1 FaSHPS同質(zhì)點(diǎn)識(shí)別

FaSHPS(fast SHP selection)算法[6]的核心思想是將假設(shè)檢驗(yàn)問(wèn)題轉(zhuǎn)化為置信區(qū)間估計(jì)，通過(guò)影像像元數(shù)據(jù)在時(shí)間序列上與參考點(diǎn)的相似性來(lái)判斷是否為同質(zhì)像元點(diǎn)。本文利用FaSHPS方法逐個(gè)計(jì)算出像元周?chē)耐|(zhì)點(diǎn)數(shù)量，以20個(gè)同質(zhì)點(diǎn)作為閾值，將同質(zhì)點(diǎn)數(shù)量高于該閾值的參考像元作為DS預(yù)選點(diǎn)。

(1)

式中，P{·}表示概率，z1-α/2表示在標(biāo)準(zhǔn)正態(tài)分布中1-α/2分位點(diǎn)，μ(S)為像元S的期望，var(A(S))表示像元S在時(shí)間序列上振幅的真實(shí)方差，N表示SAR影像的個(gè)數(shù)。

在同質(zhì)區(qū)域可認(rèn)為SAR影像的振幅服從瑞利分布[10]，則變異系數(shù)CV可由下式計(jì)算：

(2)

式中，E(·)表示期望。假定像元的散射特性在觀測(cè)時(shí)間范圍內(nèi)基本穩(wěn)定，則式(1)可轉(zhuǎn)換為：

(3)

2.2 特征值分解相位優(yōu)化

分布式目標(biāo)像元內(nèi)地物應(yīng)具有相同的散射特性，由于受到噪聲及時(shí)空失相干的影響，時(shí)間序列SAR影像干涉圖每個(gè)像元的信號(hào)矢量中包含多種類(lèi)型的散射信號(hào)，使得分布式目標(biāo)像元的相位穩(wěn)定性較差。通過(guò)對(duì)像元的相位優(yōu)化可達(dá)到對(duì)干涉圖像元相位降噪的目的。像元相干矩陣對(duì)應(yīng)的特征值代表不同的后向散射特性，將相干矩陣進(jìn)行特征值分解[5]即可得到最大特征值及其對(duì)應(yīng)的特征向量，認(rèn)為此特征向量對(duì)應(yīng)的像元散射特性是穩(wěn)定的，并將其作為優(yōu)化后的相位。

假設(shè)在勻質(zhì)區(qū)域Ω內(nèi)部存在NP個(gè)后向散射特性相近的像元，其相干矩陣T可表示為：

(4)

式中，y=[y1,y2,…,yN]為分布式目標(biāo)的同質(zhì)點(diǎn)在N景SAR影像上的復(fù)數(shù)觀測(cè)量經(jīng)歸一化處理后的復(fù)數(shù)向量，(·)H表示矩陣的共軛轉(zhuǎn)置。上式得到的相干矩陣為半正定Hermitian矩陣，經(jīng)特征值分解可得：

(5)

式中，Λ為非負(fù)實(shí)數(shù)特征值λi的對(duì)角矩陣，U為不同特征值對(duì)應(yīng)的正交特征向量。λi越大，對(duì)應(yīng)的地物散射機(jī)制就越占優(yōu)勢(shì)，因此，可將最大特征值λ1對(duì)應(yīng)的特征向量μ1作為主散射體對(duì)應(yīng)的相位。由此可計(jì)算主散射體信號(hào)Tsignal，其余的作為去相干噪聲信號(hào)Tnoise。將去相干噪聲信號(hào)的相位分量去除，保留主散射體信號(hào)的相位分量，即可達(dá)到優(yōu)化相位的目的：

T=Tsignal+Tnoise=

(6)

相位優(yōu)化完成后，還需對(duì)相位的優(yōu)化質(zhì)量進(jìn)行評(píng)估。將優(yōu)化前后的干涉相位進(jìn)行差值擬合計(jì)算得到擬合度[4]：

(7)

式中，γDS為擬合度，還可看作為分布式目標(biāo)的時(shí)序相干性；φmn表示兩幅SAR影像優(yōu)化之前的干涉相位；φm和φn分別表示兩幅SAR影像優(yōu)化之后的相位。本文以0.4作為閾值，將相位優(yōu)化質(zhì)量高于閾值的像元選定為最終的DS點(diǎn)。

2.3 處理步驟

本文利用FaSHPS算法獲得DS候選點(diǎn)，經(jīng)特征值分解相位優(yōu)化及優(yōu)化質(zhì)量評(píng)估后確定最終DS點(diǎn)，將得到的DS點(diǎn)聯(lián)合振幅離差閾值初選得到的PS點(diǎn)，構(gòu)建用于相位分析的Delaunay三角網(wǎng)。由于DEM誤差的不確定性，利用StaMPS軟件首先對(duì)視角誤差的空間不相干部分進(jìn)行相位校正，然后在時(shí)間和空間兩個(gè)維度進(jìn)行干涉相位的3D解纏，利用時(shí)間上的高通濾波器及空間上的低通濾波器對(duì)具有時(shí)間和空間相關(guān)的相位進(jìn)行濾波，以估計(jì)剩余的空間相關(guān)干擾項(xiàng)，如大氣相位和軌道相位，在減去上述分量之后最終可提取黃河三角洲區(qū)域的地表時(shí)序形變。具體的處理流程如圖1所示。

3 InSAR監(jiān)測(cè)結(jié)果分析

3.1 地表形變監(jiān)測(cè)結(jié)果

分別采用PS-InSAR和DS-InSAR兩種方法提取黃河三角洲2019-12-19～2020-12-13時(shí)序地表沉降信息(圖2)。PS-InSAR方法共選取425 729個(gè)監(jiān)測(cè)點(diǎn)，最大沉降速率為-227 mm/a；而DS-InSAR方法共選取2 365 545個(gè)監(jiān)測(cè)點(diǎn)，約是PS-InSAR方法選點(diǎn)個(gè)數(shù)的5.56倍，最大沉降速率為-238 mm/a。可見(jiàn)，DS-InSAR方法顯著增加了地面監(jiān)測(cè)點(diǎn)的空間分布密度。圖2表明，在該研究時(shí)間段內(nèi)黃河三角洲區(qū)域存在嚴(yán)重的地表沉降現(xiàn)象，兩種方法均監(jiān)測(cè)到4處明顯的大量級(jí)沉降區(qū)域，且形變結(jié)果接近，具有較好的一致性。本文監(jiān)測(cè)得到的沉降區(qū)域與文獻(xiàn)[11]得到的2015～2018年黃河三角洲的沉降區(qū)域相吻合，且通過(guò)實(shí)地調(diào)研發(fā)現(xiàn)，這些沉降區(qū)域存在大量的鹽田池和采油機(jī)，表明本實(shí)驗(yàn)可準(zhǔn)確識(shí)別沉降區(qū)域。由于DS-InSAR方法選點(diǎn)數(shù)量增多，在沉降區(qū)域其監(jiān)測(cè)效果要優(yōu)于PS-InSAR方法，沉降邊界更加明顯。

3.2 監(jiān)測(cè)結(jié)果分析

由于黃河三角洲區(qū)域的地物覆蓋類(lèi)型多以農(nóng)田、鹽田與養(yǎng)殖池為主，導(dǎo)致PS-InSAR方法選點(diǎn)數(shù)量較少，使得地表形變的分布、特征、范圍不夠顯著；而融合分布式目標(biāo)的DS-InSAR方法有效提高了非城市區(qū)域的地面監(jiān)測(cè)點(diǎn)密度。本文監(jiān)測(cè)到的4個(gè)沉降區(qū)域主要分布在鹽田和油田附近，其中A、B、C區(qū)域均位于東營(yíng)市河口區(qū)，且分布較為集中；D區(qū)域位于東營(yíng)市墾利縣南部。

為分析PS-InSAR與DS-InSAR所獲取的地表沉降的差異性，對(duì)4個(gè)沉降區(qū)域分別選取了剖面線(圖3)，并繪制了對(duì)應(yīng)的地表沉降散點(diǎn)圖(圖4)?？梢钥闯觯珼S-InSAR有效提高了測(cè)量點(diǎn)空間分布密度，能夠更加直觀地反映沉降區(qū)域地表沉降變化趨勢(shì)；同時(shí)也明顯看出，兩者具有較高的一致性，驗(yàn)證了DS-InSAR結(jié)果的可靠性。

A沉降區(qū)中心位于東營(yíng)市河口區(qū)廣河村國(guó)星鹽場(chǎng)[12]，其最大年沉降速率為-134 mm/a。由于存在過(guò)度開(kāi)采地下鹵水的情況，導(dǎo)致鹽場(chǎng)附近出現(xiàn)較大范圍的沉降漏斗，圖4(a)顯示沿A′-A″線距A′點(diǎn)約5～12 km的范圍內(nèi)，地表沉降速率達(dá)到-75 mm/a。

B沉降區(qū)中心位于鹽田與油田混合區(qū)域，最大年沉降速率達(dá)-238 mm/a。圖4(b)為B′-B″線沉降速率散點(diǎn)圖，該區(qū)域年沉降速率較大，沉降現(xiàn)象最為嚴(yán)重。該區(qū)域被黃河故道、東營(yíng)港疏港高速及孤北水庫(kù)包圍，自2015年就已出現(xiàn)沉降現(xiàn)象，由于鹽田分布較為集中且開(kāi)采力度大，如今地表沉降速率明顯加快。

D沉降區(qū)中心位于墾利縣東營(yíng)機(jī)場(chǎng)附近的永豐鹽場(chǎng)和東營(yíng)景洪鹽化有限公司，最大年沉降速率達(dá)-122 mm/a。由圖4(e)可看出，剖面線上的3個(gè)沉降漏斗對(duì)應(yīng)于永豐鹽場(chǎng)周?chē)?個(gè)集中曬鹽區(qū)域，其中兩個(gè)區(qū)域地表沉降速率均已超過(guò)-80 mm/a。另外，D′-D″剖線端點(diǎn)D′附近位于東營(yíng)機(jī)場(chǎng)，其年沉降速率達(dá)到-30 mm/a。

3.3 DS-InSAR可靠性分析

為進(jìn)一步驗(yàn)證DS-InSAR方法的可靠性，本文通過(guò)對(duì)比分析PS-InSAR與DS-InSAR在一定范圍內(nèi)同名監(jiān)測(cè)點(diǎn)的沉降速率進(jìn)行交叉驗(yàn)證。將PS-InSAR監(jiān)測(cè)結(jié)果作為基準(zhǔn)，以每個(gè)PS點(diǎn)為中心，15 m為半徑，找到與PS點(diǎn)同名的DS點(diǎn)。一共得到144 772個(gè)PS-DS同名點(diǎn)，利用皮爾遜相關(guān)性分析算法得到同名點(diǎn)的沉降速率相關(guān)性，并進(jìn)一步計(jì)算誤差分布(圖5)。

由圖5(a)可知，本文采用的DS-InSAR與PS-InSAR方法監(jiān)測(cè)得到的地表沉降速率具有良好的一致性，其同名點(diǎn)沉降速率的皮爾遜相關(guān)系數(shù)為0.727。圖5(b)表明，沉降速率誤差分布趨近于高斯分布，且均值接近0，其中誤差值小于20 mm/a的同名點(diǎn)個(gè)數(shù)為128 149，占所選同名點(diǎn)個(gè)數(shù)的88.5%。因此，DS-InSAR在提高地面目標(biāo)點(diǎn)密度的同時(shí)，具有很好的可靠性。

4 黃河三角洲地表形變影響因素分析

黃河三角洲地區(qū)地下鹵水豐富，且含鹽量高于海水，研究區(qū)域內(nèi)的河口區(qū)、墾利縣、東營(yíng)區(qū)已建成多處鹽場(chǎng)并投入生產(chǎn)，這些工廠開(kāi)采地下鹵水用于工業(yè)制溴、制鹽和制造相關(guān)化工產(chǎn)品等，但長(zhǎng)期不合理開(kāi)采地下鹵水將會(huì)導(dǎo)致地下鹵水水位下降、局部海水入侵和地面裂縫，嚴(yán)重時(shí)則會(huì)導(dǎo)致地表沉降、發(fā)生地質(zhì)災(zāi)害。眾多企業(yè)大量開(kāi)采地下鹵水不僅致使其資源減少，而且形成淺層地下鹵水沉降漏斗。因此，過(guò)度開(kāi)采地下鹵水用于制鹽和化工產(chǎn)業(yè)是造成黃河三角洲地區(qū)地面沉降的主要原因之一。

本文監(jiān)測(cè)到的A、B、C、D沉降區(qū)中心均為鹽場(chǎng)及鹽田開(kāi)采區(qū)，其中C沉降區(qū)的中心為鹽田和油田綜合開(kāi)采區(qū)。圖6為疊加了鹽田區(qū)和居民區(qū)大致范圍的典型沉降區(qū)域的衛(wèi)星影像圖?？梢钥闯觯緦?shí)驗(yàn)監(jiān)測(cè)得到的沉降區(qū)大都分布在鹽田區(qū)域，說(shuō)明黃河三角洲區(qū)域開(kāi)采地下鹵水制鹽是造成地面沉降的主要因素之一。A沉降區(qū)域(國(guó)星鹽場(chǎng))附近有少量的居民區(qū)，未發(fā)現(xiàn)該區(qū)域地表形變對(duì)居民區(qū)造成明顯影響；最大沉降區(qū)域B緊鄰東營(yíng)疏港高速公路和仙河鎮(zhèn)居民區(qū)，長(zhǎng)期過(guò)度開(kāi)采地下鹵水必將對(duì)高速公路和居民區(qū)造成損害；C沉降區(qū)域?yàn)辂}田和油田的綜合開(kāi)采區(qū)，從五號(hào)樁以南沿興港路兩側(cè)分布有大量的鹽田，且同時(shí)存在許多采油機(jī)，在二者的共同影響下造成大面積地表沉降，且附近建構(gòu)筑物出現(xiàn)裂縫現(xiàn)象；D沉降區(qū)域緊鄰紅光新村和東營(yíng)機(jī)場(chǎng)，目前其沉降速率緩慢，但永豐鹽場(chǎng)的鹵水開(kāi)采制鹽活動(dòng)對(duì)東營(yíng)機(jī)場(chǎng)附近的地基穩(wěn)定性造成一定程度的影響。

東營(yíng)市勝利油田作為中國(guó)第二大油田，其石油開(kāi)采活動(dòng)從上世紀(jì)開(kāi)始一直在進(jìn)行，長(zhǎng)期的深層石油開(kāi)采活動(dòng)會(huì)增加儲(chǔ)油層壓力，使得儲(chǔ)油層壓實(shí)收縮，從而導(dǎo)致緩慢的地表沉降，而人為進(jìn)行地下水回注能夠起到減緩地表沉降的作用，甚至?xí)l(fā)生小幅度地面抬升現(xiàn)象[13]。此外，人工回注過(guò)程中，往往需要從油田周?chē)貐^(qū)抽取地下水，這又導(dǎo)致油田周?chē)貐^(qū)進(jìn)一步發(fā)生地面沉降。沉降中心C的孤島-孤東油田便是勝利油田的重要采油區(qū)之一，由圖6(a)可看出，油田中心區(qū)域的沉降量級(jí)較大，地面沉降影響范圍廣，周?chē)貐^(qū)也發(fā)生了緩慢的地表沉降現(xiàn)象。由此可見(jiàn)，黃河三角洲地區(qū)的地表沉降現(xiàn)象與油氣開(kāi)采活動(dòng)也有一定關(guān)系。

5 結(jié) 語(yǔ)

本文采用DS-InSAR方法對(duì)2019-12～2020-12期間26景Sentinel-1A影像進(jìn)行數(shù)據(jù)處理及結(jié)果分析，得到黃河三角洲區(qū)域的地表沉降速率，結(jié)論如下：

1)與傳統(tǒng)PS-InSAR方法相比較，融合分布式目標(biāo)的DS-InSAR方法有效提高了地面監(jiān)測(cè)點(diǎn)的密度，能夠詳細(xì)描述沉降區(qū)域特征，實(shí)現(xiàn)了地表形變的精細(xì)化監(jiān)測(cè)，并且該方法具有很好的可靠性。

2)在監(jiān)測(cè)時(shí)間段內(nèi)，監(jiān)測(cè)到黃河三角洲存在多處不同程度的地表沉降現(xiàn)象，最大沉降速率達(dá)-238 mm/a。同時(shí)利用DS-InSAR方法探測(cè)到黃河三角洲地區(qū)較為準(zhǔn)確的時(shí)序地表形變趨勢(shì)，可為地表形變規(guī)律及地表沉降誘發(fā)因素的研究與分析提供參考。

3)導(dǎo)致黃河三角洲地表形變的因素主要有過(guò)度開(kāi)采地下鹵水、超負(fù)荷開(kāi)采油氣資源等。由于缺乏實(shí)測(cè)數(shù)據(jù)及水文地質(zhì)資料等，本文對(duì)黃河三角洲地表形變影響因素的分析稍有欠缺，之后可根據(jù)實(shí)地調(diào)研情況結(jié)合相關(guān)資料對(duì)地表形變機(jī)理作進(jìn)一步分析與完善。