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

(1.江蘇省地質(zhì)調(diào)查研究院，南京 210018；2.地裂縫地質(zhì)災(zāi)害重點實驗室，南京 210018)

1 引 言

蘭西縣位于黑龍江省中西部，松嫩平原東部，隸屬于黑龍江省綏化市。全縣共轄18個鄉(xiāng)、鎮(zhèn)，總?cè)丝?6.1萬人。地方工業(yè)有雄厚的基礎(chǔ)，形成了食品、紡織、機械、服裝、醫(yī)藥、造紙、化工等十個行業(yè)為主體的工業(yè)格局，尤其是亞麻系列產(chǎn)品的開發(fā)已走到了全國的前列。但蘭西縣是否存地下水開采等原因?qū)е碌牡孛娉两瞪胁坏枚?，為了獲取蘭西縣城區(qū)地表形變幅度和分布、采取合理措施控制地面沉降，本文采用InSAR技術(shù)對2012年蘭西縣城區(qū)的地面沉降的狀況進行了監(jiān)測。

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

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

2 數(shù)據(jù)選取與計算方法

2.1 數(shù)據(jù)的選取

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

表1 SAR圖像獲取信息表

圖1 SBAS干涉圖組合方式

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

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

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

圖2 初選的相干點分布圖

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

Dv=δφ

(1)

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

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

v=D+δφ

(2)

D+為D的偽逆矩陣：

(3)

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

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

Dv+C·Δq=δφ

(4)

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

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

3 蘭西縣城區(qū)地面沉降分析

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

圖3 研究區(qū)地表形變散點圖

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

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

4 相鄰軌道重復(fù)區(qū)的精度驗證

圖5 相鄰軌道重復(fù)區(qū)地面沉降散點圖

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

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

5 結(jié)束語

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

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

參考文獻：

[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，F(xiàn)ORNARO 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，F(xiàn)OMARO 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技術(shù)監(jiān)測太原市2003～2009年地表形變場[J].地球物理學(xué)報，2011，54(3)：673-680.

[12] 尹宏杰，朱建軍，李志偉，等.基于SBAS 的礦區(qū)形變監(jiān)測研究[J].測繪學(xué)報，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，F(xiàn)rance[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.