短基線集形變模型反演的正則化解算方法

姜兆英，劉國林，陶秋香

1. 山東科技大學(xué)測繪科學(xué)與工程學(xué)院，山東 青島 266590; 2. 青島農(nóng)業(yè)大學(xué)理學(xué)與信息科學(xué)學(xué)院，山東 青島 266109

?

短基線集形變模型反演的正則化解算方法

姜兆英1,2，劉國林1，陶秋香1

1. 山東科技大學(xué)測繪科學(xué)與工程學(xué)院，山東 青島 266590; 2. 青島農(nóng)業(yè)大學(xué)理學(xué)與信息科學(xué)學(xué)院，山東 青島 266109

Foundation support： The National Natural Science Foundation of China (Nos. 41274007; 41404003); Shandong Province Natural Science Foundation of China (No. ZR2012DM001); Specialized Research Fund for the Doctoral Program of Higher Education (No. 20123718110001); Shandong Taishan Scholar Construction Project under Special Funding (No. TSXZ201509)

針對短基線集形變模型反演中法方程系數(shù)矩陣呈病態(tài)的問題，提出一種正則化穩(wěn)健解算方法。該方法基于Tikhonov正則化理論，將形變速率求解問題轉(zhuǎn)化為極小化問題，根據(jù)L-曲線法選取正則化參數(shù)，考慮最小二乘殘差各個分量間的關(guān)系選取正則化矩陣，實(shí)現(xiàn)短基線集形變模型反演的穩(wěn)健解算。分別采用LS法、嶺估計(jì)法和Tikhonov正則化法對覆蓋北京地區(qū)的29景ENVISAT ASAR數(shù)據(jù)進(jìn)行處理，反演出研究區(qū)沉降速率圖。通過對代表不同沉降情況的21個點(diǎn)的均方誤差值和時間相干值、整個研究區(qū)的均方誤差圖等的對比分析，表明本文提出的短基線集形變模型反演的正則化穩(wěn)健解算方法可獲取更可靠的形變監(jiān)測結(jié)果。

短基線集；條件數(shù)；Tikhonov正則化；嶺估計(jì)；均方誤差

2002年，文獻(xiàn)[1]提出了短基線集(small baseline subset,SBAS)技術(shù)。該技術(shù)在一定程度上削弱了傳統(tǒng)DInSAR技術(shù)的空間基線失相干和大氣效應(yīng)影響，得到了廣泛的應(yīng)用。文獻(xiàn)[2]將該技術(shù)應(yīng)用于Mono盆地和Long Valley火山形變監(jiān)測；文獻(xiàn)[3]將其應(yīng)用于礦區(qū)形變監(jiān)測；文獻(xiàn)[4]將其應(yīng)用于青藏高原季節(jié)性凍土形變監(jiān)測；文獻(xiàn)[5—7]將其應(yīng)用于城市地表形變監(jiān)測。在上述應(yīng)用中，當(dāng)時序干涉對序列屬于同一集合時，SBAS形變模型的系數(shù)陣是列滿秩的，其法方程系數(shù)陣滿秩，主要采用最小二乘法(least squares,LS)[8-10]解算形變量和形變速率。當(dāng)時序干涉對序列屬于不同集合時，SBAS形變模型的系數(shù)陣不是列滿秩，其法方程系數(shù)陣秩虧，主要采用奇異值分解法(singular value decomposition，SVD)[11-12]進(jìn)行形變量和形變速率的反演。文獻(xiàn)[13]還提出利用加權(quán)最小二乘反演算法進(jìn)行SBAS形變模型求解來提高地表形變監(jiān)測的能力。文獻(xiàn)[14]提出基于L1范數(shù)估計(jì)準(zhǔn)則反演地表形變，通過模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)說明了基于L1范數(shù)估計(jì)準(zhǔn)則所得到的解比基于L2范數(shù)估計(jì)準(zhǔn)則所得到的解更準(zhǔn)確。

對于方程組(1)

AX=b

(1)

式中，A∈Rm×n；b∈Rm為已知量；X∈Rn為待求量。若它的解滿足①存在；②唯一；③穩(wěn)定，則稱方程組(1)是適定的。反之，上述3個條件中至少有一個不滿足，則稱為不適定的。而不適定性是反問題本身固有的一種特性[15]。在SBAS數(shù)據(jù)處理過程中，干涉對的組合方式在一定程度上決定了模型的設(shè)計(jì)矩陣。按照事先設(shè)置的相干系數(shù)、時間基線、空間基線閾值，干涉對自由組合，往往存在某一景影像和多景影像形成干涉像對(如本次試驗(yàn)中第9景數(shù)據(jù)與其他7景數(shù)據(jù)形成干涉對)，這將會導(dǎo)致模型設(shè)計(jì)矩陣的列向量間產(chǎn)生近似線性相關(guān)，使其法方程的系數(shù)陣雖滿秩但卻產(chǎn)生了病態(tài)，即求解短基線集形變模型為不適定問題。求解不適定問題的一般方法是用一組與原不適定問題相“鄰近”的適定問題的解去逼近原問題的解,這種方法稱為正則化方法。如何建立有效的正則化方法是反問題領(lǐng)域中不適定問題研究的重要內(nèi)容。通常的正則化方法有基于變分原理的Tikhonov正則化法[16]、各種迭代方法[17-18]以及其他的一些改進(jìn)方法。

在SBAS處理過程中當(dāng)時序干涉對序列屬于同一個集合，設(shè)計(jì)矩陣是列滿秩，法方程的系數(shù)陣滿秩但卻呈病態(tài)，此時最小二乘估計(jì)的性質(zhì)顯著變壞，形變量和形變速率反演精度將會變差。本文針對這一問題，將Tikhonov正則化方法引入到短基線集形變模型的求解中。基于Tikhonov正則化理論，將形變速率求解問題轉(zhuǎn)化為極小化問題，根據(jù)L-曲線法，選取了合理的正則化參數(shù)；考慮最小二乘殘差的各個分量間的關(guān)系選取了正則化矩陣，實(shí)現(xiàn)了SBAS形變模型的穩(wěn)健解算。在簡要介紹短基線集技術(shù)原理的基礎(chǔ)上，詳細(xì)分析了Tikhonov正則化算法，并給出正則化參數(shù)和正則化矩陣選取的準(zhǔn)則；利用ASAR真實(shí)數(shù)據(jù)進(jìn)行試驗(yàn)，并與LS估計(jì)法、嶺估計(jì)法進(jìn)行了對比分析，結(jié)果表明利用Tikhonov正則化方法反演地表形變，能提高解的穩(wěn)定性和精確度，可獲取更可靠的形變監(jiān)測結(jié)果。

1　短基線集形變模型解算的Tikhonov正則化算法

1.1短基線集形變模型

假設(shè)有N+1幅覆蓋同一區(qū)域且已配準(zhǔn)好的單視復(fù)(SLC)SAR影像，影像獲取時間依次為t0、t1、…、tN。通過限制最大的時間基線和空間基線，生成M個多視的短基線差分干涉圖。對這M個差分干涉圖進(jìn)行相位解纏，然后將所有的像元點(diǎn)參考到一個具有高相干性或沉降量已知的點(diǎn)。假設(shè)第j幅干涉圖是由tA、tB時刻獲取的兩幅SAR影像進(jìn)行干涉生成的，那么其高相干像元的差分相位值可表示為

(2)式中，λ為雷達(dá)波長；φ(tB)和φ(tA)分別為tB、tA時刻相對于t0(設(shè)為參考點(diǎn)，即φ(t0)=0)的形變引起的相位；4πB⊥jΔz/λrsinθ是可能的高程誤差Δz引起的殘余地形相位部分；φatm(tB)-φatm(tA)為獲取影像時由于大氣不一致而引起大氣相位部分；Δnj為包含的失相關(guān)影響及其他噪聲影響。

忽略大氣相位和噪聲相位部分，式(2)用矩陣可表示為

Aφ+CΔz=δφ

(3)

(4)

將式(4)帶入式(3)得

Bv+CΔz=δφ

(5)

式中，B∈RM×N；v∈RN。由文獻(xiàn)[1]所提到的方法，借助相位的三次模型，形成一個具有M個方程(這里δφ為解纏的相位值)，4個未知量的線性系統(tǒng)，借助最小二乘估計(jì)出干涉相位中由DEM誤差引起的相位部分CΔz。從解纏的相位減去估計(jì)的DEM誤差，得到如下線性系統(tǒng)

Bv=δφ′

(6)

式中，δφ′=δφ-CΔz。求出v后，再在時間段上求積分即可得到最后的相位φ。

1.2形變模型反演的Tikhonov正則化算法

對于線性系統(tǒng)(式(6))的典型反演問題，可將其改寫成

Bv=δφ′+n

(7)

式中，n為誤差向量，包含殘留的DEM誤差、大氣不均勻引起的誤差以及噪聲誤差等。

若B是列滿秩的，可用最小二乘法進(jìn)行解算，即使誤差向量的L2范數(shù)最小

(8)

從而可得其法方程為

(9)

由于B是列滿秩的，所以BTB是可逆陣，故解式(9)即得式(7)的最小二乘解為

(10)

最小二乘解在線性無偏估計(jì)類中是唯一的具有最小方差的估計(jì)。正因?yàn)檫@一點(diǎn)，最小二乘估計(jì)在線性系統(tǒng)模型的估計(jì)理論與實(shí)際應(yīng)用中占有非常重要的地位。但是在用最小二乘求解反問題時，其結(jié)果并不總是令人滿意。如文獻(xiàn)[19]的47頁給出的例子中，所求參數(shù)的真值為0，但用LS法求出該參數(shù)的估計(jì)值是9.363 1，說明在病態(tài)情況下，LS法會使某些參數(shù)估計(jì)值遠(yuǎn)遠(yuǎn)偏離真值。研究表明，在有些情況下，LS估計(jì)易受數(shù)據(jù)中的異常值的影響，即系數(shù)矩陣或常數(shù)項(xiàng)的微小變化可能會引起方程組解的巨大變化，這種方程組稱為“病態(tài)”方程組。判斷方程組是否“病態(tài)”或“病態(tài)”的程度，通常采用式(11)所示系數(shù)矩陣的條件數(shù)來衡量

(11)

(12)

式中，λmax、λmin表示矩陣N的最大、最小特征值。

(13)

則有[20]

(14)

式(14)表明方程組解的相對精度不僅取決于觀測值相對精度，還依賴于系數(shù)矩陣的條件數(shù)。當(dāng)矩陣N的條件數(shù)很大時，由式(14)知，觀測值很少的擾動都會引起方程組解很大的誤差。而在短基線集線性系統(tǒng)中，一般所采用的干涉圖組合方式，都會導(dǎo)致法方程的系數(shù)矩陣的條件數(shù)大于1000，此時式(9)是一個病態(tài)方程。在這種情況下，雖然最小二乘估計(jì)的方差在線性無偏估計(jì)類中最小，但其值卻很大，使LS估計(jì)的精度較差，表現(xiàn)出相當(dāng)?shù)牟环€(wěn)定。

為了獲取問題的最優(yōu)穩(wěn)定解，根據(jù)Tikhonov正則化理論，將地面沉降的平均速率求解問題轉(zhuǎn)化為求解下述的極小化問題，即

(15)

(BTB+αR)v=BTδφ′

(16)

2　正則化參數(shù)和矩陣的確定

Tikhonov正則化的效果主要取決于正則化參數(shù)α和正則化矩陣R的選取。正則化參數(shù)α的選擇對問題的解起著關(guān)鍵的作用。如果α太小，則對問題病態(tài)性的改善沒有起到什么作用，即解的不穩(wěn)定性依然存在；反之，如果α太大，則問題(16)雖然可以穩(wěn)定地求解，但該解與原方程(9)的解相差甚遠(yuǎn)，是一個相當(dāng)糟糕的逼近(poorapproximation)。所以最優(yōu)的正則化參數(shù)選取應(yīng)當(dāng)兼顧這兩種情況。通常，參數(shù)α的選取有先驗(yàn)(apriori)和后驗(yàn)(aposteriori)兩種策略。先驗(yàn)選法是在求正則化解之前，根據(jù)正則化參數(shù)α具有的數(shù)學(xué)特性把它確定出來，主要是依據(jù)經(jīng)驗(yàn)進(jìn)行，這種方法多具有理論上的分析意義，不具有實(shí)際應(yīng)用價值。后驗(yàn)方法是指在計(jì)算正則解的過程中根據(jù)一定的數(shù)學(xué)原則確定參數(shù)α。常見的后驗(yàn)正則化參數(shù)選取準(zhǔn)則主要有：①M(fèi)orozov的偏差原理；②Engl準(zhǔn)則；③廣義偏差原理；④Arcangeli準(zhǔn)則與廣義Arcangeli準(zhǔn)則；⑤擬最優(yōu)準(zhǔn)則；⑥廣義交叉效驗(yàn)(GCV)準(zhǔn)則；⑦L-曲線準(zhǔn)則，它們的數(shù)學(xué)理論見文獻(xiàn)[21]。本文以L-曲線準(zhǔn)則來確定正則化參數(shù)α。

令

(17)

圖1　L-曲線圖Fig.1　L-curve image

該曲線的曲率計(jì)算公式為

(18)

正則化矩陣的選取需根據(jù)參數(shù)的物理意義來確定，它的適當(dāng)選取有助于抑制偽目標(biāo)，增強(qiáng)解的可信性，即合適的正則化矩陣可以降低較小的特征值所對應(yīng)的模型參數(shù)中不可靠的信息對結(jié)果的影響。若選取不當(dāng)，會造成參數(shù)的過度平滑。正則化矩陣R改進(jìn)的是設(shè)計(jì)矩陣的相關(guān)性，文獻(xiàn)[22]對于正則化矩陣的選取問題進(jìn)行了試驗(yàn)研究，提出了事后確定方法，在一定程度上反映了參數(shù)間的數(shù)值關(guān)系，參數(shù)估計(jì)效果顯著。在本文的試驗(yàn)處理中，依據(jù)文獻(xiàn)[22]的方法，首先算出最小二乘估計(jì)的均方誤差陣(meansquarederrormatrix，記為MSEM)，求其逆后再取對角線上的元素作為式(17)中的正則化矩陣R。

確定正則化參數(shù)α和正則化矩陣R后，方程(16)的解為

(19)

其相應(yīng)的均方誤差陣[23]為

3　試驗(yàn)結(jié)果及分析

本文所用數(shù)據(jù)是覆蓋北京地區(qū)的29景ENVISAT ASAR降軌數(shù)據(jù)，時間跨度為2007年1月至2010年10月。試驗(yàn)?zāi)康氖菫榱吮容^不同反演算法所得結(jié)果的穩(wěn)定性和可靠性，不討論研究區(qū)的沉降情況。所以在數(shù)據(jù)覆蓋區(qū)域，僅選取中心地理坐標(biāo)為(116.369 8°E、39.914 5°N)面積約為900 km2的一塊區(qū)域。在數(shù)據(jù)處理過程中，采用了荷蘭Delft大學(xué)發(fā)布的DORIS精軌數(shù)據(jù)和SRTM的90 m分辨率的外部DEM數(shù)據(jù)。

首先根據(jù)三基線和最小原則[25]，選取2009年1月26日獲取的SAR影像為主影像，將其他輔影像配準(zhǔn)到主影像系統(tǒng)下。按照選取的最大垂直基線距400 m，最大時間基線距1050 d，最小的相干系數(shù)0.5，進(jìn)行多主影像干涉組合，最終得到55個干涉對，如圖2所示。通過設(shè)置0.6的幅度差分離差閾值，最終選取了69 791個高相干點(diǎn)。圖3為短基線集技術(shù)所得到的試驗(yàn)區(qū)的年平均沉降速率圖。

圖2　干涉圖的組合方式Fig.2　Combinations of interferograms

根據(jù)上述干涉對的組合方式，得到此時對應(yīng)式(9)的系數(shù)矩陣的條件數(shù)為cond(BTB)=3.056 5×103，呈現(xiàn)病態(tài)。從圖3中可以看出，經(jīng)度為116.457 8°的紅色直線跨過了試驗(yàn)區(qū)的未沉降區(qū)、沉降區(qū)和抬升區(qū)，為了更好驗(yàn)證SBAS形變模型反演的正則化解算方法的有效性，在這條直線上選取了如圖3黑色點(diǎn)所示的21個點(diǎn)，這21個點(diǎn)分別代表了試驗(yàn)區(qū)的3種不同沉降情況(未沉降、沉降和抬升)，按照由北向南的順序依次記為1、2、…、20、21。圖4給出了這21個點(diǎn)的年平均沉降速率散點(diǎn)圖。

從表1可見，嶺估計(jì)所得結(jié)果的均方誤差稍微小于最小二乘所得的均方誤差，說明嶺估計(jì)雖然改善了法方程的病態(tài)問題，但是所得結(jié)果的精確度與最小二乘相差不大。而基于Tikhonov正則化估計(jì)法所得結(jié)果的MSE遠(yuǎn)小于LS、嶺估計(jì)所得結(jié)果，圖5直觀地顯示了這3種算法所得結(jié)果的均方誤差的大小變化情況。表1和圖5中3種解算方法均方誤差MSE值的大小表明Tikhonov正則化估計(jì)法可以改善LS估計(jì)和嶺估計(jì)的結(jié)果，提高計(jì)算結(jié)果的精確度。

表1　21個選取點(diǎn)用不同方法所得的)的結(jié)果對比

圖6給出了所選取的21個點(diǎn)在短基線集技術(shù)處理過程中所獲得的相干系數(shù)值的曲線圖，圖7顯示了這21個點(diǎn)基于Tikhonov正則化估計(jì)法所得均方誤差的曲線圖。通過對比圖6和圖7，得到在所選的21個點(diǎn)當(dāng)中，大部分點(diǎn)具有這樣的關(guān)系：相干系數(shù)值大的點(diǎn)對應(yīng)的MSE小，相干系數(shù)值小的點(diǎn)對應(yīng)的MSE大。但在個別點(diǎn)如點(diǎn)6、8和11不符合這一關(guān)系，是因?yàn)椴B(tài)有Ⅰ、Ⅱ類病態(tài)。即使對某一類病態(tài)問題，采用不同處理方法如嶺估計(jì)、正則化等，其結(jié)果也會不完全相同。甚至對同一個問題，采用不同的處理方法，其效果也會不相同。

分析圖4、6、7，發(fā)現(xiàn)點(diǎn)1和點(diǎn)2，點(diǎn)2和點(diǎn)3相干系數(shù)差別不大，但均方誤差相差顯著，這主要是因?yàn)椋孩傧噜忺c(diǎn)1、2、3相對于其他18個點(diǎn)的年平均沉降速率相對較大，分別是-35mm/a、-26mm/a、-15mm/a，這3點(diǎn)均處于低相干區(qū)，且均基于選定窗口內(nèi)的灰度(或幅度)值計(jì)算出的相干系數(shù)，該值的大小會受到窗口內(nèi)鄰近像元灰度值的影響，所以其值相差不大；②在點(diǎn)1、2、3中，點(diǎn)2的相干系數(shù)值最大，說明點(diǎn)2受到除沉降以外的大氣延遲、植被等的影響最小，其去相干噪聲也最小，對噪聲反應(yīng)較敏感的點(diǎn)1、2、3的均方誤差值中，點(diǎn)2的均方誤差值就最小，且與其他兩點(diǎn)(1、3點(diǎn))相差比較顯著。

為了進(jìn)一步對比Tikhonov正則化法與LS法所得結(jié)果，本文引入文獻(xiàn)[26]定義的時間相干值(temporalcoherencevalues)

(21)

表2　所選點(diǎn)用LS法、Tikhonov正則化法計(jì)算所得的時間相干值

雖然對該研究區(qū)、相同時間段，缺少相應(yīng)的GPS、水準(zhǔn)等外部數(shù)據(jù)，但已有學(xué)者[10]對北京地區(qū)采用類似數(shù)據(jù)開展過監(jiān)測。利用與文獻(xiàn)[10]研究結(jié)果對比和內(nèi)部精度驗(yàn)證指標(biāo)(均方誤差)這兩種方式說明在SBAS形變模型反演中，本文提出的Tikhonov正則化方法優(yōu)于常用的LS法。

文獻(xiàn)[10]中，利用兩個相鄰軌道的ENVISATASAR數(shù)據(jù)監(jiān)測北京市地面沉降，得到朝陽區(qū)的三間房鄉(xiāng)、管莊鄉(xiāng)一帶在2007年6月至2010年10月期間的沉降速率達(dá)到了100mm/a。而本文借助軌道490的ENVISATASAR數(shù)據(jù)，由LS法解算得到該地區(qū)在2007年1月至2010年10月期間的最大沉降速率為82mm/a，由Tikhonov正則化方法解算得到的最大沉降速率為90mm/a，這說明在相同數(shù)據(jù)類型和相同景數(shù)的前提下，SBAS處理過程中采用Tikhonov正則化方法解算的結(jié)果比采用LS解算的結(jié)果更準(zhǔn)確。

圖9給出3種方法解算所得到整個研究區(qū)的均方誤差圖。

從圖9可以看出， 圖9(c)的值在對應(yīng)點(diǎn)處遠(yuǎn)小于(a)、(b)的值，說明Tikhonov正則化法所求解的均方誤差遠(yuǎn)小于LS法、嶺估計(jì)法所求的，進(jìn)一步說明Tikhonov正則化法求解精度高于其他兩種方法。

4　總　結(jié)

在短基線集技術(shù)形變模型的反演中，若時序干涉對序列屬于同一集合時，法方程系數(shù)矩陣的條件數(shù)大于103是存在的，從而呈現(xiàn)病態(tài)性。為了克服法矩陣的病態(tài)使估計(jì)值有較大偏差問題，本文以均方誤差為估計(jì)準(zhǔn)則，根據(jù)Tikhonov正則化法將地表形變的平均速率反演問題轉(zhuǎn)化為求目標(biāo)函數(shù)的極小值問題。根據(jù)L-曲線準(zhǔn)則確定正則化參數(shù)α，最小二乘估計(jì)的均方誤差矩陣的逆的對角線上的元素構(gòu)成正則化矩陣R，實(shí)現(xiàn)SBAS形變模型反演的穩(wěn)健解算。分別采用LS法、嶺估計(jì)法和Tikhonov正則化法對覆蓋北京地區(qū)的29景ENVISAT ASAR數(shù)據(jù)進(jìn)行處理，反演出研究區(qū)沉降速率圖。通過對代表不同沉降情況的21個點(diǎn)的均方誤差值和時間相干值、整個研究區(qū)的均方誤差圖等的對比分析，表明基于Tikhonov正則化法估計(jì)的解更穩(wěn)定，具有更高的精確度。

圖3　短基線集技術(shù)得到的研究區(qū)在2007—2010年期間的年平均沉降速率圖Fig.3　Annual mean subsidence rate of the test area is derived by SBAS approach from 2007 to 2010

圖4　21個點(diǎn)的年平均沉降速率散點(diǎn)圖Fig.4　Scatter diagram of annual mean subsidence rate of the selected 21 points

圖5　3種算法所得21個點(diǎn)MSE的曲線圖Fig.5　Mean squares error calculated by LS, ridge and Tikhonov regularization estimation of the 21 selected points

圖6　21個點(diǎn)的相干系數(shù)值Fig.6　Coherence image of the selected 21 points

圖7　21個點(diǎn)的基于Tikhonov正則化估計(jì)法所得均方誤差圖Fig.7　Mean squares errors of the selected 21 points based on Tikhonov regularization estimation

圖8　所選點(diǎn)的時間相干值散點(diǎn)圖Fig.8　Temporal coherence values of selected 21 points

圖9　3種方法所得結(jié)果的均方誤差Fig.9　The mean square error calculated by the three methods

[1]BERARDINO P, FORNARO G, LANARI R, et al. A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(11): 2375-2383.

[2]TIZZANI P, BERARDINO P, CASU F, et al. Surface Deformation of Long Valley Caldera and Mono Basin, California, Investigated with the SBAS-InSAR Approach[J]. Remote Sensing of Environment, 2007, 108(3): 277-289.

[3]尹宏杰, 朱建軍, 李志偉, 等. 基于SBAS的礦區(qū)形變監(jiān)測研究[J]. 測繪學(xué)報(bào), 2011, 40(1): 52-58.

YIN Hongjie, ZHU Jianjun, LI Zhiwei, et al. Ground Subsidence Monitoring in Mining Area Using DInSAR SBAS Algorithm[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(1): 52-58.

[4]李珊珊, 李志偉, 胡俊, 等. SBAS-InSAR技術(shù)監(jiān)測青藏高原季節(jié)性凍土形變[J]. 地球物理學(xué)報(bào), 2013, 56(5): 1476-1486.

LI Shanshan, LI Zhiwei, HU Jun, et al. Investigation of the Seasonal Oscillation of the Permafrost over Qinghai-Tibet Plateau with SBAS-InSAR Algorithm[J]. Chinese Journal of Geophysics, 2013, 56(5): 1476-1486.

[5]吳宏安, 張永紅, 陳曉勇, 等. 基于小基線DInSAR技術(shù)監(jiān)測太原市2003—2009年地表形變場[J]. 地球物理學(xué)報(bào), 2011, 54(3): 673-680.

WU Hongan, ZHANG Yonghong, CHEN Xiaoyong, et al. Ground Deformation Monitoring Using Small Baseline DInSAR Technique: A Case Study in Taiyuan City from 2003 to 2009[J]. Chinese Journal of Geophysics, 2011, 54(3): 673-680.

[6]許文斌, 李志偉, 丁曉利, 等. 利用InSAR短基線技術(shù)估計(jì)洛杉磯地區(qū)的地表時序形變和含水層參數(shù)[J]. 地球物理學(xué)報(bào), 2012, 55(2): 452-461.

XU Wenbin, LI Zhiwei, DING Xiaoli, et al. Application of Small Baseline Subsets D-InSAR Technology to Estimate the Time Series Land Deformation and Aquifer Storage Coefficients of Los Angeles Area[J]. Chinese Journal of Geophysics, 2012, 55(2): 452-461.

[7]張永紅,張繼賢,龔文瑜,等.基于SAR干涉點(diǎn)目標(biāo)分析技術(shù)的城市地表形變監(jiān)測[J]. 測繪學(xué)報(bào), 2009, 38(6): 482-487, 493. DOI: 10.3321/j.issn:1001-1595.2009.06.003.

ZHANG Yonghong, ZHANG Jixian, GONG Wenyu, et al. Monitoring Urban Subsidence Based on SAR Interferometric Point Target Analysis[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(6): 482-487, 493. DOI: 10.3321/j.issn:1001-1595.2009.06.003.

[8]張永紅, 吳宏安, 孫廣通. 時間序列InSAR技術(shù)中的形變模型研究[J]. 測繪學(xué)報(bào), 2012, 41(6): 864-869, 876. ZHANG Yonghong, WU Hongan, SUN Guangtong. Deformation Model of Time Series Interferometric SAR Techniques[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(6): 864-869, 876.

[9]FERRETTI A, PRATI C, ROCCA F. Nonlinear Subsidence Rate Estimation Using Permanent Scatterers in Differential SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(5): 2202-2212.

[10]李永生, 張景發(fā), 李振洪, 等. 利用短基線集干涉測量時序分析方法監(jiān)測北京市地面沉降[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2013, 38(11): 1374-1377.

LI Yongsheng, ZHANG Jingfa, LI Zhenhong, et al. Land Subsidence in Beijing City from InSAR Time Series Analysis with Small Baseline Subset[J]. Geomatics and Information Science of Wuhan University, 2013, 38(11): 1374-1377.

[11]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.

[12]BONANO M, MANUNTA M, MARSELLA M, et al. Long-term ERS/ENVISAT Deformation Time-series Generation at Full Spatial Resolution via the Extended SBAS Technique[J]. International Journal of Remote Sensing, 2012, 33(15): 4756-4783.

[13]AKBARI V, MOTAGH M. Improved Ground Subsidence Monitoring Using Small Baseline SAR Interferograms and a Weighted Least Squares Inversion Algorithm[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(3): 437-441.

[14]LAUKNES T R, ZEBKER H A, LARSEN Y. InSAR Deformation Time Series Using an L1-norm Small-baseline Approach[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(1): 536-546.

[15]王家映. 地球物理反演理論[M]. 北京: 高等教育出版社, 2002.

WANG Jiaying. Inverse Theory in Geophysics[M]. Beijing: Higher Education Press, 2002.

[16]TIKHONOV A N, ARSENIN V Y. Solutions of Ill-posed Problems[M]. New York: John Wiley and Sons, 1977.

[17]HE Guoqiang, LIU Linxian. A Kind of Implicit Iterative Methods for Ill-posed Operator Equation[J]. Journal of Computational Mathematics, 1999, 17(3): 275-284.

[18]STEIHAUG T.The Conjugate Gradient Method and Trust Regions in Large Scale Optimization[J]. SIAM Journal on Numerical Analysis, 1983, 20(3): 626-637.

[19]郭建峰. 測量平差系統(tǒng)病態(tài)性的診斷與處理[D]. 鄭州: 信息工程大學(xué), 2002.

GUO Jianfeng. Diagnostics and Processing of Ill-posed in Surveying Adjustment System[D]. Zhengzhou: Information Engineering University, 2002.

[20]王松桂, 賈忠貞. 矩陣論中不等式[M]. 合肥: 安徽教育出版社, 1994.

WANGSonggui,JIAZhongzhen.InequalityinMatrixTheory[M].Hefei:AnhuiEducationPress, 1994.

[21]王彥飛. 反演問題的計(jì)算方法及其應(yīng)用[M]. 北京: 高等教育出版社, 2007.

WANGYanfei.ComputationalMethodsforInverseProblemsandTheirApplications[M].Beijing:HigherEducationPress, 2007.

[22]王振杰. 測量中不適定問題的正則化解法[M].北京:科學(xué)出版社, 2006.

WANGZhenjie.RegularizationSolutionsofIll-posedProblemsinGeodesy[M].Beijing:SciencePress, 2006.

[23]楊文采. 地球物理反演和地震層析成像[M]. 北京: 地質(zhì)出版社, 1989.

YANGWencai.GeophysicalInversionandSeismicTomography[M].Beijing:GeologicalPress, 1989.

[24]張方仁, 汪曉慶. 平差參數(shù)的嶺估計(jì)和壓縮估計(jì)[J]. 武漢測繪科技大學(xué)學(xué)報(bào), 1989, 14(3): 46-58.

ZHANGFangren,WANGXiaoqing.TheRidgeEstimationandtheSteinEstimationofAdjustmentParameters[J].JournalofWuhanTechnicalUniversityofSurveyingandMapping, 1989, 14(3): 46-58.

[25]陶秋香, 劉國林.PSInSAR公共主影像優(yōu)化選取的一種新方法[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2011, 36(12): 1456-1460.

TAOQiuxiang,LIUGuolin.ANewMethodtoOptimizeandSelectCommonMasterImagesinPSInSAR[J].GeomaticsandInformationScienceofWuhanUniversity, 2011, 36(12): 1456-1460.

[26]PEPEA,LANARIR.OntheExtensionoftheMinimumCostFlowAlgorithmforPhaseUnwrappingofMultitemporalDifferentialSARInterferograms[J].IEEETransactionsonGeoscienceandRemoteSensing, 2006, 44(9): 2374-2383.

(責(zé)任編輯：叢樹平)

Regularization solution of Small Baseline Subset Deformation Model Inversion

JIANG Zhaoying1,2，LIU Guolin1，TAO Qiuxiang1

1. Geomatics College, Shandong University of Science and Technology, Qingdao 266590, China; 2. College of Science and Information, Qingdao Agricultural University, Qingdao 266109, China

For the coefficient matrix of the normal equation is ill-conditioned during inverting deformation model of small baseline subset (SBAS) InSAR technique, a regularization robust method is proposed. Based on Tikhonov regularization theory, this method converts the problem of how to solve the deformation rate into minimization problem. According to L-curve method to choose regularization parameter, considering the relationship between the individual components of least-squares residuals to choose regularization matrix, thus it achieves robust solution of SBAS deformation model inversion. We adopt respectively least-squares estimation, ridge estimation and Tikhonov regularization method to deal with 29 ENVISAT ASAR dataset relevant to the Beijing area, achieving the subsidence rate map of the study area. Through comparative analysis among the mean square error (MSE) of 21 points on behalf of the different subsidence, temporal coherence values and MSE maps of the entire study area, we confirm that Tikhonov regularization robust method in inverting SBAS deformation model can obtain more reliable results of deformation monitoring.

small baseline subset; condition number; Tikhonov regularization; ridge estimation; mean squared error

10.11947/j.AGCS.2016.20150143.

JIANG Zhaoying，LIU Guolin，TAO Qiuxiang.Regularization solution of Small Baseline Subset Deformation Model Inversion [J]. Acta Geodaetica et Cartographica Sinica,2016,45(5):566-573. DOI:10.11947/j.AGCS.2016.20150143.

P226

A

1001-1595(2016)05-0566-08

國家自然科學(xué)基金(41274007;41404003);山東省自然科學(xué)基金(ZR2012DM001);高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金(20123718110001);山東省泰山學(xué)者建設(shè)工程專項(xiàng)經(jīng)費(fèi)(TSXZ201509)

引文格式：姜兆英，劉國林，陶秋香.短基線集形變模型反演的正則化解算方法[J].測繪學(xué)報(bào)，2016,45(5)：566-573.