夏斯雨，陸藩藩，晏王波

(1.江蘇省天地圖地理信息工程技術公司，江蘇 南京 210000; 2.江蘇省測繪研究所，江蘇 南京 210000)

0 引 言

地面沉降又被稱為地面下沉或地陷，是在自然條件和人為因素作用下，因地殼表層土壓縮而導致區(qū)域性地面標高降低的一種環(huán)境地質現(xiàn)象[1-2]。

為了防治地面沉降帶來的災害，多種技術被用于監(jiān)測地表形變，如傳統(tǒng)的大地測量(GPS和精密水準)，合成孔徑雷達干涉測量(InSAR)技術[3-4]等。常規(guī)的差分干涉測量(D-InSAR)技術在實際運用中容易受空間、時間失相干及大氣延遲的制約，該方法主要用于大尺度地表突變的監(jiān)測，而對于長時間的緩慢的地表形變過程監(jiān)測受到很大的限制。近年來發(fā)展的永久散射體干涉測量(PS-InSAR)技術和短基線(SBAS-InSAR)[6]技術，在城市地面沉降監(jiān)測領域逐步取代了傳統(tǒng)的D-InSAR技術，成為利用SAR數(shù)據(jù)獲取地表形變的主要手段，能有效提高干涉圖的相干性，獲取地表形變的沉降規(guī)律，大幅提高監(jiān)測的應用場景和效益[10]。

淮安市地處江蘇省北部中心地域，北接連云港市，東毗鹽城市，南連揚州市和安徽省滁州市，西鄰宿遷市。東西最大直線距離132 公里，南北最大直線距離150 公里，面積10 072平方公里。項目獲取了研究區(qū)淮安市的數(shù)據(jù)，選用加拿大約25米分辨率的2012-2015年25景RADARSAT-2雷達影像作為數(shù)據(jù)源。選取2012年2月獲取的影像為主影像，其余24景影像為輔影像進行精密配準，計算所有像對的垂直基線。DEM采用的是公開的30 m分辨率DEM數(shù)據(jù)。

1 SBAS-InSAR法求解沉降量原理

Berardino和Lanari等在2002年提出SBAS方法，即通過設定一定的時間、空間閾值自由組合干涉對，削弱時間、空間失相關的影響，該方法是InSAR監(jiān)測中常見的方法之一。在SAR影像數(shù)量一定的情況下，可以形成更多的干涉像對，提高相位觀測量，增加多余觀測量，有利于提高最終解算結果的精度。

假設存在覆蓋同一地區(qū)的N+1幅SAR影像圖，通過設定空間基線與時間基線的閾值，組合干涉對，從干涉圖中去除參考橢球面和地形的相位貢獻，得到M幅差分干涉圖， M滿足下列條件：

(N+1)/2≤M≤(N*(N+1))/2

(1)

對于第j幅干涉圖上的一點，其干涉相位為：

δΦj=Φ(tB)-Φ(tA)≈4π[d(tB)-d(tA)]/λ

(2)

式(2)中：λ為雷達波波長；Φ(tA)和Φ(tB)分別為tA和tB時刻對應的相位；d(tA)和d(tB)是相對于基準時刻t0的雷達視線向的地表形變量。

假設IE=[IE1,……,EM]和IS=[IS1,……,ISM]分別為干涉數(shù)據(jù)處理時根據(jù)時間順序排列的主影像時間序列和從影像時間序列，且滿足IEj>ISj。則所有的差分干涉相位可以組成如下觀測方程：

δΦj=Φ(tIEj)-Φ(tISj)，j=1,……K

(3)

上式是由N個未知數(shù)的M個方程式組成的方程式，即：

δΦj=AΦ

(4)

式(4)中，A是一個M×N的系數(shù)矩陣，如果δΦj=Φ(t4)-Φ(t2)，δΦ2=Φ(t3)-Φ(t2)，則A表示如下：

當影像分多組時，上式是秩虧的，方程組沒有唯一解。此時采用奇異值分解法對矩陣A進行分解，求解出方程的最小范數(shù)意義的最小二乘解。

采用奇異值分解方法直接對相位進行求解，得到的形變結果在時間上出現(xiàn)跳躍，即不連續(xù)，此時沒有物理意義。鑒于此，將未知數(shù)轉變?yōu)橄噜徲跋瘾@取時間內像元點沿LOS向的平均速率，即：

(5)

帶入式δΦj=AΦ得：

(6)

Bv=δΦ

(7)

B是一個M×N的矩陣，對于第j行，位于主從影像之間的列B(j,i)=(ti-ti-1)，其他情況下B(j,i)=0。此時，對B進行奇異值分解，求取最小范數(shù)意義上的LOS方向的相位平均速率，當知道地表形變的先驗信息時，可以利用形變相位在時間域上滿足某些約束條件(即形變模型)來簡化相位平均速率計算的復雜度。

考慮到高程誤差對相位的影響，可以使用沉降模型來聯(lián)合估算DEM誤差，再根據(jù)各時段的沉降速率，對它們進行時間域上的積分，就可以得到各時間段的形變量。

2 淮安地區(qū)沉降速率計算與驗證

2.1 淮安地區(qū)InSAR數(shù)據(jù)處理

SBAS方法的處理步驟主要包括(圖1)：

①首先將監(jiān)測區(qū)的SAR影像配準到同一的坐標空間下，并生成干涉圖集；②通過濾波的方式去除由于失相干引起的噪聲等；③利用輔助的DEM數(shù)據(jù)模擬相位去除地形和平地效應；④利用干涉圖進行相位解纏并進行大氣相位的去除工作；⑤得出淮安監(jiān)測區(qū)沉降的平均位移速率和位置時間序列成果。

圖1 SBAS方法處理步驟

為了保證解纏的效果，確定設定研究區(qū)內相干性閾值為0.35。在選取點的基礎上，根據(jù)相干點構建Delaunay三角網(wǎng)，以三角網(wǎng)的邊為單位進行形變建模，采用最小費用流量法進行差分干涉的解纏，將解纏相位校正到相干性較高的參考點上，求解得出2012—2015年內的灌河口周邊地面沉降的平均形變速率(圖2)。

圖2 淮安市年平均沉降速率圖(2012年-2015年)

通過在時間域上對各時間段的形變速率進行積分即得到各個時間點的累積形變量，并最終得到時間序列形變?？梢钥闯觯夯窗彩袧i水縣存在多個沉降漏斗且已經連接成片，最為嚴重的地區(qū)位于朱碼鎮(zhèn)，最大沉降速率達-42 mm/a，漣城鎮(zhèn)、陳師鎮(zhèn)也有較大沉降。

2.2 淮安地區(qū)CORS數(shù)據(jù)驗證

為驗證SBAS方法得到的結果，搜集了12個灌河口周邊2012-2015年逐年觀測的CORS數(shù)據(jù)，比較CORS實測點與最鄰近點InSAR點目標測量得到的形變值，根據(jù)兩者的差異評估InSAR地表沉降監(jiān)測結果的精度，具體的CORS點分布圖如圖3所示。

由于獲取InSAR監(jiān)測點覆蓋淮安市，CORS站點也分布在研究區(qū)范圍內，所以在驗證時，認為CORS站點與最鄰近InSAR點的地質條件及受人類活動影響相同，故不做額外的誤差分析。

將實測結果與SBAS方法所得監(jiān)測結果進行比較(表1)：

圖3 CORS點分布圖

表1 CORS點與InSAR監(jiān)測結果對比mm/a

二者插值最大值為1.85 mm/a、最小值為-1.34 mm/a。通過SBAS-InSAR技術獲取的地表形變量與已有的CORS點檢測的沉降值進行比較，利用均方根誤差(Root-Mean-Square Error，RMSE)作為評定精度的指標，其公式如下：

(8)

其中，N為站點的個數(shù)；ΔHi為CORS點監(jiān)測到第i個點的沉降結果；Δhi為InSAR檢測到第i個CORS點的沉降結果。計算得到RMSE為3.04 mm/a。

3 結 語

本文利用短基線方法(SBAS)監(jiān)測了淮安市2012-2015年的地表形變，獲取了時間形變的特征、累積沉降量和年平均沉降速率。研究區(qū)內年平均沉降速率最大達-42 mm/a。

(1)從時間角度來看，在監(jiān)測周期內，2012-2015年，沉降漏斗逐漸形成且有成區(qū)連片趨勢，主要原因是各化工園區(qū)企業(yè)規(guī)模數(shù)量不斷增加，生產效能提高，導致地面沉降，故地面沉降逐年累加。

(2)從空間角度來看，沉降區(qū)域相對較為集中，主要分布在漣水縣朱碼鎮(zhèn)、漣城鎮(zhèn)、陳師鎮(zhèn)，大規(guī)模的人類生活生產活動導致了地面的沉降塌陷。

(3)從實驗方法來看，在以平原地貌為主的灌河口周邊利用短基線方法(SBAS)進行地面沉降的監(jiān)測依然有效，這也為大規(guī)模進行全省地面沉降監(jiān)測奠定了基礎，有益于推進全省減災防災工作。