韋建超 李志偉 段夢 曹云夢 馮光財 胡俊

摘要軌道誤差和長波大氣延遲組成的系統(tǒng)誤差是影響InSAR形變監(jiān)測精度的重要因素之一.傳統(tǒng)方法在空間域?qū)Ω缮鎴D的系統(tǒng)誤差建模，容易導致長波形變和系統(tǒng)誤差相混淆.本文在時空域利用附加系統(tǒng)參數(shù)對系統(tǒng)誤差建模，同時根據(jù)觀測值質(zhì)量對差分相位觀測值定權，采用附加系統(tǒng)參數(shù)的加權最小二乘法估計形變參數(shù)和系統(tǒng)誤差，實現(xiàn)了長波形變和系統(tǒng)誤差的分離.模擬實驗結果表明，在形變與系統(tǒng)誤差的空間變化特性完全一致的極端情況下，本文方法能實現(xiàn)兩者的有效分離，估計的形變速率均方根誤差比傳統(tǒng)方法降低了98.8%.ASAR數(shù)據(jù)實驗顯示當形變尺度較小且分散分布時，本文方法和傳統(tǒng)方法得到的結果相似;當形變在研究區(qū)內(nèi)表現(xiàn)為長波變化時，本文方法比傳統(tǒng)方法估計的形變結果更為穩(wěn)健.關鍵詞附加系統(tǒng)參數(shù);多時相InSAR;時空建模;形變估計

中圖分類號P237

文獻標志碼A

0引言

合成孔徑雷達干涉（InSAR）為研究大范圍地表形變提供了強有力的工具，近年來被廣泛應用于板塊運動[1-2]、地震[3-6] 、地下水抽取[7]等各種自然或人為因素導致的地表沉降研究中.當InSAR用于形變監(jiān)測時，不同誤差源會使測量結果產(chǎn)生偏差，本文主要討論軌道誤差和長波大氣延遲情況下的形變估計.如果構成差分干涉圖的兩個觀測時刻軌道參數(shù)含有誤差，在干涉圖中軌道誤差會以系統(tǒng)性的殘差條紋存在.完全去除干涉圖軌道誤差條紋要求軌道絕對精度優(yōu)于1 mm [8]，而目前ERS和Envisat衛(wèi)星精密軌道數(shù)據(jù)在法向、沿軌向和交軌向的軌道誤差分別為7、24和18 cm [9]，ALOS和Radarsat衛(wèi)星沒有提供精密軌道數(shù)據(jù)，因此基于衛(wèi)星定軌數(shù)據(jù)的軌道誤差改正遠不能滿足測量精度要求.此外，構成干涉圖的兩個觀測時刻大氣延遲不一致，導致干涉圖中存在附加的大氣延遲相位，其主要包含三部分[10]：1）由氣壓、溫度和相對濕度隨高度變化引起的地形相關的大氣延遲;2）由氣壓、溫度和濕度橫向變化而導致的長波（10千米以上）大氣延遲;3）由大氣邊界層湍流導致的小尺度（幾千米以下）大氣延遲.其中高程相關的大氣延遲可以通過高程相關的模型進行建模[11-13]，湍流延遲可通過建立隨機模型進行處理[14-15].長波大氣延遲在空間上的變化與軌道誤差相似，均呈系統(tǒng)性的梯度變化，為InSAR觀測的系統(tǒng)性誤差.

目前針對InSAR觀測的系統(tǒng)誤差建模方法主要分為三類：1）頻率域建模法，包括通過傅里葉變換[16]和小波分析[17]來改正干涉圖的系統(tǒng)誤差.該方法難以區(qū)分長波形變和系統(tǒng)誤差.2）利用衛(wèi)星軌道狀態(tài)矢量信息來獲取基線參數(shù)[18-19].該方法在高相干區(qū)域能估計軌道參數(shù)并在一定程度上改正軌道誤差，但長波大氣延遲依然存在.3）空間域建模法.該方法在空間域上建立干涉圖相位和空間位置的一階或高階函數(shù)關系[8，20-21].該方法同樣存在長波形變和系統(tǒng)誤差不能區(qū)分的問題.而多時相InSAR除了能克服傳統(tǒng)InSAR技術的失相干問題外，也為系統(tǒng)誤差的建模和估計提供了可能[22].本文用附加系統(tǒng)參數(shù)對每景SAR影像的系統(tǒng)誤差進行建模，對受湍流延遲和失相干噪聲影響的觀測值定權，用附加系統(tǒng)參數(shù)的加權最小二乘方法同時估計系統(tǒng)誤差和形變參數(shù).與現(xiàn)有方法相比，本文方法對研究區(qū)的形變類型不敏感，在長波形變時也能正確估計系統(tǒng)誤差和形變，增加了形變參數(shù)估計的精度和可靠性.

1多時域InSAR的時空函數(shù)模型

根據(jù)合成孔徑雷達干涉原理，點P在解纏繞干涉圖中的差分相位可表示為

φp=φdef，p+φtopo_res，p+φorb，p+φatm，p+φdec_noi，p， ?（1）

公式（1）中φdef，p為形變相位.假設形變的時空模型為

ldef=f（x，y，t|aτ）， （2）

其中x，y表示P的空間位置，t為時間，aτ為形變模型的τ個待求參數(shù)向量.因變量x，y可選，當形變在空間上的變化未知時，f（x，y，t|aτ）可簡化為f（t|aτ）.在tm時刻的主影像和ts時刻的從影像構成的干涉圖中，P的形變相位可表示為

φdef，p=4πλ（f（xp，yp，tm）-f（xp，yp，ts））aτ. ?（3）

式（1）中地形殘余相位φtopo_res，p可以寫成：

φtopo_res，p=4π·Bperp，pλ·R·sinθ·ΔhP， （4）

其中ΔhP為P點的地形改正參數(shù)，Bperp，p為干涉圖的垂直基線，θ為入射角，R為主影像雷達天線與地面目標之間的距離.式（1）中的φorb，p為軌道誤差相位、φatm，p為大氣延遲相位.φatm，p包含了地形相關的垂直分層延遲項φatm_ver，p

、長波延遲項φatm_long，p

和小尺度的湍流延遲項φatm_tur，p

[10]：

φatm，p=φatm_ver，p+φatm_long，p+φatm_tur，p， ?（5）

其中垂直分層延遲用高程回歸模型進行建模和改正[11];湍流延遲和失相干噪聲φdec_noi，p用隨機模型建模[14，15];而長波大氣延遲項φatm_long，p

和軌道誤差φorb，p在干涉圖上具有系統(tǒng)性的長波變化，構成多時相InSAR的系統(tǒng)誤差項.假設各成像時刻系統(tǒng)誤差項互相獨立，基于軌道誤差的非線性特性[23]，各個獨立觀測時刻的系統(tǒng)誤差用以下二維多項式建模：

trend=b1x+b2y+b3xPyP+…. （6）

為避免高階多項式帶來的龍格現(xiàn)象[21]，式（6）最多取至二次項，令m=[xyxyx2y2]，r=[b1b2b3b4b5]T，式（6）寫成矩陣形式有：

trend=mr， ?（7）

則p在干涉圖中的系統(tǒng)誤差相位為

φtrend，ms=trend，m-trend，s=m（rm-rs）. （8）

將式（3）、（4）和（8）代入式（1），得到包含形變模型參數(shù)、地形殘差參數(shù)和附加系統(tǒng)參數(shù)的觀測相位方程：

φ=4π·Bperpλ·R·sinθ·Δh+4πλ（f（x，y，tm）-

f（x，y，ts））aτ+m（rm-rs）+Δp， ?（9）

其中Δp主要包含湍流延遲和失相干噪聲的隨機誤差相位.

如果多時相InSAR包含由N景SAR影像組成的M幅干涉圖，同時對干涉圖中的S個目標點進行分析（如圖1所示），則共有M×S個觀測方程.寫成矩陣形式有：

ΨMS×1=AMS×（S+τ）X（S+τ）×1+BMS×uN YuN×1+WMS×1， ?（10）

其中Ψ為差分干涉相位觀測值向量;W是包括大氣湍流延遲和失相干噪聲在內(nèi)的觀測噪聲;X為包含地形改正向量Δh和形變模型參數(shù)向量aτ的待求參數(shù)向量;Y為由各個SAR影像的附加系統(tǒng)參數(shù)向量r構成的附加系統(tǒng)參數(shù)向量，下標u為式（7）所示的系統(tǒng)誤差向量r中參數(shù)的個數(shù);A和B分別為待求參數(shù)向量和附加系統(tǒng)參數(shù)向量的系數(shù)矩陣.如果隨機模型為：

DΨ=σ20P-1Ψ， （11）

σ0為單位權中誤差，PΨ為觀測值的權陣.如果僅考慮相干性對觀測質(zhì)量的影響，可以根據(jù)相干性估計各差分相位觀測值的方差[24]：

σ2i=（1-μ2i）/（2Lμ2i）， ?（12）

其中μi為i處的相干性，L為相干性估計的多視數(shù)，再根據(jù)方差定權.如果考慮湍流延遲以及失相干噪聲對觀測質(zhì)量的影響，可以建立顧及湍流延遲和失相干噪聲的方差協(xié)方差[14-15].最后按附加系統(tǒng)參數(shù)的加權最小二乘準則平差，則有：

ATPψA

ATPψB

BTPψA

BTPψB

=

ATPψΨ

BTPψΨ. （13）

令N11=ATPΨA，N12=ATPΨB，N21=BTPΨA，N22=BTPΨB，H=N22-N21N-111N12，可以求解形變參數(shù)和附加系統(tǒng)參數(shù)：

=

N-111+N-111N12H-1N21N-111

-N-111N12M-1

-H-1N21N-111H-1

ATPΨΨ

BTPΨΨ. ?（14）

2實例分析

2.1模擬實驗

形變在研究區(qū)內(nèi)為長波形變的情形常見于地質(zhì)構造引起的地表形變，或者研究區(qū)位于地震斷裂帶某一側(cè)時等情形.本文模擬了極端情況下，即形變和系統(tǒng)誤差在空間上的變化完全一致，都可以用式（6）表示的長波變化時的形變.在時序上，形變隨時間線性變化時的情形.形變的時空模型為

f（x，y，t|aτ）=（a1x+a2y+a3xPyP）t. ?（15）

對應的形變相位為

φdef=x（tm-ts）y（tm-ts）xy（tm-ts）·aτ， ?（16）

而時序上互相獨立的系統(tǒng)誤差相位可表示為

φtrend=[…x…-x…y…-y…xy…-xy]·r. ?（17）

可見即使形變參數(shù)向量aτ和附加系統(tǒng)參數(shù)向量r在空間上的變化完全一致，但兩者在時空函數(shù)模型中的系數(shù)矩陣并不相同.當目標點S足夠多時，即可同時估計出形變參數(shù)和系統(tǒng)誤差系數(shù).此外，本文模擬了湍流大氣延遲、地形殘差、軌道誤差及失相干噪聲，最后得到包含各種誤差的兩組完全一樣的纏繞干涉圖序列（圖1）.接著按圖2所示的實驗流程，分別用空間域建模法（傳統(tǒng)方法）和附加系統(tǒng)參數(shù)的時空函數(shù)模型（本文方法）估計系統(tǒng)誤差系數(shù)和形變速率參數(shù)，并分別計算系統(tǒng)誤差系數(shù)和形變速率參數(shù)的均方根誤差（RMSE），作為參數(shù)估計精度的評價，在模擬實驗中，對所有觀測值取相同的權.由于傳統(tǒng)方法在空間域建模時只能求出各干涉圖的差分系統(tǒng)誤差系數(shù)，為便于比較，先將本文方法得到的各個SAR影像的系統(tǒng)誤差系數(shù)按干涉相對進行差分，得到與干涉圖對應的差分系統(tǒng)誤差系數(shù)后再與傳統(tǒng)方法的結果比較.圖3顯示了干涉圖在距離向x上的差分系統(tǒng)誤差參數(shù)估值Δ1，干涉圖在方位向y上的差分系統(tǒng)誤差參數(shù)估值Δ2，以及二次項xy上的差分系統(tǒng)誤差參數(shù)估值Δ3.本文方法得到的Δ1，Δ2和Δ3的RMSE值比傳統(tǒng)方法分別降低了97.54%、92.45%和12.90%.圖4顯示了沉降參數(shù)的估計結果，傳統(tǒng)方法得到的形變速率在空間變化趨勢上與模擬值明顯不符，RMSE為60.92 mm/a;而本文方法得到的形變速率與模擬值更接近，RMSE為0.67 mm/a.

2.2真實數(shù)據(jù)實驗

本文選擇中心位置位于39°50′N，116°38′E，面積約為4 000 km2的實驗區(qū)（圖5），該區(qū)域平均海拔為43.5 m，地勢較為平坦.對覆蓋實驗區(qū)的33景時間跨度從2007年1月3日至2010年9月29日的ASAR影像，利用GAMMA軟件進行如下處理：

1）以2009年1月7日獲取的影像為公共影像，將其他所有影像進行配準并重采樣至公共影像.

2）以382 m為垂直基線閾值，280 d為時間基線閾值對多時相InSAR進行構網(wǎng)，優(yōu)化網(wǎng)形最終生成61對小基線差分干涉對（圖6）.

3）用30 s空間分辨率的SRTM DEM數(shù)據(jù)消除地形相位;在方位向和距離向分別進行20和4的多視;采用改進的Goldstein濾波方法[25]對干涉圖進行濾波，并采用最小費用流法對相干性高于0.7的區(qū)域進行解纏，最后得到多時相InSAR的解纏相位序列.

4）對解纏繞的差分相位觀測值分別用傳統(tǒng)方法和本文方法估計其形變速率參數(shù)，同時考慮了湍流延遲和失相干噪聲對觀測質(zhì)量的影響，用文獻[14]介紹的方法估計觀測值的時空方差協(xié)方差陣用于參數(shù)估計的定權.由于相干性大于0.7的區(qū)域絕大部分位于平原地帶（圖7），且最大高差僅為438 m，因此不考慮隨高度變化的大氣垂直分層延遲影響.

最終得到的形變速率結果如圖7所示，正值（紅色）表示從影像相對于主影像在雷達視線方向的上升，負值（綠色）表示從影像相對于主影像下降.兩種方法估計的形變速率都處在0～145 mm/a的區(qū)間，其中沉降區(qū)都主要分布于潮白河、溫榆河和泃河流域的沖積和洪積扇平原上，空間分布不均勻.位于昌平、順義、朝陽和通州的沉降區(qū)有逐漸連成一片同時向東擴張的趨勢，初步解釋為與城市擴張相關.中心城區(qū)、海淀區(qū)、房山區(qū)、豐臺區(qū)與石景山區(qū)的地表沉降速率都在10 mm/a，而且較為穩(wěn)定;昌平區(qū)的沙河鎮(zhèn)與上莊鎮(zhèn)出現(xiàn)了沉降漏斗，年沉降速率超過70 mm/a;朝陽區(qū)金盞鄉(xiāng)一帶沉降明顯，最大平均沉降速率接近80 mm/a.圖7c為兩種方法估計的形變速率之差，差值在10 mm/a以內(nèi)，兩種方法得到形變速率在量級和空間格局上都非常接近，

本文接著選擇位于沉降區(qū)一側(cè)的一個局部區(qū)域作為局部研究區(qū)（圖5和圖7c中紅色矩形區(qū)域），假設只有該區(qū)域能獲得高相干觀測值，其在圖7a和7b中的形變速率向量分別用向量V1，V2表示.由于V1和V2的差值在6 mm/a以內(nèi)，本文取=（V1+V2）/2作為該研究區(qū)的形變速率參考值（圖8a），分別用傳統(tǒng)方法和本文方法估計該區(qū)域的形變速率向量V′1 和V′2 （圖8b、圖8c）.可見當以該局部區(qū)域作為研究區(qū)時，傳統(tǒng)方法估計的沉降速率V′1 與參考值存在一個長波變化的偏差.而本文方法得到的V′2 與參考值更為接近.如果用Vi和V′i 的相關性評估形變估計的穩(wěn)健性，用V′i 的RMSE值來評估兩種方法的精度，結果見圖8.本文方法比空間域建模法在穩(wěn)健性方面提高了100%，在RMSE值上降低了59.73%.

3討論

模擬實驗結果表明，當研究區(qū)域內(nèi)的形變在空間上表現(xiàn)為長波緩慢形變時，本文方法得到的結果更加穩(wěn)健，即使在形變和系統(tǒng)誤差在空間上的變化完全一致的極端情況下，本文方法也能實現(xiàn)形變和系統(tǒng)誤差的有效分離.而在真實數(shù)據(jù)實驗中，當以整個干涉圖覆蓋區(qū)域作為研究對象時，由于形變相對于研究區(qū)域的尺度更小，且分散的分布在研究區(qū)的不同位置，

在整個研究區(qū)域內(nèi)不是長波變化，使用傳統(tǒng)方法和本文方法得到的形變速率較為接近.而當以標識的局部區(qū)域作為研究區(qū)時，形變相對于研究區(qū)域為大尺度的長波變化，傳統(tǒng)方法在估計系統(tǒng)誤差參數(shù)時明顯受到形變的影響，一部分長波形變被誤當成系統(tǒng)誤差進行了補償.而本文方法得到的形變速率在空間分布和數(shù)值上與參考值更為接近，說明本文方法受形變類型的影響更小，其形變參數(shù)的估計結果更穩(wěn)健.造成這種結果的主要原因如下：

1）空間域建模法在估計系統(tǒng)參數(shù)時，只是在空間維度對每個干涉圖利用二維多項式擬合，將形變、地形誤差項、大氣延遲和失相干等分量視為隨機誤差.而實際上這些量在空間上不一定滿足零均值隨機分布的特點，尤其是干涉圖中的長波形變會對系統(tǒng)誤差系數(shù)的估計造成干擾，易將長波形變誤當成系統(tǒng)誤差，造成系統(tǒng)誤差的過估計.而附加系統(tǒng)參數(shù)的時空函數(shù)模型在時空維度上對系統(tǒng)誤差和長波形變分別建模，即使兩者在空間上的變化完全一致，也能根據(jù)兩者在時間上的不同特性進行區(qū)分，從而實現(xiàn)將長波形變從系統(tǒng)誤差中有效分離，避免形變的錯誤估計.

2）空間域建模法對每個干涉圖單獨估計系統(tǒng)參數(shù)時，假設各干涉圖互相獨立.而實際上干涉圖之間包含有公共影像，不滿足獨立條件.附加系統(tǒng)參數(shù)法對每一個觀測時刻的系統(tǒng)誤差賦予一組獨立的系統(tǒng)參數(shù)，差分組合的誤差系數(shù)由觀測時刻的系數(shù)計算得到，在數(shù)據(jù)的解算過程中增加了閉合差約束，因此得到的各干涉圖的系統(tǒng)參數(shù)符合性更好.

3）空間域建模法沒有考慮不同相干性區(qū)域差分相位的觀測質(zhì)量對結果的影響，而本文方法可以根據(jù)相干性、大氣湍流延遲和失相干噪聲對觀測質(zhì)量的影響來給相位觀測值合理定權.因此，其得到的結果和實際更加相符.

總而言之，隨著覆蓋范圍更小的高分辨率SAR影像（如TerraSAR）數(shù)據(jù)的應用越來越廣，或者在多時相InSAR應用中可能存在可用于分析建模的高相干區(qū)域只覆蓋局部的形變區(qū)域（如本文實例中選的紅色矩形區(qū)域），此時形變在研究區(qū)域內(nèi)可能呈現(xiàn)出長波的變化趨勢.這種情況下傳統(tǒng)方法會導致系統(tǒng)誤差的過估計和形變的低估計，因此，有必要采用能區(qū)分長波形變和系統(tǒng)誤差，對形變類型不敏感的本文方法來獲得穩(wěn)健的形變估計結果.

4結論

本文系統(tǒng)介紹了附加系統(tǒng)參數(shù)的多時相InSAR時空建模和形變估計方法.該方法給每個觀測時刻的軌道誤差和長波大氣延遲賦予一組獨立的附加系統(tǒng)參數(shù)建立模型，對形變建立包含形變參數(shù)的時空函數(shù)模型，在時空維建立多時相InSAR差分相位的觀測模型，利用附加系統(tǒng)參數(shù)的加權最小二乘平差方法同時估計系統(tǒng)誤差系數(shù)和形變參數(shù).通過模擬和真實數(shù)據(jù)實驗驗證了附加系統(tǒng)參數(shù)的多時相InSAR時空函數(shù)模型在有效估計系統(tǒng)參數(shù)的同時，很好地保留了長波趨勢的形變信號.與傳統(tǒng)多時相InSAR去軌道誤差方法相比，它增加了系統(tǒng)參數(shù)閉合差的約束，同時也減少了長波形變對系統(tǒng)參數(shù)估計的影響，還考慮了不同觀測質(zhì)量差分相位的權重.雖然本文方法還存在著只針對解纏繞相位、增加計算復雜性等不足，但隨著覆蓋率更小的高分辨率SAR影像的廣泛應用，以及在多時相InSAR應用中可能會出現(xiàn)的較小范圍的高相干區(qū)域，導致研究區(qū)域內(nèi)形變在空間上呈長波變化的可能性增加，本文的研究增加了利用InSAR估計長波形變的穩(wěn)健性.

參考文獻

References

[1]Champenois J，F(xiàn)runeau B，Pathier E，et al.Monitoring of active tectonic deformations in the Longitudinal Valley （Eastern Taiwan） using Persistent Scatterer InSAR method with ALOS PALSAR data[J].Earth and Planetary Science Letters，2012，337/338：144-155

[2]Bürgmann R，Hilley G，F(xiàn)erretti A，et al.Resolving vertical tectonics in the San Francisco Bay Area from permanent scatterer InSAR and GPS analysis[J].Geology，2006，34（3）：221-224

[3]Massonnet D，Rossi M，Carmona C，et al.The displacement field of the Landers earthquake mapped by radar interferometry[J].Nature，1993，364（6433）：138-142

[4]Massonnet D，F(xiàn)eigl K，Rossi M，et al.Radar interferometric mapping of deformation in the year after the Landers earthquake[J].Nature，1994，369（6477）：227-230

[5]Zebker H A，Rosen P A，Goldstein R M，et al.On the derivation of coseismic displacement fields using differential radar interferometry：the Landers earthquake[J].Journal of Geophysical Research：Solid Earth，1994，99（B10）：19617-19634

[6]Feng G C，Hetland E A，Ding X L，et al.Coseismic fault slip of the 2008 Mw7.9 Wenchuan earthquake estimated from InSAR and GPS measurements[J].Geophysical Research Letters，2010，37（1）：73-78

[7]Lu Z，Danskin W R.InSAR analysis of natural recharge to define structure of a ground-water basin，San Bernardino，California[J].Geophysical Research Letters，2001，28（13）：2661-2664

[8]Hanssen R F.Radar interferometry：data interpretation and error analysis[M].Springer Science & Business Media，2001

[9]Scharroo R，Visser P.Precise orbit determination and gravity field improvement for the ERS satellites[J].Journal of Geophysical Research：Oceans，1998，103（C4）：8113-8127

[10]Bekaert D P S，Walters R J，Wright T J，et al.Statistical comparison of InSAR tropospheric correction techniques[J].Remote Sensing of Environment，2015，170：40-47

[11]Onn F，Zebker H A.Correction for interferometric synthetic aperture radar atmospheric phase artifacts using time series of zenith wet delay observations from a GPS network[J].Journal of Geophysical Research Atmospheres，2006，111（B9）：B09102

[12]Li Z W，Xu W B，F(xiàn)eng G C，et al.Correcting atmospheric effects on InSAR with MERIS water vapour data and elevation-dependent interpolation model[J].Geophysical Journal International，2012，189（2）：898-910

[13]Bekaert D P S，Hooper A，Wright T J.A spatially variable power law tropospheric correction technique for InSAR data[J].Journal of Geophysical Research：Solid Earth，2015，120（2）：1345-1356

[14]Cao Y M，Li Z W，Wei J C，et al.Stochastic modeling for time series InSAR：with emphasis on atmospheric effects[J].Journal of Geodesy，2018，92（2）：185-204

[15]Agram P S，Simons M.A noise model for InSAR time series[J].Journal of Geophysical Research：Solid Earth，2015，120（4）：2752-2771

[16]朱珺，丁曉利，許兵，等.基于條紋頻率檢測的干涉圖去線性趨勢精化方法[J].大地測量與地球動力學，2011，31（3）：138-141

ZHU Jun，DING Xiaoli，XU Bing，et al.Refine method for removing the linear trend in interferogram based on fringe frequency test[J].Journal of Geodesy and Geodynamics，2011，31（3）：138-141

[17]Shirzaei M，Walter T R.Estimating the effect of satellite orbital error using wavelet-based robust regression applied to InSAR deformation data[J].IEEE Transactions on Geoscience and Remote Sensing，2011，49（11）：4600-4605

[18]Kohlhase A O，F(xiàn)eigl K L，Massonnet D.Applying differential InSAR to orbital dynamics：a new approach for estimating ERS trajectories[J].Journal of Geodesy，2003，77（9）：493-502

[19]Bhr H，Hanssen R F.Reliable estimation of orbit errors in spaceborne SAR interferometry[J].Journal of Geodesy，2012，86（12）：1147-1164

[20]Ferretti A，Montiguarnieri A，Prati C，et al.InSAR principles：guidelines for SAR interferometry processing and interpretation[J].Journal of Financial Stability，2007，10（10）：156-162

[21]Xu B，Li Z W，Wang Q J，et al.A refined strategy for removing composite errors of SAR interferogram[J].IEEE Geoscience and Remote Sensing Letters，2014，11（1）：143-147

[22]Zhang L，Ding X L，Lu Z，et al.A novel multitemporal InSAR model for joint estimation of deformation rates and orbital errors[J].IEEE Transactions on Geoscience and Remote Sensing，2014，52（6）：3529-3540

[23]Liu G，Hanssen R F，Guo H D，et al.Nonlinear model for InSAR baseline error[J].IEEE Transactions on Geoscience and Remote Sensing，2016，54（9）：5341-5351

[24]Rodriguez E，Martin J M.Theory and design of interferometric synthetic aperture radars[J].IEE Proceedings F Radar and Signal Processing，1992，139（2）：147

[25]Li Z W，Ding X L，Huang C，et al.Improved filtering parameter determination for the Goldstein radar interferogram filter[J].ISPRS Journal of Photogrammetry and Remote Sensing，2008，63（6）：621-634

Spatio-temporal modeling and deformation estimation of multi-temporal

InSAR with additional systematical parameters

WEI Jianchao1，2LI Zhiwei2DUAN Meng2CAO Yunmeng2FENG Guangcai2HU Jun2

1School of Resource & Environment and Safety Engineering，Hunan University of Science and Technology，Xiangtan411201

2School of Geosciences and Info-Physics，Central South University，Changsha410083

AbstractThe systematic error of InSAR is composed of orbit error and long-wavelength atmospheric delay，which is one of the most important factors that affect the accuracy of InSAR deformation monitoring.Traditional method models the systematic error of each individual interferogram in spatial domain，so it is difficult to distinguish between the long-wavelength deformation and systematic error.In this study，we modeled the systematic errors in spatio-temporal domain，determined the weight matrix of the observations according to its quality，and estimated the deformation parameters and systematic error by the weighted least square method.The simulation results show that，even in the extreme cases where the spatial characteristics of deformation and systematic error are completely consistent，the proposed method can separate the deformation and systematic errors effectively.The deformation rate RMSE of the proposed method is 98.8% lower than that of traditional method.The experiment with ASAR shows that，the proposed method and traditional method obtain similar results when the deformation scale is small and scattered over the study area.However，when large scale long-wavelength deformation occurs in the study area，the results obtained by the proposed method are more robust than those obtained by traditional method.

Key wordsadditional systematical parameters;multi-temporal InSAR;spatio-temporal modeling;deformation estimation

收稿日期2019-10-15

資助項目國家自然科學基金（41474007，41404013，41222027）

作者簡介韋建超，男，博士，研究方向為合成孔徑雷達干涉測量的誤差建模與參數(shù)反演.wjchao1608@163.com

1湖南科技大學資源環(huán)境與安全工程學院，湘潭，411201

2中南大學地球科學與信息物理學院，長沙，410083