焦廣棋，孫 玉,楊 琰

(1. 福州大學 空間數(shù)據(jù)挖掘與信息共享教育部重點實驗室， 福建 福州 350108； 2. 地理空間信息技術國家地方聯(lián)合工程研究中心， 福建 福州 350108；3. 衛(wèi)星空間信息技術綜合應用國家地方聯(lián)合工程研究中心， 福建 福州 350108)

星載合成孔徑雷達(synthetic aperture radar，SAR)信號穿透大氣層時受到壓強、溫度和總電子含量等影響，其傳播速率及路徑發(fā)生改變而導致雷達信號相位延遲，稱為大氣延遲效應[1]。大氣延遲包含電離層延遲和對流層延遲：電離層延遲由自由電子密度分布不均引起，可顯著影響極地、高緯地區(qū)和長波段(L波)信號，但對C波段信號影響較??；對流層中的溫度、濕度和水汽等變化引起對流層延遲，該延遲可影響任何波段和地區(qū)，水汽含量豐富的地區(qū)影響更大[2]。本研究僅討論對流層延遲對C波段信號的影響，對電離層延遲未做修正。因此，所使用的Sentinel 1A影像干涉原始相位中包含不同程度的對流層延遲，嚴重干擾了形變監(jiān)測及其機制的正確解譯，是合成孔徑雷達干涉測量技術(interferometric SAR, InSAR)數(shù)據(jù)處理中主要的誤差源之一[3]。因此，估計并改正對流層延遲對基于InSAR技術的高精度形變監(jiān)測至關重要。

對流層延遲可分為由地形起伏引起的垂直分層延遲和具有隨機性的大氣對流過程引起的湍流混合延遲。目前，國內(nèi)外學者已針對對流層延遲改正開展了諸多研究，主要分為3類：①時空濾波法。整個大氣延遲被假設為具有時空相關的隨機分布性質，時空濾波可以有效抑制湍流混合延遲，但選擇合適的濾波器和濾波窗口大小成為難點。時空濾波常被用于時序InSAR技術處理流程中去除對流層延遲，比如永久散射體方法(persistent scatterer InSAR，PSInSAR)、短基線子集方法(small baseline subset，SBAS)[4]。②函數(shù)模型。以干涉相位與地形高程之間的關聯(lián)性建立函數(shù)模型，包括線性模型和冪律模型，能夠有效去除垂直分層延遲[5]。③外部輔助數(shù)據(jù)集。利用外部數(shù)據(jù)集提供氣象參數(shù)計算改正對流層延遲，但光學遙感數(shù)據(jù)受限于云、雨天氣及夜間環(huán)境；全球導航定位衛(wèi)星系統(tǒng)(global navigation satellites system，GNSS)數(shù)據(jù)受GNSS站點密度的影響導致改正效果不佳；全球大氣模型(global atmospheric models，GAM)是近年來發(fā)展的改正對流層延遲的主要數(shù)據(jù)源，以全球覆蓋、數(shù)據(jù)免費、適用性強等優(yōu)勢在單幅、時序干涉圖對流層延遲改正中應用廣泛，包括現(xiàn)代衛(wèi)星應用生產(chǎn)再分析資料(Modern-Era retrospective analysis for research and applications, v2，MERRA2)、歐洲中期天氣預報中心(European centre for medium-range weather forecasts，ECMWF)發(fā)布的再分析資料(ECMWF reanalysis interim/v5，ERA-Interim/5)[6]。InSAR通用性大氣改正在線服務(generic atmospheric correction online service，GACOS)融合了ECMWF數(shù)據(jù)集和連續(xù)全球定位系統(tǒng)(global positioning system，GPS)對流層延遲估計數(shù)據(jù)，采用迭代對流層分解模型處理后向用戶免費提供全球性的對流層延遲數(shù)據(jù)[7]。

針對GAM以及GACOS改正時序InSAR干涉圖對流層延遲的可靠性，李鵬等[1]采用線性模型、GACOS及ERA-Interim改正覆蓋阿爾金斷裂帶地區(qū)Sentinel 1A和ASAR影像生成干涉圖中的對流層延遲，結果表明兩類SAR影像生成的干涉圖經(jīng)GACOS改正后的相位標準差削減量分別可達28.8%和14.8%，當去除軌道平面后可達54.5%和68.1%。Wang等[3]以青藏高原地區(qū)的寬幅Sentinel 1A影像為例，對比評估ERA-Interim、ERA5、GACOS及MERRA2對流層延遲改正的效果，結果表明ERA5與GACOS具有相近的改正效果，驗證了2018年新發(fā)布時間分辨率最高的ERA5應用于青藏高原InSAR對流層延遲改正的巨大潛力。ERA5改正對流層延遲的效果在內(nèi)陸地區(qū)已被證明，但在沿海地勢平坦地區(qū)相比其他GAM的表現(xiàn)尚需進行對比研究。日本關東平原位于沿海地勢平坦、地質構造復雜地區(qū)，2011年Mw9.0地震引起地表滲透率改變導致關東平原內(nèi)出現(xiàn)地表隆升現(xiàn)象，更需進一步開展該地區(qū)的地殼形變監(jiān)測研究[8]。

不同對流層延遲改正數(shù)據(jù)集的改正效果具有差異性，需要可靠的統(tǒng)計指標去衡量改正效果最優(yōu)的數(shù)據(jù)集。目前，最常用的統(tǒng)計指標相位標準差能快速識別最佳對流層延遲改正數(shù)據(jù)集，但未能從空間尺度的角度去衡量。半變異函數(shù)和顧及距離變化的相位標準差不僅可表達對流層延遲改正量的空間關系，還可從空間尺度上評估對流層延遲改正模型的效果[9]。本研究利用覆蓋日本關東平原中部2016年12月—2020年8月的57景Sentinel 1A影像獲取時序干涉原始相位，采用2種統(tǒng)計指標(顧及距離變化的相位標準差和半變異函數(shù))定性與定量評估2種GAM(ERA5與MERRA2)和GACOS對流層延遲改正的效果。

1 研究區(qū)域與數(shù)據(jù)源

1.1 研究區(qū)域概況

日本關東平原西、北部接山地和丘陵，中部包含臺地和洼地，東部和南部臨太平洋和東京灣，面積約為1.6萬 km2，大部分海拔位于100 m以下。該地區(qū)屬于亞熱帶季風氣候，受海洋作用影響，具有明顯的海洋性氣候特征，氣溫年較差較小；沿岸的日本暖流導致年降水偏多，降雨量達1 500 mm以上；夏、秋季多臺風，對流層中水汽含量豐富[10]。世界上類似關東平原地形、緯度的地區(qū)如中國的華北平原。如圖1所示，研究區(qū)位于35.03°N～36.67°N，139.09°E～140.40°E(紅色矩形所示)，面積約13 300 km2。該區(qū)域內(nèi)中部地形相對開闊，林地分布較少，人口密度大，擁有世界上經(jīng)濟高度發(fā)達的東京都市圈。Morishita使用2015—2020年Sentinel 1A影像監(jiān)測東京周圍地表隆升達5 mm/a，在Tochigi地區(qū)出現(xiàn)季節(jié)性的沉降[11]。為了利用InSAR技術更加精細地反映關東平原地表形變的時空演化規(guī)律，需要最大限度減弱或去除對流層延遲對形變結果的干擾。

圖1 研究區(qū)示意圖

1.2 數(shù)據(jù)簡介

“哨兵一號”(Sentinel 1A/B)衛(wèi)星向公眾免費提供高質量的SAR影像，寬幅模式的地面覆蓋范圍約250 km×160 km，分辨率為5 m×20 m(距離向×方位向)[12]。本研究獲取2016年12月8日—2020年8月25日覆蓋研究區(qū)的Sentinel 1A影像，具體參數(shù)信息見表1。歐空局提供與影像日期對應的精密定軌星歷數(shù)據(jù)；美國國家航空航天局(national aeronautics and space administration，NASA)提供30 m分辨率數(shù)字高程模型(shuttle radar topography mission)[13]。

表1 Sentinel 1A影像參數(shù)信息

本研究采用的3種對流層延遲改正數(shù)據(jù)集包括：①ERA5是ECMWF于2018年發(fā)布的時間分辨率最高的模型，采用了含有衛(wèi)星數(shù)據(jù)改正、濕度分析等方法的四維同化技術(4D-VAR)處理后提供的全球氣象再分析資料，時間分辨率1 h，空間分辨率0.25°×0.25°，壓力層37層[14]；②GACOS是采用迭代對流層分解模型，融合全球連續(xù)性GPS的對流層延遲估計數(shù)據(jù)、地形高程數(shù)據(jù)與高分辨率的ECMWF產(chǎn)品，提供1979年至今的全球對流層延遲產(chǎn)品，時間分辨率6 h，空間分辨率0.125°×0.125°，壓力層137層，是目前最高的模型[7]；③MERRA2是NASA利用大氣數(shù)據(jù)同化系統(tǒng)對衛(wèi)星、地面觀測站、無線電探空及艦船的觀測數(shù)據(jù)進行質量控制、增量分析及同化，提供1980年至今的全球氣象再分析資料，時間分辨率6 h，空間分辨率0.5°×0.625°，壓力層42層[15]。

2 研究方法

2.1 GAM對流層延遲改正原理

GAM在全球尺度上以均勻分布的網(wǎng)格點形式提供相同時間間隔(1 h或6 h)的氣象參數(shù)，利用Sentinel 1A衛(wèi)星在研究區(qū)的成像時間提取GAM提供的溫度、氣壓和位勢高等氣象要素計算對流層延遲。利用大氣折射率N來表示對流層延遲，由于云層中水滴引起的液態(tài)水延遲對C波段的影響小于0.1～0.4 mm/km，可忽略不計。N由靜力延遲和濕延遲組成，可用式(1)表示[16]，利用式(2)對N積分計算得到對流層延遲[7]：

(1)

(2)

式中：Nhydro和Nwet指靜力延遲、濕延遲；P指大氣壓，hPa；T指溫度，K；e指水汽分壓，hPa；k1、k2、k3系數(shù)值分別是0.776、0.716和3.75×103K2Pa-1；θ指雷達信號入射角，°；λ指雷達波長，m；z指地面高程，m；zref指參考高度值15 km，認為超過15 km以上的高度，大氣折射率不隨高程改變。

根據(jù)影像獲取時間t，提取GAM中網(wǎng)格點上的氣象參數(shù)，利用相對濕度計算水汽分壓，再用式(3)～(5)計算GAM網(wǎng)格點上的天頂絕對對流層延遲[16]：

(3)

(4)

(5)

式中：P=Pd+e，為總氣壓；Rd和Rv分別為靜力延遲氣體和濕延遲氣體常數(shù)，287.05、461.495 Jkg-1K-1；gm為重力加速度。

由于GAM的分辨率較低，在垂直向上采用三次樣條插值及水平向上采用雙線性插值得到每景SAR影像中像素點的天頂對流層延遲。由于SAR衛(wèi)星系統(tǒng)側視成像，InSAR技術只能監(jiān)測視線向的形變信息，因此需要利用每個像素的入射角信息將天頂方向的對流層延遲投影到視線向，以主、輔影像中對應像素的對流層延遲相減得到干涉圖的對流層延遲，并將對流層延遲改正數(shù)從原始干涉圖中扣除來完成對流層延遲改正。Sentinel 1A衛(wèi)星為C波段衛(wèi)星，波長為5.6 cm，利用相位和距離之間的轉換關系，采用式(6)將對流層延遲轉換成對流層延遲相位：

(6)

式中，φ為對流層延遲相位，d為對流層延遲距離。

2.2 統(tǒng)計指標

采用不同的數(shù)據(jù)集改正對流層延遲可能會引入數(shù)據(jù)集本身的誤差，嚴重影響形變解譯的可靠性和精度，合適的評估指標對于識別最佳改正數(shù)據(jù)集尤為重要。統(tǒng)計指標相位標準差基于全局或局部范圍內(nèi)的干涉原始相位中無較大形變信號的準則來評估對流層延遲改正前后的效果[1]。因此，篩選未進行對流層延遲改正而解算的地表形變速率±0.5 mm/a之內(nèi)的數(shù)據(jù)點為評估數(shù)據(jù)集。本研究采用以下2種評價指標衡量對流層延遲改正前后的效果：

1) 考慮距離變化的相位標準差(standard deviation,STD)。隨機選擇100個半徑從1～100 km的窗口，計算每個相同大小的窗口內(nèi)未進行對流層延遲改正和每種數(shù)據(jù)集進行對流層延遲改正后所有數(shù)據(jù)點的相位標準差，如式(7)所示。隨著窗口距離增大會包含更多的數(shù)據(jù)點，然后計算100個大小相同窗口的相位標準差的平均值，作為識別時序InSAR對流層延遲改正最佳數(shù)據(jù)集的指標：

(7)

式中，干涉圖數(shù)量M，干涉原始相位、3種數(shù)據(jù)集改正對流層延遲后的相位均表示為xi。

2) 半變異函數(shù)。該指標表示設定固定的間隔滯后距離h，計算變量在x和x+h處的樣本方差值。距離相近的數(shù)據(jù)點的性質較為相似，半方差值隨空間距離的增大而增大，如式(8)所示[9]。為提高計算效率，對前期篩選的數(shù)據(jù)點進行降采樣(1∶100)，設置滯后距離h為1 km，從0～100 km，計算數(shù)據(jù)點對的半方差值：

(8)

式中，干涉原始相位、3種數(shù)據(jù)集改正對流層延遲后相位在點x處的值為Z(x)，h為滯后距離，Var[·]為計算采樣點方差函數(shù)。

2.3 數(shù)據(jù)處理流程

本研究時序InSAR干涉圖的對流層延遲估計與改正的具體技術路線見圖2。①基于SNAP軟件進行時序差分干涉：選擇時間2018年8月12日的影像為主影像，其余56景輔影像分別與主影像進行配準，將配準后的影像進行差分干涉處理得到時序差分干涉圖；②識別永久散射體點(persistent scatters，PS)：采用StaMPS軟件根據(jù)振幅離差指數(shù)(0.4)和相干點的相位空間選擇長時間內(nèi)保持穩(wěn)定雷達散射特性的PS點[17]，得到時序PS點集；③估計與改正對流層延遲：以TRAIN平臺計算GACOS、MERRA2和ERA5估計的對流層延遲[18]；④評估最優(yōu)的對流層延遲改正數(shù)據(jù)集：采用顧及距離變化的相位標準差和半變異函數(shù)評估對流層延遲改正前后的效果，得到最佳改正數(shù)據(jù)集。

圖2 評估3種數(shù)據(jù)集改正對流層延遲技術路線圖

3 結果與討論

3.1 對流層延遲改正結果

對預處理得到的56幅時序InSAR干涉圖進行對流層延遲改正，采用3種數(shù)據(jù)集改正后給出不同的改正表現(xiàn)，選擇2期單幅干涉圖進行說明。單幅2018年7月19日—2018年8月12日干涉圖時間基線24 d，空間基線-49.4 m，圖3(a)給出了3種數(shù)據(jù)集估計的對流層延遲差分相位與干涉圖改正前后的結果對比，對流層延遲差分相位呈明顯的空間變化特征，從相模灣、東京灣沿岸向內(nèi)陸地區(qū)(從東南向西北)遞增。GACOS模型估計的對流層延遲相位在空間上表現(xiàn)出更平滑且有利于捕捉小尺度變化，但GACOS改正后的結果顯示對流層延遲相位仍殘留在西北山區(qū)地帶。圖4(a)、4(c)顯示在25～100 km的空間尺度上，3種數(shù)據(jù)集均有明顯的對流層延遲改善效果。如表2所示，計算了整幅干涉圖GACOS、ERA5及MERRA2改正后的相位標準差分別減少了33.5%、27.9%及25.8%。

單幅2018年8月12日—2019年4月9日干涉圖時間基線216 d，空間基線-15.7 m，采用3種數(shù)據(jù)集估計的對流層延遲差分相位具有明顯的差異，其中圖3(b)顯示MERRA2在相模灣沿岸表現(xiàn)異常。圖4(b)、4(d)顯示MERRA2在13～100 km的空間尺度上，相比另外2種數(shù)據(jù)集改正結果出現(xiàn)了較大的偏離，可能與MERRA2的空間分辨率有關。在50～100 km的空間尺度上，GACOS與ERA5表現(xiàn)出較為明顯的改善效果。如表2所示，整幅干涉圖GACOS、ERA5改正后的相位標準差分別減少了41.7%和12.9%，而MERRA2改正之后相位標準差增加了65.9%。

表2 干涉圖對流層延遲改正后結果對比

圖3 3種數(shù)據(jù)集對InSAR干涉圖對流層延遲的估計與改正結果

圖4 統(tǒng)計指標評估2幅干涉圖對流層延遲改正的效果

采用上述2種統(tǒng)計指標評估了56幅時序InSAR干涉圖對流層延遲改正前后的效果，并以評估結果的平均值作為總體評價指標。由圖5可知2種統(tǒng)計指標的結果基本一致，不同的顏色帶表示57景影像生成56幅時序干涉圖，在相同的空間尺度上，每幅干涉圖經(jīng)過對流層延遲改正后會對應56個評估值，然后計算56個評估值的標準差，表明每種數(shù)據(jù)集改正時序干涉圖對流層延遲后評估值的離散程度。相位標準差是子區(qū)域內(nèi)對流層改正后相位值變化的總體指標，顧及距離變化的相位標準差是在100個子區(qū)域內(nèi)各自計算相位標準差后再取所有子區(qū)域標準差的平均值。半變異函數(shù)只判斷計算滯后距離內(nèi)目標點的相位值，該指標能夠在離散的空間尺度上有效捕捉改正后相位值的空間變化。在大于30 km的空間尺度上，半變異函數(shù)值明顯表現(xiàn)出ERA5比GACOS能夠更有效改正干涉圖中對流層延遲的優(yōu)勢。56幅干涉圖對流層延遲改正前后的評估結果平均值顯示0～30 km的空間尺度上，3種數(shù)據(jù)集均不能有效去除對流層延遲的影響，其中在大于30 km的空間尺度上MERRA2比GACOS和ERA5出現(xiàn)了較大的偏離。在30～100 km的空間尺度上，ERA5能夠有效改正對流層延遲，但GACOS與原始干涉相位的曲線相近也取得了微弱的改正效果。

圖5 統(tǒng)計指標評估56幅InSAR干涉圖對流層延遲的平均結果

另外，從定量角度比較了3種數(shù)據(jù)集改正對流層延遲的效果，表3給出了統(tǒng)計整幅干涉圖相位標準差的結果，可看出GACOS、ERA5及MERRA2對流層延遲改正后干涉圖的相位標準差數(shù)量分別降低了38、36和12幅，對應的干涉圖平均降低了21.6%、27.3%及17.6%。ERA5改正后的相位標準差平均值為1.18 cm，相比原始相位標準差平均值(1.43 cm)降低了17.5%，GACOS改正后相比原始值降低了13.3%。因此，相比GACOS及MERRA2，采用ERA5進行估計與改正時序干涉圖對流層延遲取得了較好的效果。

表 3 3種數(shù)據(jù)集時序InSAR對流層延遲改正結果

3.2 地表形變解算結果

根據(jù)上述3種數(shù)據(jù)集進行對流層延遲改正后的評估結果可知ERA5改正表現(xiàn)較好，利用ERA5去除對流層延遲后，采用PSInSAR方法獲取了研究區(qū)雷達視線向形變速率，如圖6(a)所示。研究區(qū)內(nèi)的地表形變機制可分為自然因素和人類活動2個方面：復雜的地質環(huán)境主要影響隆升區(qū)，發(fā)生在東京周圍((圖6(b))。 由于2011年Mw9.0東北大地震可能導致該地區(qū)的地表滲透率增強和孔隙壓力改變， 最大的隆起區(qū)東南部邊界清晰，可能受地下地質邊界的影響，這與前人研究結果一致[8]。人類活動的影響主要體現(xiàn)在沉降區(qū)域，在古賀-野木町周圍(圖6(c))和相模灣沿岸(圖6(d))發(fā)生了地表下沉，古賀-野木町周圍擁有易收縮的土壤，可能由于夏季農(nóng)業(yè)地下水使用引起；而相模灣沿岸沉降主要由獨特的地質構造背景和圍填海工程引起[11]。

(a)形變速率場;(b)隆升區(qū)域分布在東京周圍;(c)沉降1區(qū)域分布在古賀-野木周圍;(d)沉降2分布在相模灣沿岸

利用衛(wèi)星平臺的方位角和每個PS點的入射角信息，將GPS的NEU方向的觀測值投影到SAR衛(wèi)星系統(tǒng)的雷達視線向并計算GPS站點的形變速率以驗證InSAR監(jiān)測結果的可靠性。以每個GPS站點為中心，半徑100 m內(nèi)的InSAR觀測值進行反距離加權插值得到該站點的InSAR觀測值并計算形變速率。以G171站點為參考點計算其他GPS站點的相對形變速率，發(fā)現(xiàn)InSAR與GPS的較差小于2.5 mm/yr，具有較高的一致性(見表4)。

表4 InSAR與GPS結果比較

3.3 討論

日本關東平原地區(qū)的時序干涉圖中對流層延遲明顯，受海洋的影響導致該地區(qū)大氣對流活動頻繁，更容易產(chǎn)生空間小尺度且無規(guī)則的湍流混合延遲。研究區(qū)內(nèi)大部分地表高程變化較小，僅在地形起伏劇烈的西北山區(qū)地帶少量存在由大氣厚度變化引起的垂直分層延遲。對流層延遲改正前后的相位標準差均值和削減量反映了不同改正數(shù)據(jù)集對InSAR對流層延遲的敏感性具有差異，可能是由于不同改正數(shù)據(jù)集與Sentinel 1A影像的采集時間存在不同的時間差。GAM依賴于地面觀測氣象資料的采集，關東平原地區(qū)毗鄰太平洋，易受臺風等事件影響引起水汽的快速變化，難以被現(xiàn)有精度的GAM所捕捉。不同時空分辨率、氣壓分層數(shù)的數(shù)據(jù)集提供的氣象參數(shù)也存在差異性，也導致了改正效果不同，圖7表明了ERA5和MERRA2提供的溫度、位勢和相對水汽含量在第二壓力層的顯著差異性。MERRA2的時空分辨率比GACOS、ERA5更低，提供的大氣參數(shù)難以計算研究區(qū)內(nèi)對流層的真實延遲。GACOS融合了全球連續(xù)GPS天頂延遲數(shù)據(jù)進行增強處理，提高了空間與時間上的改正精度，研究區(qū)內(nèi)GACOS略低于ERA5的改正效果。ERA5在時間上更接近Sentinel 1A影像的獲取時間，進而在時序InSAR對流層延遲改正中降低了對流層估計的不確定性。

圖7 ERA5與MERRA2在第2壓力層的溫度、水汽及位勢參數(shù)比較(白色區(qū)域為空值)

4 結論

本研究以覆蓋日本關東平原的57景Sentinel 1A影像為例，采用考慮距離變化的相位標準差和半變異函數(shù)評估分析了GACOS、ERA5及MERRA2這3種數(shù)據(jù)集改正時序InSAR干涉圖對流層延遲的效果。結果表明，考慮距離變化的相位標準差和半變異函數(shù)可用于不同空間尺度對流層延遲改正效果的評估，二者評估結果一致。對于單幅干涉圖來說，3種數(shù)據(jù)集的改正效果具有差異；對于56幅時序InSAR干涉圖改正效果綜合評估來看，在0～30 km空間尺度上，3種數(shù)據(jù)集在該研究區(qū)內(nèi)的對流層延遲改正效果均不明顯；在30～100 km空間尺度上，ERA5與GACOS具有一定改正效果，MERRA2效果較差。

以日本關東平原作為典型的沿海平原地區(qū)，采用InSAR技術在海岸帶地區(qū)開展大范圍的地表形變監(jiān)測受對流層延遲誤差的嚴重影響，結果顯示，時間分辨率最高的ERA5數(shù)據(jù)集在該地區(qū)的對流層延遲改正中具有較好的效果，能夠為InSAR技術在國內(nèi)外沿海平原地區(qū)的應用提供參考。