韋建超 李志偉 段夢 曹云夢 馮光財(cái) 胡俊

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

中圖分類號(hào)P237

文獻(xiàn)標(biāo)志碼A

0引言

合成孔徑雷達(dá)干涉（InSAR）為研究大范圍地表形變提供了強(qiáng)有力的工具，近年來被廣泛應(yīng)用于板塊運(yùn)動(dòng)[1-2]、地震[3-6] 、地下水抽取[7]等各種自然或人為因素導(dǎo)致的地表沉降研究中.當(dāng)InSAR用于形變監(jiān)測時(shí)，不同誤差源會(huì)使測量結(jié)果產(chǎn)生偏差，本文主要討論軌道誤差和長波大氣延遲情況下的形變估計(jì).如果構(gòu)成差分干涉圖的兩個(gè)觀測時(shí)刻軌道參數(shù)含有誤差，在干涉圖中軌道誤差會(huì)以系統(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ù)的軌道誤差改正遠(yuǎn)不能滿足測量精度要求.此外，構(gòu)成干涉圖的兩個(gè)觀測時(shí)刻大氣延遲不一致，導(dǎo)致干涉圖中存在附加的大氣延遲相位，其主要包含三部分[10]：1）由氣壓、溫度和相對濕度隨高度變化引起的地形相關(guān)的大氣延遲;2）由氣壓、溫度和濕度橫向變化而導(dǎo)致的長波（10千米以上）大氣延遲;3）由大氣邊界層湍流導(dǎo)致的小尺度（幾千米以下）大氣延遲.其中高程相關(guān)的大氣延遲可以通過高程相關(guān)的模型進(jìn)行建模[11-13]，湍流延遲可通過建立隨機(jī)模型進(jìn)行處理[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ū)域能估計(jì)軌道參數(shù)并在一定程度上改正軌道誤差，但長波大氣延遲依然存在.3）空間域建模法.該方法在空間域上建立干涉圖相位和空間位置的一階或高階函數(shù)關(guān)系[8，20-21].該方法同樣存在長波形變和系統(tǒng)誤差不能區(qū)分的問題.而多時(shí)相InSAR除了能克服傳統(tǒng)InSAR技術(shù)的失相干問題外，也為系統(tǒng)誤差的建模和估計(jì)提供了可能[22].本文用附加系統(tǒng)參數(shù)對每景SAR影像的系統(tǒng)誤差進(jìn)行建模，對受湍流延遲和失相干噪聲影響的觀測值定權(quán)，用附加系統(tǒng)參數(shù)的加權(quán)最小二乘方法同時(shí)估計(jì)系統(tǒng)誤差和形變參數(shù).與現(xiàn)有方法相比，本文方法對研究區(qū)的形變類型不敏感，在長波形變時(shí)也能正確估計(jì)系統(tǒng)誤差和形變，增加了形變參數(shù)估計(jì)的精度和可靠性.

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

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

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

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

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

其中x，y表示P的空間位置，t為時(shí)間，aτ為形變模型的τ個(gè)待求參數(shù)向量.因變量x，y可選，當(dāng)形變在空間上的變化未知時(shí)，f（x，y，t|aτ）可簡化為f（t|aτ）.在tm時(shí)刻的主影像和ts時(shí)刻的從影像構(gòu)成的干涉圖中，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點(diǎn)的地形改正參數(shù)，Bperp，p為干涉圖的垂直基線，θ為入射角，R為主影像雷達(dá)天線與地面目標(biāo)之間的距離.式（1）中的φorb，p為軌道誤差相位、φatm，p為大氣延遲相位.φatm，p包含了地形相關(guān)的垂直分層延遲項(xiàng)φatm_ver，p

、長波延遲項(xiàng)φatm_long，p

和小尺度的湍流延遲項(xiàng)φatm_tur，p

[10]：

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

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

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

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

為避免高階多項(xiàng)式帶來的龍格現(xiàn)象[21]，式（6）最多取至二次項(xiàng)，令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主要包含湍流延遲和失相干噪聲的隨機(jī)誤差相位.

如果多時(shí)相InSAR包含由N景SAR影像組成的M幅干涉圖，同時(shí)對干涉圖中的S個(gè)目標(biāo)點(diǎn)進(jìn)行分析（如圖1所示），則共有M×S個(gè)觀測方程.寫成矩陣形式有：

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

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

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

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

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

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

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實(shí)例分析

2.1模擬實(shí)驗(yàn)

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

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

對應(yīng)的形變相位為

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

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

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

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

2.2真實(shí)數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

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

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

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

3討論

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

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

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

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

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

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

4結(jié)論

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

參考文獻(xià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].大地測量與地球動(dòng)力學(xué)，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

資助項(xiàng)目國家自然科學(xué)基金（41474007，41404013，41222027）

作者簡介韋建超，男，博士，研究方向?yàn)楹铣煽讖嚼走_(dá)干涉測量的誤差建模與參數(shù)反演.wjchao1608@163.com

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

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