基于SBAS的蘭西縣城區地面沉降監測研究

朱葉飛，陳火根，李向前，詹雅婷，蘇一鳴

(1.江蘇省地質調查研究院，南京 210018；2.地裂縫地質災害重點實驗室，南京 210018)

1 引 言

蘭西縣位于黑龍江省中西部，松嫩平原東部，隸屬于黑龍江省綏化市。全縣共轄18個鄉、鎮，總人口46.1萬人。地方工業有雄厚的基礎，形成了食品、紡織、機械、服裝、醫藥、造紙、化工等十個行業為主體的工業格局，尤其是亞麻系列產品的開發已走到了全國的前列。但蘭西縣是否存地下水開采等原因導致的地面沉降尚不得而知，為了獲取蘭西縣城區地表形變幅度和分布、采取合理措施控制地面沉降，本文采用InSAR技術對2012年蘭西縣城區的地面沉降的狀況進行了監測。

合成孔徑雷達差分干涉測量(DInSAR)是近年來發展起來的一種監測地表形變的手段[1-4]，能夠對地表進行全天候、大范圍、厘米甚至毫米級精度的監測。然而對于長時間緩慢地面沉降而言，由于SAR干涉像對的時間和空間基線距較大，導致影像時間、空間失相干嚴重[5-6]，利用傳統DInSAR技術難以反演其形變序列。近年來國際上一些學者提出小基線集(SBAS)方法，并先后開展了基于SBAS的DInSAR監測地表沉降的應用和研究[7-12]，通過采用SBAS組合進行干涉測量，有效地減弱空間基線引起的失相干問題，提高了DInSAR技術形變計算的精度。SBAS方法利用奇異值分解(SVD)方法將多個小基線集聯合起來求解，有效地解決不同SAR數據集之間空間基線過長造成的時間不連續問題，提高監測的時間分辨率。

本文引入SBAS方法，采用C波段的Radarsat-2數據對蘭西縣城區2012年2月～11月的地面沉降進行監測，獲取該地區時間序列形變場，發現和定位沉降漏斗，對研究區沉降漏斗的時空狀況進行分析，揭示沉降漏斗隨時間演化的規律。

2 數據選取與計算方法

2.1 數據的選取

本文共收集從2012年2月～2012年11月覆蓋研究區的11景Radarsat-2衛星寬模式的單視數據，各圖像具體參數信息如表1所示。該地區的空間垂直基線和時間基線分布較為理想，時間基線最大間隔288天，最大的垂直基線為483m。本文要根據SBAS的原則計算出小基線集劃分的情況，并估算出比較合理的小基線集組合方式(圖1)，共選擇了53個干涉組合，每個組合的垂直基線都小于400m。

表1 SAR圖像獲取信息表

圖1 SBAS干涉圖組合方式

2.2 基于相干點的SBAS計算方法

由于SBAS方法相干目標的提取是基于各個時期的單視數據，相干目標選擇前需要對單視數據進行配準一樣，本文選用2012年6月28日獲取的SAR影像作為參考主影像，將所有影像以復相干系數作為測度配準至主影像，配準的精度控制在1/10個像元內。

在利用SBAS時序分析模型解算形變速率之前，首先要進行相干目標的篩選。目前相干目標像元的篩選方法主要有：基于像元強度的穩定性[13-14]、基于相位穩定性[15]、基于空間相干性[16]、點目標檢測[17-18]等。像元強度穩定性檢測算法是根據像元強度的穩定性來判斷是否為相干點，一般要求SAR影像不少于30景。相位穩定性方法算法復雜，需采用大量模擬試驗來確定閾值。結合本文數據量小的特點并考慮到算法的復雜性，采用空間相干性閾值和點目標檢測的方法來選取相干目標點，這樣既提高相干目標識別準確性，也增大相干目標的數量。本文共選取研究區初始相干點6365個參與運算，初選的相干點分布如圖2所示。

圖2 初選的相干點分布圖

根據SBAS組合原則對收集到的11幅Radarsat-2圖像進行組合，共選擇了53個干涉對，利用研究區工作區DEM去除地形效應后得到差分干涉圖。依據相干點的差分相位δφ可建立矩陣方程：

Dv=δφ

(1)

假設由N+1幅SAR圖像得到M個干涉圖，D是一個M×N矩陣，對第j行，位于主輔圖像獲取時間之間的列，D(j，k)=tk+1-tk，其他的D(j，k)=0；υ[N×1]為相干點的平均相位速度。

將SVD分解應用于矩陣方程(1)，就可以得到速度矢量v的最小范數解。地表形變最小范數意義上的最小二乘相位估計可表示為：

v=D+δφ

(2)

D+為D的偽逆矩陣：

(3)

式中ui為DDT的特征向量；vi為DTD的特征向量；λi為DDT的特征值。

從差分相位的組成出發，除了形變相位分量外，還有高程誤差Δq的相位分量，建立方程組如下：

Dv+C·Δq=δφ

(4)

其中，C[MX1]是與基線距相關的系數矩陣，由此可以得到DEM誤差。

由于大氣相位在空間上表現為低頻、在時間上表現為高頻，因此對以上線性模型的殘余相位在空間和時間上進行適當濾波，分離出大氣相位和非線性形變相位。將非線性形變量與線性形變量相加，從而得到各個時間段累計的沉降量。

3 蘭西縣城區地面沉降分析

將研究區域的SAR數據采用SBAS方法進行時序分析得到蘭西縣城區2012年2月5日～2012年11月19日的高相干點的平均沉降速率(圖3)。通過圖4可以看出，2012年蘭西縣城區存在一定范圍的地面沉降，主要呈漏斗形分布，另外在城區東側存在局部小范圍的沉降。沉降中心位于蘭西縣城區南側5號點附近，最大地面沉降速率達到2.43cm/年，以此為中心地面沉降幅度向四周逐漸趨緩。筆者對研究區進行實地野外調查，發現沉降主要原因為地下水的開采導致的地面沉降，包括各類企業(如食品廠等)利用深水井進行地下水開采、開發區大量新建高層建筑基坑排水、城鎮居民生活用地下水開采等，在局部地區表現為墻體開裂等沉降跡象。

圖3 研究區地表形變散點圖

在反演線性形變速率的基礎上，通過空間濾波和時間濾波相結合的方式去除大氣延遲相位的影響，反演得到非線性形變量。將線性形變量和非線性形變量結果疊加，得到研究區各時間段的累計形變量，選取散點圖中沉降漏斗南北向剖面線上7個典型相干點，各相干點的地面沉降累計沉降序列如圖4所示。從時間序列上看雖然不同點的沉降曲線有所差異，但蘭西縣城區沉降漏斗的地面沉降基本上是呈現線性的沉降。

圖4 研究區典型高相干點累計沉降序列圖

4 相鄰軌道重復區的精度驗證

圖5 相鄰軌道重復區地面沉降散點圖

由于沒有研究區同期的水準觀測數據對InSAR結果進行驗證，本文選取東側相鄰軌道的9景Radarsat-2數據，該幅圖完全覆蓋研究區范圍。對相鄰軌道的研究區進行SBAS方法的獨立計算，并對重復研究區的InSAR計算結果進行相互對比驗證。選取的9景SAR數據獲取時間在2012年1月29日～2012年10月19日之間，依據SBAS組合原則共選擇了28個干涉組合，計算得到的2012年線性形變速率如圖5所示。在重復區域共選取研究區7個典型的同名高相干點，通過比較同名點的地表形變速率(表2)發現：兩個相鄰軌道獨立計算的同名點2012年地表形變線性速率的相對誤差都在0.4cm/年以內，誤差的均方差為0.132cm/年。由于不同軌道圖像獲取時間的差異、軌道誤差及大氣延遲相位的差異，則計算的結果存在一定的差異，但上述精度驗證表明相鄰軌道重復區上的計算結果非常一致，表明本文的計算結果準確地反映了研究區地面沉降的幅度和分布狀況。

表2 高相干同名點的相互驗證(單位：cm/年)

5 結束語

本文引入SBAS方法監測蘭西縣城區地面沉降，突破InSAR觀測周期以及空間和時間基線的限制，對DInSAR單次離散的觀測進行時間序列分析，準確地獲取研究區的地表形變速率圖，發現了該研究區沉降漏斗中心速率達到2.43cm/年，并獲得了該研究區沉降漏斗地表形變累積量的時序變化，發現該漏斗區地面沉降基本呈現線性沉降的規律。經實地調查發現研究區沉降的原因主要是地下水的過度開采。傳統的以“點”為基礎的研究區水準、GPS測量，由于空間采樣率低，很難準確發現和定位沉降漏斗的中心、范圍、大小和演化規律，而本文研究結果對預防蘭西縣城區地質災害和環境破壞具有重要的意義。

在沒有同期的水準測量數據的情況下，通過選用相鄰軌道SAR數據進行InSAR獨立計算并進行精度的對比驗證，發現兩者的計算結果非常一致，表明本文計算結果準確地反映了研究區的地表形變的狀況。下一步工作是獲取更多的SAR圖像參與SBAS的運算以提高InSAR計算的精度，并開展外業水準測量工作與InSAR計算的結果進行對比。

參考文獻：

[1] GABRIEL A K，GOLDSTEIN R M，ZEBKER H A.Mapping small elevation changes over large areas：Differential radar interferometry[J].Journal of Geophysical Research，1989，94(B7)：9183-9191.

[2] MASSONNET D，ROSSI M，CARMONA C，et al.The displacement field of the landers earthquake mapped by radar interferometry[J].Nature，1993，364(8)：138-142.

[3] DING X L，LIU G X，LI Z W，et al.Ground subsidence monitoring in Hong Kong with satellite SAR interferometry[J].Photogrammetry Engineering and Remote Sensing，2004，70(10)：1151-1156.

[4] HU J，LI Z W，DING X L，et al.Two-dimensional co-seismic surface displacements field of the chi-chi earthquake inferred from SAR image matching[J].Sensors，2008，8(10)：6484-6495.

[5] GATELLI F，GUAMIERI A M，PARIZZI F，et al.The wavenumber shift in SAR interfereometry[J].IEEE Transactions on Geosciences and Remote Sensing，1994，32(4)：855-864.

[6] ZEBKER H A，ROSEN P A，HENSLEY S.Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps[J].Journal of Geophysical Research，1997，102(B4)：7547-7563.

[7] BERARDINO P，FORNARO G，LANARI R，et al.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].IEEE Transcations on Geoscience and Remote Sensing，2002，40 (11)：2375-2383.

[8] LANARI R，MORA O，MANUNTA M，et al.A small baseline approach for investigating deformations on full resolution differential SAR interferograms[J].IEEE Transactions on Geoscience and Remote Sensing，2004，42(7)：1377-1386.

[9] CASU F，MANZO M，LANARI R.A quantitative assessment of the SBAS algorithm performance for surface deformation retrieval from DInSAR data[J].Remote Sensing of Environment，2006，102：195-210.

[10] BERARDINO P，CASU F，FOMARO G，et al.Small baseline DIFSAR techniques for earth surface deformation analysis[C].Istituto Per Il Rilevamento Elettromagnetico Dell’Ambiente，IREA-National Research Council of Italy (CNR)，Napoli：2003.

[11] 吳宏安，張永紅，陳曉勇，等.基于小基線DInSAR技術監測太原市2003～2009年地表形變場[J].地球物理學報，2011，54(3)：673-680.

[12] 尹宏杰，朱建軍，李志偉，等.基于SBAS 的礦區形變監測研究[J].測繪學報，2011，40(1)：52-58.

[13] FERRETTI A，PRATI C，ROCCA F.Nonlinear subsidence rate estimation using permanent scatters in differential SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing，2000，38(5)：2202-2212.

[14] KAMPES B M.Displacement parameter estimation using permanent scatterer interferometry[D].Delft：Delft University of Technology，2005.

[15] HOOPER A，ZEBKER H A，SEGALL P，et al.A new method for measuring feformation on volcanoes and other natural terrains using InSAR persistent scatterers[J].Geophysical Research Letters，2004，31 (23)：1-5.

[16] MORA O，MALLORQUI J J，BROQUETAS A.Linear and nonlinear terrain deformation maps from a reduced det of interferometric SAR images[J].IEEE Transactions on Geoscience and Remote Sensing，2003，41(10)：2243-2253.

[17] WNMER C，WEGMULLER U，STROZZI T，et al.Interferometric point target analysis for deformation mapping[A].IGARSS’03，Toulouse，France[C].21-25，July，2003，7：4362-4364.

[18] ANDREAS W，WEMER C，STROZZI T，et al.Combination of point and extended target based interferometric techniques[A].IGARSS’04，Anchorage，Alaska，USA[C].20-24，Sep，2004，2：989-991.