王輝 曾琪明 焦健 陳繼偉
結(jié)合序貫平差方法監(jiān)測(cè)地表形變的InSAR時(shí)序分析技術(shù)
王輝 曾琪明?焦健 陳繼偉
北京大學(xué)遙感與地理信息系統(tǒng)研究所, 北京大學(xué)地球與空間科學(xué)學(xué)院, 北京 100871; ?通信作者, E-mail: qmzeng@pku.edu.cn
基于短重訪周期 SAR 衛(wèi)星影像, 對(duì)黃河三角洲地表形變進(jìn)行高效和持續(xù)監(jiān)測(cè)的 SBAS-InSAR 時(shí)序分析。首先對(duì)研究區(qū)已有的 SAR 影像集進(jìn)行干涉處理, 得到干涉圖, 并進(jìn)行大氣效應(yīng)校正和軌道誤差去除, 然后利用傳統(tǒng)的 SBAS (small baseline subset)方法獲取地表形變。在此基礎(chǔ)上, 當(dāng)增加新的 SAR 數(shù)據(jù)時(shí), 采取漸進(jìn)式 SBAS 方法處理, 反演新時(shí)刻的地表形變。漸進(jìn)式 SBAS 方法融合序貫平差的思想, 以已有的解算結(jié)果為基礎(chǔ), 結(jié)合新的觀測(cè)數(shù)據(jù)進(jìn)行增量解算, 可以達(dá)到整體解算的效果。相對(duì)于傳統(tǒng)的 SBAS 方法, 每次增加新影像都要采用重新全部解算的方式, 能夠減少冗余運(yùn)算, 提高計(jì)算效率。實(shí)驗(yàn)證明, 基于 2018 年 5 月—2016 年 8 月在黃河三角洲地區(qū)獲取的 Sentinel-1 衛(wèi)星 SAR 數(shù)據(jù), 利用漸進(jìn)式 SBAS 方法反演的地表形變與地表實(shí)測(cè)控制點(diǎn)結(jié)果相近, 相關(guān)系數(shù)()為 0.82, 且與傳統(tǒng)的 SBAS 方法反演得到的地表形變速率差異在 1mm/a內(nèi), 解算時(shí)間縮短約 40%, 能夠持續(xù)高效地監(jiān)測(cè)地表形變。
InSAR時(shí)序分析; 小基線集(SBAS); 漸進(jìn)式SBAS; 序貫平差; 地表形變; 持續(xù)監(jiān)測(cè)
雷達(dá)差分干涉測(cè)量(differential interferometry synthetic aperture radar, D-InSAR)是通過(guò)對(duì)兩景或多景 SAR 影像進(jìn)行干涉處理, 獲取高精度和高空間分辨率微小形變的技術(shù)[1-3]。傳統(tǒng)的 D-InSAR 技術(shù)因時(shí)空失相干和大氣效應(yīng)等因素的干擾[4-5], 測(cè)量精度受到影響。SBAS (small baseline set, 小基線集)-InSAR 時(shí)序分析技術(shù)[6-8]通過(guò)對(duì)時(shí)域上一系列差分干涉影像的處理, 能夠消除或減弱上述因素的影響, 得到毫米級(jí)的沉降速率。目前, SBAS-InSAR已廣泛應(yīng)用于對(duì)地面沉降、滑坡和泥石流等災(zāi)害的研究中[9-11], 成為長(zhǎng)時(shí)間緩慢地表形變監(jiān)測(cè)的主流方法之一。
隨著測(cè)量技術(shù)的發(fā)展, 以歐洲航天局(European Space Agency, ESA) Sentinel-1 為代表的 SAR 衛(wèi)星重訪周期不斷縮短, 可以穩(wěn)定地獲取 SAR 數(shù)據(jù), 為InSAR 時(shí)序分析提供豐富的數(shù)據(jù)基礎(chǔ), 也給高效、持續(xù)地反演地表形變帶來(lái)挑戰(zhàn)。傳統(tǒng)的 SBAS 方法通常是在某區(qū)域積累一定量的數(shù)據(jù)后進(jìn)行處理, 獲得數(shù)據(jù)積累時(shí)間段內(nèi)的形變監(jiān)測(cè)結(jié)果[12-13]; 若新增幾景 SAR 影像, 就要“從頭開(kāi)始”, 重復(fù)全部的處理流程來(lái)獲取新時(shí)刻的形變結(jié)果。隨著區(qū)域 SAR 數(shù)據(jù)的增加, 運(yùn)算的復(fù)雜度也快速加大。因此, 這種“從頭開(kāi)始”的處理思想難以滿足動(dòng)態(tài)、高效和持續(xù)的地表形變監(jiān)測(cè)的需求。
陳繼偉[14]將序貫平差思想[15]融合到 InSAR 時(shí)序分析中, 并探索在已有的數(shù)據(jù)集中新增一景 SAR影像后反演新時(shí)刻地表沉降量的可能性, 在保證精度的同時(shí), 提高計(jì)算效率。Ansari 等[16]提出結(jié)合序貫估計(jì)的高效 InSAR 時(shí)序算法, 采用遞歸估計(jì)和分析數(shù)據(jù)協(xié)方差矩陣, 將數(shù)據(jù)劃分為小數(shù)據(jù)集, 然后對(duì)數(shù)據(jù)集進(jìn)行壓縮, 避免對(duì)整個(gè)數(shù)據(jù)集進(jìn)行重新處理, 能夠極大地提高時(shí)間效率。胡俊等[17]研究新增兩張或多景 SAR 影像后, 利用結(jié)合序貫平差方法的InSAR 時(shí)序分析技術(shù)反演地表形變的一般性規(guī)律, 并用模擬數(shù)據(jù)驗(yàn)證方法的可行性。
本研究提出漸進(jìn)式 SBAS 時(shí)序分析方法融合序貫平差的思想, 在已有解算結(jié)果的基礎(chǔ)上, 結(jié)合新的觀測(cè)數(shù)據(jù)進(jìn)行增量解算, 便可達(dá)到整體解算的效果, 實(shí)現(xiàn)對(duì)地表形變的高效和持續(xù)監(jiān)測(cè)。選取 2015年 5 月— 2016 年 8 月覆蓋黃河三角洲地區(qū)的 34 景Sentinal-1 衛(wèi)星 SAR 影像作為實(shí)驗(yàn)數(shù)據(jù), 開(kāi)展?jié)u進(jìn)式SBAS 時(shí)序分析。通過(guò)與實(shí)測(cè)控制點(diǎn)數(shù)據(jù)以及傳統(tǒng)SBAS 時(shí)序分析方法反演的地表形變對(duì)比, 驗(yàn)證漸進(jìn)式 SBAS 方法反演地表形變的可靠性和在時(shí)間效率上的優(yōu)越性。
如圖 1 所示, 漸進(jìn)式 SBAS 方法分為兩個(gè)階段。階段 1: 已有 SAR 數(shù)據(jù)的傳統(tǒng) SBAS 方法求解; 階段 2: 新增 SAR 數(shù)據(jù)后的漸進(jìn)式 SBAS 方法求解。階段 1 和 2 各自包含兩個(gè)步驟, 其中步驟 1 是干涉圖的大氣效應(yīng)校正[18]和軌道誤差去除[19], 兩個(gè)階段的步驟 1 相似, 僅處理的干涉圖不相同。階段 1 的SBAS 求解(步驟 2)為已有 SAR 數(shù)據(jù)對(duì)應(yīng)的地表形變量解算, 為階段 2 的漸進(jìn)式 SBAS 求解(步驟 2)提供數(shù)據(jù)基礎(chǔ), 使得階段 2 通過(guò)增量結(jié)算獲得新時(shí)刻的地表形變, 同時(shí)更新已有的地表形變量。
圖1 漸進(jìn)式SBAS方法流程
衛(wèi)星對(duì)同一區(qū)域連續(xù)觀測(cè), 在0–t時(shí)間內(nèi)獲得+1 景 SAR 影像, 在給定時(shí)間基線閾值t和空間基線閾值b后, 得到張干涉圖。將經(jīng)過(guò)大氣效應(yīng)校正和軌道誤差去除的差分干涉圖作為觀測(cè)值, 地表形變量作為待求值進(jìn)行 SBAS 求解, 得到的觀測(cè)方程為
將式(1)和(2)整理成如下矩陣形式:
亦即
為了驗(yàn)證漸進(jìn)式 SBAS 方法反演地表形變的精度和時(shí)間效率, 選擇實(shí)驗(yàn)區(qū)內(nèi)包含一定量數(shù)據(jù)的SAR 數(shù)據(jù)集, 分別用漸進(jìn)式 SBAS 方法和傳統(tǒng)SBAS 方法處理, 進(jìn)行 3 個(gè)方面的對(duì)比: 1)對(duì)比漸進(jìn)式 SBAS 方法反演的地表形變結(jié)果與實(shí)測(cè)控制點(diǎn)測(cè)量結(jié)果; 2)對(duì)比兩種方法的反演結(jié)果; 3)統(tǒng)計(jì)兩種方法反演地表形變所需時(shí)間。
選取山東半島東營(yíng)市境內(nèi)黃河三角洲為研究區(qū)。如圖 2 所示, 該區(qū)域是以墾利縣寧海為頂點(diǎn), 向東展開(kāi)的扇形區(qū)域[21], 受地質(zhì)構(gòu)造及人類(lèi)活動(dòng)的影響, 存在明顯的地面持續(xù)沉降。已采取多種手段對(duì)該區(qū)域的地面沉降進(jìn)行監(jiān)測(cè), 以傳統(tǒng)的水準(zhǔn)測(cè)量[22]和 GPS 監(jiān)測(cè)[23]為主。
選取 34 景 Sentinel-1 衛(wèi)星 SAR 數(shù)據(jù), 影像的時(shí)間跨度為 2015 年 5 月 19 日—2016 年 8 月 29 日。采用覆蓋實(shí)驗(yàn)區(qū) 30m 分辨率的 SRTM (shuttle radar topography mission) DEM 作為參考高程數(shù)據(jù), 輔助去除地形相位。同時(shí)獲取東營(yíng)市雙王城水庫(kù) 2015年 9 月底首測(cè)和 2016 年 8 月復(fù)測(cè)的 16 個(gè)控制點(diǎn), 用于對(duì)比驗(yàn)證。
2.3.1 干涉網(wǎng)絡(luò)設(shè)計(jì)及干涉處理
以 2015 年 11 月 15 日的圖像作為主圖像, 時(shí)間基線閾值為 72 天, 空間基線閾值為 100m, 得到 115組干涉影像對(duì), 構(gòu)建時(shí)空基線網(wǎng)絡(luò)(圖 3)。用基于地理位置的方法進(jìn)行配準(zhǔn), 在重采樣后, 用 GMTSAR軟件進(jìn)行干涉處理。結(jié)合外部 DEM 去除地形相位, 并進(jìn)行相位解纏, 最后得到去除大氣效應(yīng)和軌道誤差的干涉圖。此外, 為了抑制干涉圖噪聲, 進(jìn)行多視處理(方位向 4 視, 距離向 12 視)。
圖2 研究區(qū)地理位置
圖3 干涉對(duì)的時(shí)空基線網(wǎng)絡(luò)
2.3.2 漸進(jìn)式SBAS和傳統(tǒng)SBAS處理
已有 SAR 數(shù)據(jù)的傳統(tǒng) SBAS 處理(圖 4): 2015年 5 月 19 日―2016 年 6 月 6 日共 30 景圖像, 生成的106 對(duì)干涉對(duì)(圖 3 中時(shí)空基線網(wǎng)絡(luò)黑色部分)作為已有SAR 數(shù)據(jù)集, 反演已有 SAR 數(shù)據(jù)的對(duì)應(yīng)的地表形變量。
新增 SAR 數(shù)據(jù)的傳統(tǒng) SBAS 方法處理(圖 4): 在 2015 年 5 月 19 日―2016 年 6 月 6 日共 30 景 SAR影像的基礎(chǔ)上, 新增 1~4 景 SAR 影像后, 利用傳統(tǒng)SBAS 方法, 對(duì) 2015 年 5 月 19 日―2016 年 8 月 29 日所有 SAR 影像重新整體解算, 獲取新時(shí)刻對(duì)應(yīng)的地表形變量。
新增 SAR 數(shù)據(jù)后漸進(jìn)式 SBAS 方法處理(圖 4):以第一部分傳統(tǒng) SBAS 方法的產(chǎn)生中間結(jié)果為基礎(chǔ)(虛線框, 不再處理), 使用 2015 年 6 月 30 日―2016年 8 月 29 日時(shí)間段內(nèi) 1–4 景 SAR 影像(實(shí)線框), 結(jié)合漸進(jìn)式 SBAS 方法進(jìn)行增量處理, 反演新增影像對(duì)應(yīng)時(shí)刻的地表形變量。
圖 5 為兩種方法反演得到的黃河三角洲地區(qū)2015 年 5 月 19 日—2016 年 8 月 29 日沿視線方向的形變速率場(chǎng), 可以看出, 該區(qū)域存在大范圍的地表形變, 空間分布不均勻, 東營(yíng)市、濱州市和廣饒縣等主要城鎮(zhèn)存在不同程度的沉降, 區(qū)域東部存在明顯的沉降漏斗(圖 5 中藍(lán)-紫色區(qū)域)。
選取雙王城水庫(kù) 2015 年 9 月底首測(cè), 2016 年 8月復(fù)測(cè)的 16 個(gè)控制點(diǎn)(圖 6)作為驗(yàn)證數(shù)據(jù), 其高程坐標(biāo)由二等水準(zhǔn)測(cè)量獲得。InSAR 測(cè)量結(jié)果為 LOS方向, 地面控制點(diǎn)獲得的是水平和高程的形變。為了統(tǒng)一形變測(cè)量方向, 將控制點(diǎn)形變量投影至 LOS方向。
圖4 新增影像后漸進(jìn)式SBAS和傳統(tǒng)SBAS獲取新時(shí)刻地表形變的實(shí)驗(yàn)設(shè)計(jì)
(a)傳統(tǒng)SBAS; (b)漸進(jìn)式SBAS
用漸進(jìn)式 SBAS 得到的形變速率與地面測(cè)量獲得的形變速率的相關(guān)系數(shù)()為 0.82, 差值絕對(duì)值最大為 8.6mm/a, 最小為 0.48mm/a。兩種方法得到的結(jié)果顯示水庫(kù)東西兩側(cè)形變規(guī)律有較高的一致性, 形變速率雖然有差異, 但處在同一個(gè)數(shù)量級(jí)(圖 7)。
對(duì)兩種方法反演得到的整個(gè)研究區(qū)形變速率(2015 年 5 月 19 日— 2016 年 8 月 29 日)進(jìn)行差值運(yùn)算, 結(jié)果如圖 8(a)所示。對(duì)形變速率的差值進(jìn)行統(tǒng)計(jì), 結(jié)果如下: 絕大多數(shù)差值<0.3mm/a, 平均值為0.08mm/a, 標(biāo)準(zhǔn)差為 0.168mm/a, 兩種方法反演得出的地表形變速率差異極小。
圖6 地面控制點(diǎn)分布
圖7 漸進(jìn)式SBAS反演形變速率與地面實(shí)測(cè)形變速率對(duì)比
選取研究區(qū)東部和南部沉降漏斗區(qū), 更詳細(xì)地比較兩種方法反演典型區(qū)域地表形變的差異, 結(jié)果如圖 8(b)和(c)所示。在沉降漏斗附近, 兩種方法解算出的累積地表形變量幾乎相同, 絕大多數(shù)差值接近 0, 只有極小的區(qū)域存在差值為 4~8mm 的異常點(diǎn)(可能是系統(tǒng)誤差導(dǎo)致), 與沉降漏斗幾百毫米的沉降量相比, 可以忽略不計(jì)。上述結(jié)果說(shuō)明, 即使在沉降漏斗附近, 兩種方法的反演結(jié)果也十分接近。
進(jìn)一步地, 在研究區(qū)東部、南部和東南部嚴(yán)重沉降區(qū)域(見(jiàn)圖 5), 隨機(jī)選取 4 個(gè)特征點(diǎn)位, 對(duì)比兩種方法反演的沿視線方向形變量, 結(jié)果如圖 9 所示??傮w而言, 所選特征點(diǎn)位的累計(jì)形變量差異在 1.5 mm 以內(nèi), 與累積形變量相比, 這樣的差異顯得微小, 可見(jiàn)兩種方法反演的累計(jì)形變量基本上相同; 在像元層次上, 兩種方法反演的地表形變也幾乎相同。
漸進(jìn)式 SBAS 方法是為適應(yīng)動(dòng)態(tài)增加的 SAR 數(shù)據(jù)的時(shí)序處理而提出的, 每次新增觀測(cè)數(shù)據(jù), 都可獲得該數(shù)據(jù)成像時(shí)刻的累計(jì)形變量。將新增 2~3 景SAR 影像時(shí)漸進(jìn)式 SBAS 反演的累計(jì)形變量與傳統(tǒng)SBAS 方法進(jìn)行對(duì)比, 結(jié)果見(jiàn)圖 10??傮w而言, 兩種方法反演地表累計(jì)形變量的差異均值和標(biāo)準(zhǔn)差均很小, 在 1mm 以內(nèi), 說(shuō)明漸進(jìn)式 SBAS 在處理連續(xù)的動(dòng)態(tài)的 SAR 數(shù)據(jù)時(shí), 仍與傳統(tǒng) SBAS 方法保持極高的一致性。
傳統(tǒng) SBAS 和漸進(jìn)式 SBAS 是按照相同的規(guī)則選擇高相干點(diǎn), 所以兩種方法處理的高相干點(diǎn)個(gè)數(shù)相等。為了將問(wèn)題簡(jiǎn)化, 我們以其中任意一個(gè)高相干點(diǎn)為例, 分析大氣效應(yīng)和軌道誤差糾正(步驟1)以及形變量求解(步驟 2)的時(shí)間復(fù)雜度。
(a)整個(gè)研究區(qū); (b)東部沉降漏斗; (c)南部沉降漏斗
(a1)~(a4) 傳統(tǒng)SBAS; (b1)~(b4) 兩種方法的差值
假定實(shí)驗(yàn)數(shù)據(jù)集合由景 SAR 影像構(gòu)成的已有 SAR 數(shù)據(jù)集和′景 SAR 影像構(gòu)成的新增 SAR 數(shù)據(jù)集組成, 在相同的時(shí)空基線約束下, 分別形成和′張干涉圖。對(duì)于單個(gè)高相干點(diǎn), 新增′景SAR 影像后, 傳統(tǒng) SBAS 在步驟 1 需要對(duì)(+′)張干涉圖進(jìn)行大氣效應(yīng)和軌道誤差糾正, 漸進(jìn)式SBAS 則只需對(duì)新增 SAR 數(shù)據(jù)和部分已有 SAR 數(shù)據(jù)形成的′張干涉圖進(jìn)行大氣效應(yīng)和軌道誤差糾正。每次新增 SAR 數(shù)據(jù), 傳統(tǒng) SBAS 都避不開(kāi)對(duì)已有 SAR 數(shù)據(jù)集的處理, 因此比漸進(jìn)式 SBAS 效率低很多。兩種方法處理問(wèn)題的規(guī)模與干涉圖的個(gè)數(shù)有關(guān), 時(shí)間復(fù)雜度為((+))和()。從理論上講, 當(dāng)新增 SAR 影像數(shù)不多時(shí), 已有 SAR 數(shù)據(jù)形成的干涉圖數(shù)量遠(yuǎn)大于新增 SAR 數(shù)據(jù)與部分已有 SAR 數(shù)據(jù)形成的干涉圖′, 因此漸進(jìn)式 SBAS 比傳統(tǒng) SBAS 效率高很多。在形變量求解階段(步驟 2), 傳統(tǒng) SBAS 需要求解(+)景 SAR 數(shù)據(jù)對(duì)應(yīng)的地表形變, 存在大量對(duì)已有景 SAR 數(shù)據(jù)集對(duì)應(yīng)的地表形變的重復(fù)運(yùn)算, 冗余度高, 漸進(jìn)式 SBAS 則只針對(duì)新增的′景 SAR 數(shù)據(jù)對(duì)應(yīng)的地表形變求解, 并更新已有的景 SAR 數(shù)據(jù)對(duì)應(yīng)的地表形變, 計(jì)算的復(fù)雜度大大降低。從 1.2 節(jié)的介紹可知, 形變量的求解主要涉及線性方程組(式(5)和(6))的最小二乘法計(jì)算。線性方程組求解的復(fù)雜度與待求解未知數(shù)個(gè)數(shù)有關(guān), 兩種方法待求解的未知數(shù)個(gè)數(shù)分別為(+)和, 時(shí)間復(fù)雜度分別為((+)3)和(3)。
在配置相同的計(jì)算機(jī)上分別運(yùn)行漸進(jìn)式 SBAS與傳統(tǒng) SBAS, 計(jì)算機(jī)配置為 Intel Core i5 CPU (2.8GHz), 4G 內(nèi)存。針對(duì)每次新增的 SAR 圖像, 分別記錄大氣效應(yīng)和軌道誤差糾正以及形變量求解所需的時(shí)間。
圖10 新增2~3景SAR影像時(shí)漸進(jìn)式SBAS (a)與傳統(tǒng)SBAS (b)反演累計(jì)形變量的差值
相鄰柱子中, 左側(cè)為傳統(tǒng)SBAS, 右側(cè)為漸進(jìn)式SBAS
圖 11 比較在 30 景 SAR 影像的基礎(chǔ)上, 新增 1~ 4 景 SAR 影像后兩種方法反演地表形變所需時(shí)間。在步驟 1, 傳統(tǒng) SBAS 所需時(shí)間是漸進(jìn)式 SBAS 的 3~ 5 倍, 原因是漸進(jìn)式 SBAS 只處理新增的干涉圖, 而傳統(tǒng) SBAS 對(duì)所有干涉圖都進(jìn)行處理, 因此漸進(jìn)式SBAS 能大大地減少運(yùn)行的時(shí)間, 與前面的理論分析可以較好地相互驗(yàn)證。在形變量求解階段, 雖然從理論上講漸進(jìn)式 SBAS 的時(shí)間效率應(yīng)該比傳統(tǒng)SBAS 高很多, 但統(tǒng)計(jì)結(jié)果表明兩者在前期(新增 1~ 2 景 SAR 圖像時(shí))的運(yùn)行時(shí)間差異并不大, 原因主要是在解算形變結(jié)果時(shí), 漸進(jìn)式 SBAS 一方面要引入和使用已有數(shù)據(jù)的傳統(tǒng) SBAS 求解階段保留的協(xié)方差矩陣, 另一方面還要對(duì)已有數(shù)據(jù)反演的地表形變結(jié)果進(jìn)行修正, 此外還有干涉相位相關(guān)文件的讀入與保存, 故不能像前兩個(gè)階段那樣明顯地減少運(yùn)行時(shí)間。但是, 隨著 SAR 圖像增加, 漸進(jìn)式 SBAS 的時(shí)間效率相對(duì)于傳統(tǒng) SBAS 仍得到較大的提升, 在增加第 4 景 SAR 圖像時(shí), 其形變求解時(shí)間是傳統(tǒng) SBAS 的 64.7%。
總體而言, 在 30 景 SAR 影像的基礎(chǔ)上新增加1~4 景 SAR 影像時(shí), 傳統(tǒng) SBAS 的運(yùn)行時(shí)間隨著SAR 圖像增加呈現(xiàn)線性增加趨勢(shì), 與傳統(tǒng) SBAS 相比, 漸進(jìn)式 SBAS 的時(shí)間效率不斷提高。因此, 處理動(dòng)態(tài)增加的 SAR 數(shù)據(jù)時(shí), 漸進(jìn)式 SBAS 的時(shí)間效率更高。
本文以 Sentinel-1 衛(wèi)星 SAR 影像作為實(shí)驗(yàn)數(shù)據(jù), 使用漸進(jìn)式 SBAS 方法, 獲得黃河三角洲地區(qū)地表形變信息。通過(guò)對(duì)比傳統(tǒng) SBAS 方法的處理結(jié)果與地面實(shí)測(cè)控制點(diǎn)數(shù)據(jù), 驗(yàn)證了漸進(jìn)式 SBAS 方法的可靠性。結(jié)合理論分析和定量統(tǒng)計(jì), 證明了漸進(jìn)式SBAS 方法在時(shí)間上的高效性。
使用雙王城水庫(kù)地面實(shí)測(cè)控制點(diǎn)數(shù)據(jù), 對(duì)漸進(jìn)式 SBAS 方法反演地表形變結(jié)果進(jìn)行驗(yàn)證, 兩者形變規(guī)律一致, 說(shuō)明漸進(jìn)式 SBAS 方法反演地表形變的精度是可靠的。
漸進(jìn)式 SBAS 與傳統(tǒng) SBAS 方法反演得到的地表形變速率差值在 1mm/a 以內(nèi), 累計(jì)形變量差值在1.5mm 以內(nèi), 得到的形變監(jiān)測(cè)結(jié)果基本上相同。
與傳統(tǒng) SBAS 方法相比, 針對(duì)隨時(shí)間動(dòng)態(tài)增加的數(shù)據(jù), 漸進(jìn)式 SBAS 方法具有更高的時(shí)間效率, 適用于地表形變的連續(xù)、動(dòng)態(tài)和高效監(jiān)測(cè)。
[1] Massonnet D, Rossi M, Carmona C, et al. The dis-placement field of the Landers earthquake mapped by radar interferometry. Nature, 1993, 364: 138?142
[2] Massonnet D, Thatcher W, Vadon H. Detection of postseismic fault-zone collapse following the Landers earthquake. Nature, 1996, 382: 612?616
[3] Zebker H A, Rosen P. On the derivation of coseismic displacement fields using differential radar interfero-metry: the Landers earthquake. Geosci Remote Sens Symp, 1994, 32(5): A206?A207
[4] 董玉森, Ge Linlin, Chang Hsingchun, 等. 基于差分雷達(dá)干涉測(cè)量的礦區(qū)地面沉降監(jiān)測(cè)研究. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2007, 32(10): 888?891
[5] Zebker H A, Villasenor J. Decorrelation in interfero-metric radar echoes. IEEE Transactions on Geosci-ence and Remote Sensing, 1992, 30(5): 950?959
[6] Berardino P, Fornaro G, Lanari R, et al. A new algo-rithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(11): 2375?2383
[7] Hooper A, Zebker H, Segall P, et al. A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers. Geophy-sical Research Letters, 2004, 31(23): 1?5
[8] Hooper A, Segall P, Zebker H. Persistent scatterer inteferometric synthetic aperture radar for crustal deformation analysis, with application to Vlocan Alcedo, Galapagos. Journal of Geophysical Research, 2007, 112 (B7): Doi: 10. 1029/2006JB004763
[9] 陳繼偉, 曾琪明, 焦健, 等. Sentinel-1A衛(wèi)星TOPS模式數(shù)據(jù)的SBAS時(shí)序分析方法——以黃河三角洲地區(qū)為例. 國(guó)土資源遙感, 2017, 29(4): 82?87
[10] 張金芝, 黃海軍, 畢海波, 等. SBAS時(shí)序分析技術(shù)監(jiān)測(cè)現(xiàn)代黃河三角洲地面沉降. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2016, 41(2): 242?248
[11] 李珊珊, 李志偉, 胡俊, 等. SBAS-InSAR技術(shù)監(jiān)測(cè)青藏高原季節(jié)性凍土形變. 地球物理學(xué)報(bào), 2013, 56(5): 1476?1486
[12] Casu F, Manzo M, Lanari R. A quantitative assess-ment of the SBAS algorithm performance for surface deformation retrieval from D-InSAR data. Remote Sensing of Environment, 2006, 102(3): 195?210
[13] 楊成生, 劉媛媛, 敖萌. 基于SBAS時(shí)序分析的大同地面沉降與地下水活動(dòng)研究. 國(guó)土資源遙感, 2015, 27(1): 127?132
[14] 陳繼偉. TOPS模式SAR數(shù)據(jù)的漸進(jìn)式SBAS時(shí)序分析方法[D]. 北京: 北京大學(xué), 2017
[15] 曾安敏, 楊元喜, 歐陽(yáng)桂崇. 附加約束條件的序貫平差. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2008, 33(2): 183?186
[16] Ansari H, de Zan F, Bamler R. Sequential estimator: toward efficient InSAR time series analysis. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(10): 5637?5652
[17] 胡俊, 劉計(jì)洪, 李志偉, 等. 一種基于序貫平差的InSAR地表形變監(jiān)測(cè)方法[P]. 中國(guó) 109061641A, 2018?12?21
[18] Jolivet R, Grandin R, Lasserre C, et al. Systematic InSAR tropospheric phase delay corrections from global meteorological reanalysis data. Geophysical Research Letters, 2011, 38(17): Doi: 10.1029/2011GL048757
[19] Cavalié O, Lasserre O, Doin M P, et al. Measurement of interseismic strain across the Haiyuan fault (Gansu, China), by InSAR. Earth and Planetary Science Let-ters, 2008, 275(3): 246?257
[20] 隋立芬, 劉雁雨, 王威. 自適應(yīng)序貫平差及其應(yīng)用. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2007, 32(1): 51?54
[21] 葉青超. 黃河三角洲的地貌結(jié)構(gòu)及發(fā)育模式. 地理學(xué)報(bào), 1982, 37(4): 349?363
[22] 宋波, 王德生, 王錦麗. 東營(yíng)地面沉降監(jiān)測(cè). 地礦測(cè)繪, 2004, 20(1): 34?36
[23] 毛繼軍, 孟黎, 蘇艷紅, 等. 基于GPS的東營(yíng)市地面沉降監(jiān)測(cè)研究. 工程勘察, 2014, 42(10): 56?59
InSAR Time Series Analysis Technique Combined with Sequential Adjustment Method for Monitoring of Surface Deformation
WANG Hui, ZENG Qiming?, JIAO Jian, CHEN Jiwei
Institute of Remote Sensing and Geographical Information System, School of Earth and Space Sciences, Peking University, Beijing 100871; ? Corresponding author, E-mail: qmzeng@pku.edu.cn
Based on the SAR satellite imagery with short revisit periods, this paper proposes a novel SBAS-InSAR time series analysis technique for efficient and continuous monitoring of surface deformation in the study area. First, standard interferometric processing is carryed out with the existing SAR image set in the study area to obtain the interference map. Thenthe atmospheric and orbital errors are removed, and the traditional SBAS (small baseline subset) method is used to obtain the surface deformation. On this basis, when a new SAR image is added, the progressive SBAS is adopted to invert the surface deformation at the new moment. The progressive SBAS method integrates the idea of ??sequential adjustment based on the obtained results derived from existing data set, and combines the newly acquired data to implement incremental calculations, finally achieves the equivalent effect of overall processing. Compared with the traditional SBAS method which needs to resolve all the calculations every time when a new image is added, the progressive SBAS method can reduce redundant operations and improve computing efficiency. The experiment proves that based on the Sentinel-1 satellite SAR data acquired in the Yellow River Delta from May 2018 to August 2016, the surface deformations retrieved by the progressive SBAS method are almost the same as the results of the measured ground level. The correlation coefficient () is 0.82, and the difference between the ground deformation rate and the traditional SBAS method is within 1 mm/a. The solution time is shortened by about 40%, and the ground deformation can be efficiently and continuously monitored.
InSAR time series analysis; small baseline set (SBAS); progressive SBAS; sequential adjustment; surface deformation; continuous monitoring
10.13209/j.0479-8023.2021.002
國(guó)家自然科學(xué)基金(41571337)資助
2020–02–24;
2020–05–28