利用DS-InSAR技術監測黃河三角洲地表形變

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

1 中國礦業大學自然資源部國土環境與災害監測重點實驗室，江蘇省徐州市大學路1號，221116

2 山東省國土測繪院，濟南市經十路2301號，250013

黃河流域是我國重要的生態屏障和經濟地帶，在推動我國社會經濟發展與生態建設方面都具有非常重要的作用。黃河三角洲作為我國最大的三角洲，蘊含著豐富的鹵水、石油、天然氣等自然資源。隨著人們對地下鹵水的抽取及石油的不斷開采，黃河三角洲存在多處不同程度的地表沉降區，嚴重影響當地居民的生產生活。因此，準確監測黃河三角洲地區地表沉降對于減少該地區的自然災害和經濟損失具有重要意義[1]。

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

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

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

1 研究區域及數據概況

1.1 區域概況

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

本文研究區覆蓋河口區、利津縣及墾利縣的大部分區域，包含勝利油田(中國第二大油田)、孤東油田、東營港口、東營機場及大量鹽田區。由于黃河三角洲地區特殊的地質構造特點及大量開采油田、地下鹵水等人類活動導致地層壓力發生變化，該地區存在多個沉降區，且沉降量級大，嚴重影響了當地人民的生產生活。

1.2 實驗數據

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

2 DS-InSAR方法原理

2.1 FaSHPS同質點識別

FaSHPS(fast SHP selection)算法[6]的核心思想是將假設檢驗問題轉化為置信區間估計，通過影像像元數據在時間序列上與參考點的相似性來判斷是否為同質像元點。本文利用FaSHPS方法逐個計算出像元周圍的同質點數量，以20個同質點作為閾值，將同質點數量高于該閾值的參考像元作為DS預選點。

(1)

式中，P{·}表示概率，z1-α/2表示在標準正態分布中1-α/2分位點，μ(S)為像元S的期望，var(A(S))表示像元S在時間序列上振幅的真實方差，N表示SAR影像的個數。

在同質區域可認為SAR影像的振幅服從瑞利分布[10]，則變異系數CV可由下式計算：

(2)

式中，E(·)表示期望。假定像元的散射特性在觀測時間范圍內基本穩定，則式(1)可轉換為：

(3)

2.2 特征值分解相位優化

分布式目標像元內地物應具有相同的散射特性，由于受到噪聲及時空失相干的影響，時間序列SAR影像干涉圖每個像元的信號矢量中包含多種類型的散射信號，使得分布式目標像元的相位穩定性較差。通過對像元的相位優化可達到對干涉圖像元相位降噪的目的。像元相干矩陣對應的特征值代表不同的后向散射特性，將相干矩陣進行特征值分解[5]即可得到最大特征值及其對應的特征向量，認為此特征向量對應的像元散射特性是穩定的，并將其作為優化后的相位。

假設在勻質區域Ω內部存在NP個后向散射特性相近的像元，其相干矩陣T可表示為：

(4)

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

(5)

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

T=Tsignal+Tnoise=

(6)

相位優化完成后，還需對相位的優化質量進行評估。將優化前后的干涉相位進行差值擬合計算得到擬合度[4]：

(7)

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

2.3 處理步驟

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

3 InSAR監測結果分析

3.1 地表形變監測結果

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

3.2 監測結果分析

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

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

A沉降區中心位于東營市河口區廣河村國星鹽場[12]，其最大年沉降速率為-134 mm/a。由于存在過度開采地下鹵水的情況，導致鹽場附近出現較大范圍的沉降漏斗，圖4(a)顯示沿A′-A″線距A′點約5～12 km的范圍內，地表沉降速率達到-75 mm/a。

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

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

3.3 DS-InSAR可靠性分析

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

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

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

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

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

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

5 結 語

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

1)與傳統PS-InSAR方法相比較，融合分布式目標的DS-InSAR方法有效提高了地面監測點的密度，能夠詳細描述沉降區域特征，實現了地表形變的精細化監測，并且該方法具有很好的可靠性。

2)在監測時間段內，監測到黃河三角洲存在多處不同程度的地表沉降現象，最大沉降速率達-238 mm/a。同時利用DS-InSAR方法探測到黃河三角洲地區較為準確的時序地表形變趨勢，可為地表形變規律及地表沉降誘發因素的研究與分析提供參考。

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