白希選 張兵強(qiáng) 徐 昊 湛舒文 嚴(yán) 翼

1 武漢工程大學(xué)土木工程與建筑學(xué)院，武漢市雄楚大道693號(hào)，430074 2 中國(guó)科學(xué)院精密測(cè)量科學(xué)與技術(shù)創(chuàng)新研究院，武漢市徐東大街340號(hào)，430077 3 中國(guó)科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院，北京市玉泉路19號(hào)甲，100049 4 中南民族大學(xué)資源與環(huán)境學(xué)院，武漢市民族大道182號(hào)，430074

合成孔徑雷達(dá)差分干涉測(cè)量(D-InSAR)技術(shù)是目前地表形變監(jiān)測(cè)的重要手段，但雷達(dá)信號(hào)在大氣中傳播時(shí)會(huì)受到空氣中各類介質(zhì)的影響，導(dǎo)致相位延遲或提前。對(duì)流層和電離層是大氣中造成信號(hào)延遲的主要介質(zhì)層，其中對(duì)流層延遲是D-InSAR 形變測(cè)量最主要的誤差源[1-3]。D-InSAR對(duì)流層校正方法可分為兩大類，一類是通過(guò)GPS、MERIS、MODIS等外部數(shù)據(jù)建立大氣誤差模型[4-6]，如GACOS(generic atmospheric correction online service for InSAR)系統(tǒng)[7-9]；另一類是利用InSAR對(duì)流層延遲在時(shí)空上的統(tǒng)計(jì)特性，基于統(tǒng)計(jì)學(xué)方法來(lái)削弱其影響[10]，即擬合研究區(qū)LOS向形變相位與數(shù)字高程模型(digital elevation model，DEM)[11-13]進(jìn)行誤差校正。

目前已有學(xué)者對(duì)GACOS與DEM數(shù)據(jù)在時(shí)序InSAR形變監(jiān)測(cè)中大氣校正效果進(jìn)行研究[10]，但是未考慮大氣噪聲與地形(相位-高程)函數(shù)中形變信號(hào)對(duì)擬合函數(shù)的影響，且不適用于同震形變觀測(cè)。施賀青等[13]基于DEM擬合相位-高程函數(shù)并降低大氣誤差，但是未考慮選擇不同擬合區(qū)域?qū)M合函數(shù)與校正精度的影響。本文從地震活動(dòng)頻發(fā)地區(qū)選取震級(jí)MS5.4～7.0事件。首先以精河地震為例，通過(guò)研究區(qū)數(shù)據(jù)分析得到基于外部DEM進(jìn)行大氣誤差校正時(shí)擬合區(qū)域的選取方法，然后以外部DEM數(shù)據(jù)分別擬合與地形相關(guān)的大氣誤差模型和GACOS數(shù)據(jù)建立的大氣誤差模型進(jìn)行誤差校正；最后依據(jù)實(shí)驗(yàn)結(jié)果分析兩種大氣校正方法的適用性和校正精度。

1 數(shù)據(jù)源與誤差分析

1.1 數(shù)據(jù)源

Sentinel-1A/B為C波段SAR衛(wèi)星，在高緯度地區(qū)的重訪周期最短可以達(dá)到1 d[14]。本文選取10對(duì)覆蓋地震影響區(qū)域的Sentinel-1A/B衛(wèi)星干涉寬幅模式(IW)雷達(dá)影像(表1)，基于GMTSAR軟件，采用兩軌法處理SAR數(shù)據(jù)。

表1 SAR數(shù)據(jù)基本信息

1.2 大氣誤差分析

大氣延遲誤差主要來(lái)源于電離層和對(duì)流層， C波段SAR衛(wèi)星受電離層干擾較小，且兩次成像期間電離層延遲效應(yīng)變化較小，可在數(shù)據(jù)處理過(guò)程中抵消大部分延遲信號(hào)。對(duì)流層約包含75%的大氣質(zhì)量和幾乎全部的大氣水汽，天頂對(duì)流層總延遲(ZTD)可以表示為[15]:

(1)

式中，Rd為干空氣氣體常數(shù)，ρ為干空氣總質(zhì)量密度，Zw為可壓縮系數(shù)，Pw為水汽局部氣壓，T為大氣溫度，k1、k2、k3為大氣折射常數(shù)，h0與h分別為測(cè)站與對(duì)流層頂部對(duì)應(yīng)的高程。式(1)右邊第1項(xiàng)為天頂靜力學(xué)延遲(ZHD)，第2項(xiàng)為天頂濕延遲(ZWD)，SAR影像中大氣延遲主要由對(duì)流層中濕延遲造成，引起的相位誤差可表示為：

Δφ=φ1-φ2=

(2)

2 基于DEM數(shù)據(jù)擬合大氣誤差模型

現(xiàn)有研究結(jié)果表明，大氣延遲與地形具有相關(guān)性[2,10,16]?；贗nSAR技術(shù)提取形變場(chǎng)時(shí)，通過(guò)擬合高程與形變場(chǎng)大氣噪聲可得到噪聲隨高程變化的函數(shù)模型，然后扣除與地形變化相關(guān)的大氣延遲信號(hào)。雖然同震形變信號(hào)遠(yuǎn)小于高程，但一般強(qiáng)于大氣噪聲，擬合大氣誤差函數(shù)時(shí)若未扣除形變區(qū)域，擬合函數(shù)將受到真實(shí)形變信號(hào)的顯著影響。Albino等[16]基于DEM對(duì)形變場(chǎng)進(jìn)行大氣誤差校正，但在擬合相位-高程函數(shù)時(shí)未分離形變信號(hào)與大氣噪聲，導(dǎo)致函數(shù)模型中包含真實(shí)形變信號(hào)的影響。

2.1 不同擬合區(qū)域參數(shù)分析

為避免同震形變對(duì)擬合相位-高程函數(shù)產(chǎn)生影響，在擬合大氣延遲信號(hào)時(shí)扣除干涉圖中有明顯同震形變條紋的區(qū)域，剩余相位解纏后的區(qū)域用于擬合兩次SAR成像期間大氣延遲信號(hào)。相位解纏共選擇13個(gè)區(qū)域，其形變邊界逐步擴(kuò)大，如圖1(a)所示，其中黑色矩形表示選擇的所有實(shí)驗(yàn)區(qū)相位解纏窗口，紅色矩形為擬合函數(shù)時(shí)扣除的區(qū)域，面積為1 020 km2。圖1(b)為擬合區(qū)面積與相位-高程相關(guān)系數(shù)，其中橫坐標(biāo)表示擬合區(qū)面積。從圖中可以看出，精河地震區(qū)域形變信號(hào)受大氣噪聲影響較大，包含大氣噪聲的形變信號(hào)與高程具有強(qiáng)相關(guān)性，相關(guān)系數(shù)絕對(duì)值最低為0.839 1，最高為0.915 7，擬合區(qū)域面積為2 760 km2時(shí)達(dá)到最大值，之后緩慢下降。為避免形變信號(hào)對(duì)標(biāo)準(zhǔn)差的影響，圖1(c)、1(d)、1(e)為去除同震影響區(qū)域的大氣校正前后標(biāo)準(zhǔn)差及其差值，其中橫坐標(biāo)為扣除同震影響區(qū)域后的面積。從圖中可以看出，大氣校正前后標(biāo)準(zhǔn)差隨擬合區(qū)面積增加表現(xiàn)出明顯的變化趨勢(shì)，擬合面積為2 350～2 760 km2時(shí)，校正后標(biāo)準(zhǔn)差降到最小值且校正前后標(biāo)準(zhǔn)差變化量趨近最大值，即此時(shí)校正效果最優(yōu)。同時(shí)，從圖1(g) 、1(h)可以看出，研究區(qū)域總面積小于3 800 km2時(shí)，大氣校正后均方根下降速度較快，隨后降速變緩。

2.2 擬合區(qū)域范圍

通過(guò)§2.1分析可知，基于DEM數(shù)據(jù)擬合大氣誤差函數(shù)進(jìn)行同震形變場(chǎng)大氣校正時(shí)，其結(jié)果主要受擬合區(qū)域影響。相位-高程函數(shù)可通過(guò)同震影響區(qū)域以外的數(shù)據(jù)擬合得到，擬合大氣誤差時(shí)在較小范圍內(nèi)擴(kuò)大或縮小扣除的同震影響區(qū)域(圖1(a)紅色矩形區(qū)域)對(duì)誤差函數(shù)模型與校正結(jié)果基本無(wú)影響。因此，擬合相位-高程函數(shù)時(shí)扣除的同震影響區(qū)域只需觀察干涉圖條紋進(jìn)行選擇。

根據(jù)均方根、標(biāo)準(zhǔn)差及大氣校正前后同震形變對(duì)比等綜合分析可知，擬合區(qū)域面積選擇在2 300～2 800 km2時(shí)，大氣校正前后LOS向形變值與標(biāo)準(zhǔn)差、均方根變化均較為穩(wěn)定，相關(guān)系數(shù)也接近幅值，可取得較好的大氣校正效果。大氣分布具有時(shí)空差異，擬合區(qū)面積大且與形變區(qū)域較遠(yuǎn)時(shí)，大氣延遲隨高程的變化關(guān)系與形變場(chǎng)中心大氣分布差異性越大；擬合區(qū)面積小則相關(guān)系數(shù)低，無(wú)法準(zhǔn)確擬合高程與形變之間的變化關(guān)系，校正后噪聲同樣較大。因此，根據(jù)校正前后均方根、標(biāo)準(zhǔn)差變化值及LOS向形變場(chǎng)綜合分析可知，精河地震同震形變場(chǎng)擬合區(qū)域面積為2 300 ～2 800km2時(shí)噪聲相對(duì)較小且變化穩(wěn)定。圖2為擬合區(qū)域面積為2 350 km2和2 760 km2時(shí)大氣校正前后形變場(chǎng)，從圖中可以看出，LOS向最大形變量穩(wěn)定在7.0～7.5 cm之間。

圖1 研究區(qū)與相關(guān)參數(shù)Fig.1 Study area and related parameters

紅色矩形為擬合相位-高程相關(guān)模型時(shí)去除的區(qū)域，紅色矩形之外為擬合區(qū)域圖2 不同擬合區(qū)校正結(jié)果Fig.2 Correction results of different fitting areas

3 誤差校正分析

Albino等[16]基于GACOS與DEM數(shù)據(jù)對(duì)2017年印尼火山噴發(fā)前后形變進(jìn)行大氣校正，Bekaert等[10]通過(guò)冪律函數(shù)對(duì)大氣噪聲進(jìn)行校正，但上述方法均未分析研究區(qū)形變?cè)跀M合誤差函數(shù)時(shí)的影響?？紤]到上述因素以及不同區(qū)域的地形差異，在對(duì)其余研究案例進(jìn)行大氣校正時(shí)，以精河地震誤差函數(shù)模型的擬合區(qū)為參考，根據(jù)同震影響范圍選擇擬合區(qū)域。從數(shù)據(jù)處理結(jié)果可知，利用GACOS數(shù)據(jù)校正后的形變場(chǎng)可能會(huì)存在顯著趨勢(shì)，本文分別在是否設(shè)置參考點(diǎn)的情況下利用GACOS進(jìn)行誤差校正。此外，為考慮軌道誤差造成的影響[17]，通過(guò)上述兩種方法進(jìn)行大氣誤差校正后擬合并去除形變場(chǎng)趨勢(shì)。

3.1 基于DEM的大氣誤差校正

圖3為校正前后的形變場(chǎng)、誤差模型與相位-高程線性函數(shù)。表2為每個(gè)地震事件對(duì)應(yīng)的相位-高程相關(guān)系數(shù)、標(biāo)準(zhǔn)差、均方根等參數(shù)。圖3(a)為拜城地震，對(duì)比校正前后的形變場(chǎng)可以看出，校正后幅值下降約3 cm，均方根下降2.6 cm。圖3(d)中大氣校正之前很難分辨其形變信號(hào)，然而LOS向形變與DEM相關(guān)系數(shù)達(dá)0.8，校正后噪聲大幅減小，形變量約為5 cm。圖3(e)中校正后均方根下降3.3 cm。圖3(g)中校正前斷層兩盤沿LOS向形變幅值為6 cm和-9 cm，校正后形變量為8 cm和-7 cm，噪聲下降約2 cm，比張克亮等[18]的GNSS研究結(jié)果以及楊九元等[19]的InSAR研究結(jié)果更具可靠性。圖3(i)為近海地區(qū)，與地形無(wú)關(guān)的對(duì)流層延遲信號(hào)一般較為顯著，基于DEM數(shù)據(jù)進(jìn)行校正后噪聲下降相對(duì)較小。圖3中其他事件的大氣噪聲與高程相關(guān)系數(shù)較低，基于DEM進(jìn)行校正后整體均方根與標(biāo)準(zhǔn)差變化微弱。

3.2 基于GACOS數(shù)據(jù)的大氣誤差校正

基于GACOS數(shù)據(jù)的校正結(jié)果如圖3中a5～j5所示，均方根及標(biāo)準(zhǔn)差等參數(shù)見(jiàn)表2。為提高數(shù)據(jù)處理精度，將GACOS大氣延遲信號(hào)逐像素投影至LOS向，然后分別在未設(shè)置和設(shè)置參考點(diǎn)的情況下進(jìn)行大氣校正。未設(shè)置參考點(diǎn)時(shí)部分研究案例校正后噪聲反而增加，均方根增加值最大可達(dá)9 cm左右，校正后形變值可靠性低。根據(jù)形變場(chǎng)均值設(shè)置參考點(diǎn)或擬合趨勢(shì)，校正后的形變值與現(xiàn)有較為可靠的研究結(jié)果相近[13,18-19]。

3.3 誤差校正結(jié)果分析與討論

綜合圖3中形變場(chǎng)以及圖4中校正前后參數(shù)變化趨勢(shì)可知，相關(guān)系數(shù)0.45以上的地震事件基于DEM函數(shù)模型可取得較好的校正效果，均方根降幅在1～4 cm之間；相關(guān)系數(shù)大于0.8時(shí)，基于DEM數(shù)據(jù)的校正效果整體優(yōu)于GACOS。由表2和圖4可知，在相位-高程相關(guān)系數(shù)低于0.4時(shí)，DEM無(wú)法擬合可靠函數(shù)模型，而利用GACOS數(shù)據(jù)進(jìn)行誤差校正，均方根與標(biāo)準(zhǔn)差均下降，從形變場(chǎng)也可觀察到噪聲顯著下降。相關(guān)系數(shù)為0.4～0.8時(shí)，使用GACOS校正也可達(dá)到與DEM校正相近的結(jié)果，可以使用兩種方法相互驗(yàn)證。因此，基于兩類數(shù)據(jù)進(jìn)行大氣誤差校正時(shí)可以互補(bǔ)，噪聲與高程相關(guān)性越強(qiáng)，說(shuō)明地形越復(fù)雜且存在顯著大氣噪聲，此時(shí)GACOS數(shù)據(jù)受GNSS站密度與ECMWF(European centre for medium-range weather forecasts)等數(shù)據(jù)源精度影響，校正精度相對(duì)較低；而相關(guān)性越強(qiáng)，說(shuō)明噪聲分布更依賴于高程，只使用DEM擬合大氣誤差函數(shù)可獲得更為精確的結(jié)果。

4 結(jié) 語(yǔ)

本文以大氣噪聲顯著的精河地震為例，建立基于DEM進(jìn)行大氣誤差校正時(shí)擬合區(qū)的選取方法；然后分別基于DEM和GACOS建立的誤差模型對(duì)不同地區(qū)與震級(jí)事件的LOS向形變數(shù)據(jù)進(jìn)行校正，并分析校正精度和不同方法的適用性，得到以下結(jié)果：

表2 研究區(qū)大氣校正前后相關(guān)參數(shù)

圖中a～j依次為表1中序號(hào)1～10對(duì)應(yīng)的地震事件；圖中紅色矩形為擬合相位-高程函數(shù)時(shí)扣除的主要同震影響區(qū)域，矩形以外為擬合區(qū)圖3 基于DEM和GACOS數(shù)據(jù)大氣延遲校正Fig.3 Atmospheric delay correction based on DEM and GACOS data

(c)中有無(wú)參考點(diǎn)對(duì)應(yīng)的標(biāo)準(zhǔn)差一致，因此只體現(xiàn)出3組數(shù)據(jù)圖4 大氣校正前后參數(shù)隨相關(guān)性變化趨勢(shì)Fig.4 Variation trend of parameters with correlation before and after atmospheric correction

1)考慮到不同發(fā)震區(qū)域地形差異以及擬合區(qū)與形變區(qū)大氣延遲的空間差異性，基于DEM擬合大氣誤差模型時(shí)，最佳擬合區(qū)并非固定值，在一定范圍內(nèi)均可取得較好的大氣校正結(jié)果。以精河地震為例，相位-高程相關(guān)系數(shù)絕對(duì)值臨近最大值且擬合區(qū)域面積在2 300～2 800 km2范圍內(nèi)時(shí)，校正效果整體較好。因此，本文未給出通過(guò)多次實(shí)驗(yàn)確定擬合面積的計(jì)算方程，而是給出相關(guān)性與校正精度隨擬合區(qū)面積的變化趨勢(shì)，以供擬合誤差函數(shù)時(shí)參考。

2)相位-高程相關(guān)系數(shù)絕對(duì)值大于0.4時(shí)，必須考慮大氣誤差對(duì)形變值的影響；相關(guān)系數(shù)小于0.4時(shí)，噪聲受高程影響較小，可基于GACOS進(jìn)行誤差校正；相關(guān)系數(shù)為0.4～0.8時(shí)，可使用兩種校正方法相互驗(yàn)證；相關(guān)系數(shù)大于0.8時(shí)，基于DEM校正后的均方根與標(biāo)準(zhǔn)差均低于GACOS，強(qiáng)相關(guān)性也表明形變場(chǎng)噪聲分布顯著依賴于高程，可單獨(dú)基于DEM擬合噪聲分布模型。

3)在本文所有研究案例中，基于上述兩種方法校正的數(shù)據(jù)結(jié)果在消除趨勢(shì)后均方根與標(biāo)準(zhǔn)差整體下降，因此，大氣校正后可通過(guò)對(duì)形變場(chǎng)進(jìn)行多項(xiàng)式擬合以消除部分趨勢(shì)。