孫建霖 楊忠杰 馬金輝

摘? 要：針對時序InSAR應(yīng)用過程中處理較為困難的大氣延遲誤差，文章利用StaMPS時序InSAR方法和TRAIN大氣校正工具箱，以劉家峽水庫及其以北地區(qū)為例，研究ERA-I、MERRA-2和GACOS三種常用的大氣數(shù)據(jù)在時序InSAR大氣延遲誤差校正中的差異。實驗表明，對輕度大氣影響區(qū)域，三種數(shù)據(jù)均可實現(xiàn)較好的校正效果;對中度大氣影響區(qū)域，三種數(shù)據(jù)的使用應(yīng)當區(qū)分不同的研究區(qū)域、數(shù)據(jù)覆蓋程度和大氣影響強弱而定，且三者之間并無明顯的優(yōu)劣之分;對重度大氣影響區(qū)域，在數(shù)據(jù)滿足基本條件的前提下，應(yīng)當優(yōu)先考慮剔除該影像。

關(guān)鍵詞：時序InSAR;大氣延遲;ERA-I;MERRA-2;GACOS

中圖分類號：V247.1+6? ? ? ? ?文獻標志碼：A? ? ? ? ?文章編號：2095-2945（2019）11-0023-05

Abstract： To deal with the difficult atmospheric delay error in the application of time-series InSAR， this paper uses StaMPS time series InSAR method and TRAIN atmospheric correction toolbox to study the difference of three commonly used atmospheric data ERA-I， MERRA-2 and GACOS. The example is Liujiaxia Reservoir and its northern area. The experiment shows that the three kinds of data can achieve good correction effect for the area with mild atmospheric influence. For the moderate atmospheric impact area， the three kinds of data should be used to distinguish different research areas， data coverage and atmospheric impact strength， and there is no obvious distinction between the three good and bad; for regions with severe atmospheric impacts， the removal of the image should be given priority on the premise that the data meets the basic conditions.

Keywords： Time Series InSAR; atmospheric delay; ERA-I; MERRA-2; GACOS

1 概述

對流層大氣延遲誤差校正是InSAR技術(shù)所面臨的最大挑戰(zhàn)之一[1]。其主要由兩次成像時對流層大氣壓力、溫度和相對濕度的不同，導(dǎo)致雷達信號穿過時的延遲量不等而引起。大氣延遲的存在可能掩蓋地表實際微小形變量，影響InSAR地表形變監(jiān)測精度[2]。因此，研究和選用合適的大氣誤差校正方法，成為InSAR應(yīng)用中的一項關(guān)鍵問題。對流層大氣延遲量，其主要包括三部分[2]：一是由于與近地表異常氣流有關(guān)的大氣湍流引起的分量;二是由大氣壓力、溫度和相對濕度變化引起的分量;三是由因地形起伏而引起的與地形有關(guān)的延遲分量。大氣延遲量可由輔助數(shù)據(jù)估算得到，如大氣模型數(shù)據(jù)、GPS測量數(shù)據(jù)、多光譜議等。

2001年，由Hanssen等[3]的研究表明，對特定高度h=h1和對流層頂部htop之間，沿雷達視線方向?qū)?yīng)干濕延遲大氣折射率為：

其中，P代表總大氣壓、T為溫度、e為水氣局部氣壓、？茲為入射角、？姿為雷達波長;系數(shù)k1、k2'、k3為統(tǒng)計或經(jīng)驗常數(shù)值，通常取k1=77.6KhPa-1、k2'=23.3KhPa-1、k3=3.75·105K2hPa-1[4]。SAR影像成像時大氣壓、溫度等數(shù)據(jù)可由輔助數(shù)據(jù)估算得到，由此即便可計算對應(yīng)大氣相位延遲。需要注意的是對InSAR干涉大氣延遲相位為主輔影像大氣延遲折射率差，而非整體折射率。本文選用三種較為容易獲得，且有質(zhì)量保證的大氣模型數(shù)據(jù)，對比分析基于不同數(shù)據(jù)的大氣延遲誤差校正效果，為時序InSAR應(yīng)用中的大氣校正研究者提供參考。

2 實驗區(qū)與數(shù)據(jù)處理

本文以甘肅省與青海省東南交接地區(qū)為實驗區(qū)，采用Sentinel-1A數(shù)據(jù)為例，使用DORIS（Delft object-oriented radar interferometric software）和StaMPS（Stanford Method for Persistent Scatterers）[5]開源軟件包完成相關(guān)干涉及時序處理。圖1為實驗區(qū)位置示意圖，長方形框為研究區(qū)地表覆蓋位置，大小約為150km×90km。三角為探空數(shù)據(jù)監(jiān)測站位置。

實驗區(qū)位于青海省和甘肅省交界帶，覆蓋區(qū)域主要包括臨夏回族自治州永靖縣、劉家峽水庫、蘭州河口鎮(zhèn)等地。實驗區(qū)主要受劉家峽水庫和黃河及其支流影響，大氣條件不均一且復(fù)雜多變，特別是在夏季，水庫蒸發(fā)量大，影響甚至控制區(qū)域大氣環(huán)境。實驗區(qū)最低點高程1545m，最高點高程2916m，高差1371m。

主影像依據(jù)Hooper等[5]提出的基于時空基線和多普勒基線的主影像選擇方法優(yōu)選20170320為主影像，共生成20景干涉圖，時間跨度約為2年。所用外部參考DEM數(shù)據(jù)為SRTM-3，水平分辨率90m。去除地平及地形相位影響，解纏并經(jīng)DEM誤差和軌道誤差校正后，在研究區(qū)提取共計7531632個PS點。由解纏相位可確定該區(qū)域夏季部分影像受大氣影響較重。限于本文篇幅和研究目的，僅給出4景特例，分別代表大氣影響較輕、輕、重、較重，對應(yīng)解纏并經(jīng)DEM和軌道誤差校正后的相位見圖2。4景干涉圖輔影像分別為20160206、20171022、20170823、20160816，在后文中分別對應(yīng)A、B、C、D，構(gòu)成4個對照組。大氣延遲相位計算使用由Bekaert等研究的TRAIN（Toolbox for Reducing Atmospheric InSAR Noise）[2]工具箱和式（1）、（2）完成。

3 基于ERA-I數(shù)據(jù)的大氣延遲相位校正

ERA-I（European Centre for Medium-Range Weather Forecasts Re-Analysis Interim）是由歐洲中期天氣預(yù)報中心（European Centre for Medium-Range Weather Forecasts，ECMWF）基于其2006年發(fā)布的IFS（Integrated Focecast System）circle31r2生產(chǎn)的自1989以來的中期氣象再分析資料[6]，目前提供近實時數(shù)據(jù)。數(shù)據(jù)垂直分層60層，水平分辨率約80km（采用T255格網(wǎng)），向上延伸至0.1hPa，采樣時間為每日0、6、12、18時（UTC時間）。觀測量主要包括大氣壓、濕度、溫度、風向等。經(jīng)實驗，位于本文研究區(qū)的ERA-I數(shù)據(jù)節(jié)點9個，為提高校正精度，將研究區(qū)進行擴展，實際使用數(shù)據(jù)節(jié)點共計30個。

利用由ERA-I估計得到的大氣延遲包括液態(tài)水延遲量和濕延遲量（湍流大氣延遲量）。圖3為A～D組總延遲量（第一行）和經(jīng)誤差校正后相位（第二行）。對大氣影響較輕的A組，在劉家峽水庫東側(cè)及西北側(cè)均存在不同程度的過校正問題。對大氣影響較重的D組，高海拔地區(qū)欠校正，而在沿河道的中海拔地區(qū)存在過校正問題。

4 基于MERRA/MERRA-2數(shù)據(jù)的大氣延遲相位校正

MERRA-2（the Modern-Era Retrospective Analysis for Research and Applications，version 2）[7]是由NASA Global Modeling and Assimilation Office （GMAO）為現(xiàn)代衛(wèi)星應(yīng)用生產(chǎn)的最新大氣再分析資料。相比其前身MERRA（2016年以前），MERRA-2除MERRA所有的數(shù)據(jù)外，還具有MERRA所不具備觀測數(shù)據(jù)類型，并且包括對GEOS（Goddard Earth Observing System）模型的數(shù)據(jù)支持。因此，本文所用數(shù)據(jù)為MERRA-2，時間為06時和12時數(shù)據(jù)，由NASA提供（https：//urs.earthdata.nasa.gov/home）。數(shù)據(jù)垂直分層為72層，水平分辨率為0.5°×0.625°，向上延伸至0.01hPa。經(jīng)實驗，位于本文研究區(qū)的MERRA-2數(shù)據(jù)節(jié)點有15個，同樣進行適當擴展（擴展范圍與ERA-I數(shù)據(jù)相同），共實際使用數(shù)據(jù)節(jié)點為42個。較ERA-I多12個，數(shù)據(jù)分辨率有所提高。

由MERRA-2估計得到的大氣延遲，與ERA-I類似，包括液態(tài)水延遲量和濕延遲量兩部分。圖4為基于MERRA-2數(shù)據(jù)計算得到的實驗區(qū)A～D組大氣總延遲量（第一行）和經(jīng)大氣誤差校正后相位（第二行）。在受大影響重和較重的C、D組校正相位中，由于校正誤差使位于劉家峽水庫東北側(cè)的已知沉降區(qū)完全淹沒在大氣相位中。但較ERA-I數(shù)據(jù)，MERRA-2在數(shù)據(jù)節(jié)點增多后，整體校正效果優(yōu)于ERA-I結(jié)果。因此，本文認為MERRA-2較ERA-I數(shù)據(jù)，在用于大氣延遲相位校正時大區(qū)域校正效果較好，但對地表的形變保相效果較差。

5 基于GACOS的大氣延遲相位校正

GACOS（Generic Atmospheric Correction Online Service）是由紐卡斯爾大學InSAR研究團隊提供的大氣改正數(shù)據(jù)[8]。其利用ITD（Iterative Tropospheric Decomposition）模型[9]分離對流層和湍流大氣延遲信號，并且提供用于InSAR大氣校正的高分辨率天頂總延遲圖。該數(shù)據(jù)使用高分辨率ECMWF天氣模型數(shù)據(jù)，水平格網(wǎng)大小為0.125°，垂直分層137層，6小時時間分辨率。表2為本實驗區(qū)所用GACOS數(shù)據(jù)信息。

GACOS所提供的僅為對應(yīng)日期研究區(qū)的大氣延遲量，而在InSAR中實際大氣延遲量為差分量，需在時間上相對于主影像，空間上相對于共同參考點做差分處理。在經(jīng)差分后，得到A～D組四個干涉像對對應(yīng)大氣延遲量和校正后相位，見圖5。其中第一行為大氣校正量，第二行為校正后相位。在大氣影響較輕的A、B組，可以分離出地表形變相位，但部分區(qū)域如劉家峽水庫西南區(qū)域。對C、D組，則完全不能分離出地表形變相位，且部分區(qū)域過校正度隨著大氣影響從ERA-I到MERRA-2隨著數(shù)據(jù)分辨率的提高，其大區(qū)域校正效果明顯提升。而且對較重大氣影響干涉圖的過校正問題同時得到緩解。GACOS數(shù)據(jù)較ERA-I和MERRA-2數(shù)據(jù)，在大氣影響較重的干涉圖中并未表現(xiàn)出好的校正效果，反而出現(xiàn)成片過校正區(qū)域。由于大氣校正其重點在大氣濕延遲校正，且無論何種方法或數(shù)據(jù)，都存在過校正或欠校正的問題。因此，在實際應(yīng)用中，在數(shù)據(jù)量能夠保證研究時間精度的條件下，應(yīng)當優(yōu)先考慮剔除大氣影像較重的干涉圖，僅保留大氣濕延遲未完全淹沒地表形變的干涉圖。在研究范圍較小時，ERA-I和MERRA-2數(shù)據(jù)空間分辨率會在很大程度上影響大氣校正參數(shù)估計，使校正產(chǎn)生偏差。而在研究區(qū)范圍過大時，大氣誤差又可能經(jīng)解纏傳播，影響全局。因此，對ERA-I、MERRA-2等大氣再分析數(shù)據(jù)，應(yīng)當依據(jù)實際研究范圍和地表形變特征選用。

選擇劉家峽水庫東北側(cè)一明顯沉降區(qū)，對比大氣校正后的時序形變結(jié)果。圖6為實驗區(qū)地表年均形變速率及劉家峽水庫東北側(cè)一沉降區(qū)局部?？梢?，基于不同大氣數(shù)據(jù)所得地表年均形變速率保持相對一致。圖7為圖6中放大局部對應(yīng)中心形變量曲線。

參考文獻：

[1]羅海濱，何秀鳳.用GPS改正InSAR大氣延遲誤差的研究[J].大地測量與地球動力學，2007，27（3）.

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

[3]Hanssen R F.Radar Interferometry[J].2001.

[4]Smith E K，Weintraub S.The Constants in the Equation for Atmospheric Refractive Index at Radio Frequencies[J].Proceedings of the Ire，1953，41（8）：1035-1037.

[5]Hooper A，Zebker H，Segall P，et al.A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers[J].Geophysical Research Letters，2004，31（23）：1-5.

[6]ECMWF.ERA-Interim[CP].https：//www.ecmwf.int/en/forecasts/datasets/archive-dataset

[7]Gelaro R，Mccarty W，Randles C，et al.The Modern-Era Retrospective Analysis for Research and Applications， Version 2 （MERRA-2）[J].Journal of Climate，2017，30（14）.

[8]Chen Y，Zhenhong L，Penna N T，et al.Generic atmospheric correction model for Interferometric Synthetic Aperture Radar observations[J].Journal of Geophysical Research：Solid Earth，2018.

[9]Yu C ，Penna N T，Li Z.Generation of real-time mode high-resolution water vapor fields from GPS observations[J].Journal of Geophysical Research：Atmospheres，2017，122（3）：2008-2025.