曹建濤 鄭翔元 范洪冬 李國華 黃 晨

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

2 山東省國土測繪院，濟(jì)南市經(jīng)十路2301號，250013

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

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

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

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

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

1.1 區(qū)域概況

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

本文研究區(qū)覆蓋河口區(qū)、利津縣及墾利縣的大部分區(qū)域，包含勝利油田(中國第二大油田)、孤東油田、東營港口、東營機(jī)場及大量鹽田區(qū)。由于黃河三角洲地區(qū)特殊的地質(zhì)構(gòu)造特點(diǎn)及大量開采油田、地下鹵水等人類活動導(dǎo)致地層壓力發(fā)生變化，該地區(qū)存在多個(gè)沉降區(qū)，且沉降量級大，嚴(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è)干涉對。

2 DS-InSAR方法原理

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

FaSHPS(fast SHP selection)算法[6]的核心思想是將假設(shè)檢驗(yàn)問題轉(zhuǎn)化為置信區(qū)間估計(jì)，通過影像像元數(shù)據(jù)在時(shí)間序列上與參考點(diǎn)的相似性來判斷是否為同質(zhì)像元點(diǎn)。本文利用FaSHPS方法逐個(gè)計(jì)算出像元周圍的同質(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(·)表示期望。假定像元的散射特性在觀測時(shí)間范圍內(nèi)基本穩(wěn)定，則式(1)可轉(zhuǎn)換為：

(3)

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

分布式目標(biāo)像元內(nèi)地物應(yīng)具有相同的散射特性，由于受到噪聲及時(shí)空失相干的影響，時(shí)間序列SAR影像干涉圖每個(gè)像元的信號矢量中包含多種類型的散射信號，使得分布式目標(biāo)像元的相位穩(wěn)定性較差。通過對像元的相位優(yōu)化可達(dá)到對干涉圖像元相位降噪的目的。像元相干矩陣對應(yīng)的特征值代表不同的后向散射特性，將相干矩陣進(jìn)行特征值分解[5]即可得到最大特征值及其對應(yīng)的特征向量，認(rèn)為此特征向量對應(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ù)觀測量經(jīng)歸一化處理后的復(fù)數(shù)向量，(·)H表示矩陣的共軛轉(zhuǎn)置。上式得到的相干矩陣為半正定Hermitian矩陣，經(jīng)特征值分解可得：

(5)

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

T=Tsignal+Tnoise=

(6)

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

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

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

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

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

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

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

A沉降區(qū)中心位于東營市河口區(qū)廣河村國星鹽場[12]，其最大年沉降速率為-134 mm/a。由于存在過度開采地下鹵水的情況，導(dǎo)致鹽場附近出現(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ū)域被黃河故道、東營港疏港高速及孤北水庫包圍，自2015年就已出現(xiàn)沉降現(xiàn)象，由于鹽田分布較為集中且開采力度大，如今地表沉降速率明顯加快。

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

3.3 DS-InSAR可靠性分析

為進(jìn)一步驗(yàn)證DS-InSAR方法的可靠性，本文通過對比分析PS-InSAR與DS-InSAR在一定范圍內(nèi)同名監(jiān)測點(diǎn)的沉降速率進(jìn)行交叉驗(yàn)證。將PS-InSAR監(jiān)測結(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)測得到的地表沉降速率具有良好的一致性，其同名點(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ū)、墾利縣、東營區(qū)已建成多處鹽場并投入生產(chǎn)，這些工廠開采地下鹵水用于工業(yè)制溴、制鹽和制造相關(guān)化工產(chǎn)品等，但長期不合理開采地下鹵水將會導(dǎo)致地下鹵水水位下降、局部海水入侵和地面裂縫，嚴(yán)重時(shí)則會導(dǎo)致地表沉降、發(fā)生地質(zhì)災(zāi)害。眾多企業(yè)大量開采地下鹵水不僅致使其資源減少，而且形成淺層地下鹵水沉降漏斗。因此，過度開采地下鹵水用于制鹽和化工產(chǎn)業(yè)是造成黃河三角洲地區(qū)地面沉降的主要原因之一。

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

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

5 結(jié) 語

本文采用DS-InSAR方法對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)測點(diǎn)的密度，能夠詳細(xì)描述沉降區(qū)域特征，實(shí)現(xiàn)了地表形變的精細(xì)化監(jiān)測，并且該方法具有很好的可靠性。

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

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