基于實測氣象參數(shù)的InSAR大氣校正及其應用研究

黃嘉艇 孔 揚 王 盼 占 果 唐 旭

1 南京信息工程大學遙感與測繪工程學院,南京市寧六路219號,210044 2 寧波市氣象局,浙江省寧波市氣象路118號,315012 3 華設設計集團股份有限公司,南京市紫云大道9號,210014

合成孔徑雷達干涉測量(InSAR)技術因其監(jiān)測范圍大、時空分辨率高和監(jiān)測周期長等優(yōu)勢[1],在地面高程獲取和地表形變監(jiān)測領域得到廣泛應用,但由于受對流層延遲效應的嚴重影響,其在高精度形變監(jiān)測領域的應用有限[2-3]。時序InSAR技術采用時空濾波算法,能在很大程度上削弱對流層延遲誤差,但僅能消除大氣湍流混合延遲效應,無法削弱大氣垂直分層延遲效應[4]。因此,提高時序InSAR技術的大氣延遲校正效果成為國內(nèi)外研究的重點[5-6]。

目前我國利用地面氣象數(shù)據(jù)校正InSAR對流層延遲的研究較少。基于此,本文對浙江省寧波市鄞州區(qū)2018～2020年58景哨兵影像作干涉處理,得到相位解纏后的干涉圖,利用傳感器獲取的地面氣象數(shù)據(jù)和美國國家環(huán)境預報中心(NCEP)提供的氣象再分析資料,結合大氣校正模型進行對流層延遲誤差校正,并利用相位標準差STD評估2種數(shù)據(jù)的校正效果,以分析不同數(shù)據(jù)對大氣延遲校正模型的適用性。

1 研究區(qū)概況和實驗數(shù)據(jù)

鄞州區(qū)東起北侖港、寧波保稅區(qū),西部與余姚接壤,南瀕奉化地區(qū),北臨長江三角洲,屬于寧波市重要經(jīng)濟樞紐。研究區(qū)地理范圍為121.13°～121.90°E、29.61°～29.95°N,占地面積約為1 381 km2。

本文選取歐空局提供的2018-01～2020-12共58景極化方式為VV的Sentinel-1A影像,條帶為寬幅模式,覆蓋范圍約為250 km×160 km,空間分辨率為20 m×5 m。為更好地進行干涉測量,外部DEM數(shù)據(jù)選取美國國家航空航天局提供的空間分辨率為30 m的數(shù)字高程模型,去除地形相位的影響;選取歐空局提供的精密軌道數(shù)據(jù)對軌道信息進行修正,去除平地相位的影響。

選用2種數(shù)據(jù)集削弱對流層延遲:1)傳感器每分鐘記錄的氣壓、溫度等地面氣象數(shù)據(jù)。測站均勻分布在監(jiān)測區(qū)域,相鄰兩測站間的距離最大為38.47 km、最小為5.31 km;2)NCEP提供的中國區(qū)域大氣再分析資料,包含氣壓、溫度、水汽混合比等氣象參數(shù),水平分辨率為0.5°×0.5°,垂直方向的氣壓層共分為37層,時間分辨率為6 h,數(shù)據(jù)集時間為UTC 00:00、06:00、12:00和18:00。

2 數(shù)據(jù)處理

2.1 SAR數(shù)據(jù)干涉處理

本文采用SBAS-InSAR技術監(jiān)測鄞州區(qū)地表形變,數(shù)據(jù)處理流程如圖1所示。以覆蓋研究區(qū)的58景SAR影像為數(shù)據(jù)源,對影像進行裁剪。綜合考慮多普勒質心頻譜差等因素后,選取2019-05-29最優(yōu)SAR影像為主影像,其余57景SAR數(shù)據(jù)為副影像,將副影像分別與主影像進行配準。設置時間基線閾值為50 d,空間基線閾值為臨界基線的45%,經(jīng)SAR影像干涉處理后共生成149對干涉對。為更好地進行干涉測量,分別引入精密軌道數(shù)據(jù)和外部DEM數(shù)據(jù)以去除平地相位和地形相位的影響。選用Goldstein自適應濾波去除干涉處理時存在的噪聲相位。為更好地處理相干性高的區(qū)域,采用Delaunay MCF方法進行相位解纏。

圖1 數(shù)據(jù)處理流程

由于鄞州區(qū)臨近海域,對流層空氣濕度較大,與濕度相關的大氣延遲較為嚴重,傳統(tǒng)的時空濾波算法只能去除大氣湍流混合延遲,對于垂直方向的大氣延遲并不適用。因此,必須引入外部氣象數(shù)據(jù)并建立大氣校正模型來求解天頂方向的對流層總延遲ZTD。由于微波信號的傳播路徑為雷達視線方向,故需將雷達視線向的大氣延遲投影到天頂方向,對解纏后的干涉圖引入改正量天頂總延遲ZTD,以削弱大氣延遲誤差的影響,提高解纏圖的質量,從而更好地反演地面年平均形變速率和地表累積形變量。

2.2 大氣延遲相位去除

雷達信號在傳播過程中會受大氣折射的影響,產(chǎn)生路徑延遲和路徑彎曲,使得InSAR監(jiān)測過程中形成大氣延遲誤差。利用每幅干涉圖對應的2期ZTD改正量進行計算,得到對應的大氣相位,將解纏后的干涉圖與之差分相減,以達到減弱大氣延遲誤差的目的。天頂總延遲ZTD可表示為:

ZTD=ZHD+ZWD

(1)

式中,ZHD為天頂靜力學延遲,ZWD為天頂濕延遲。

天頂靜力學延遲ZHD可通過Saastamoinen模型求解得到[7]。Saastamoinen模型是目前常用的對流層延遲改正模型之一,預測精度可達mm級,模型公式如下:

ZHD=

(2)

式中,es為氣象站點處的氣壓,φ為站點位置的緯度,h為站點大地高。

(3)

(4)

式中,di為氣象站到最近4個格網(wǎng)點的距離,Pi為格網(wǎng)點處的壓強,Ti為格網(wǎng)點處的溫度。NCEP再分析數(shù)據(jù)集時間分辨率為6 h,每個氣象站點每日只能記錄4個壓強數(shù)據(jù),為獲取衛(wèi)星過境時站點附近的氣象數(shù)據(jù),采用三次樣條插值法(cubic spline interpolation,CSI)獲取特定時刻的壓強。同理,站點附近的溫度也可通過上述方法求得。

NCEP記錄的溫度為地球表面數(shù)據(jù),而記錄的氣壓為平均海平面數(shù)據(jù),式(2)中氣象站周圍的實際氣壓值es可通過下式求得:

es=PresMSL×

(5)

天頂濕延遲ZWD與大氣加權平均溫度Tm、大氣可降水量PWV密切相關,通過一個轉換系數(shù)可求得天頂方向的大氣水汽量,即

(6)

式中,Rv為水汽氣體常數(shù),取值461J·(kg·K)-1;k1、k2、k3均為大氣折射常數(shù),分別取值77.6 K·(hPa)-1、71.98 K·(hPa)-1、3.754×105K2·(hPa)-1;mv和md分別為水蒸氣和干空氣的摩爾質量,取值18.015 2 g·(mol)-1、28.964 4 g·(mol)-1;Tm為大氣加權平均溫度。

大氣加權平均溫度Tm是求解ZWD的關鍵參數(shù)。目前使用最廣泛的是Bevis提出的實用模型Tm=70.2+0.72Ts,即Tm與地表溫度Ts呈線性關系[8]。但Tm模型受氣候與地形等因素的限制,多適用于歐美區(qū)域,在中國區(qū)域的應用精度較低。本文利用分布在中國地區(qū)的84個探空站資料和對應的氣象數(shù)據(jù),建立研究區(qū)的Tm回歸經(jīng)驗模型。同時,為更好地表示加權平均溫度,在模型中引入水汽壓[9],計算公式如下:

Tm=α0+α1×T+α2×e

(7)

式中,T為氣象站點處的絕對溫度,e為水汽壓。利用84個探空站數(shù)據(jù)估算多項式系數(shù)α0、α1、α2,分別取值92.61、0.634、0.279 7。

水汽壓e需通過飽和水汽壓公式并結合露點溫度間接求得,即

(8)

式中,Td為露點溫度。

結合美國國家環(huán)境預報中心提供的氣象再分析資料可知,PWV可通過對地面上空各氣壓層中的可降水量進行數(shù)值積分求得,計算公式為[10]:

(9)

式中,g為研究區(qū)的重力加速度,取值9.793 6 m/s2;Ps、Pt分別為對流層頂氣壓和地面氣壓;r為水汽混合比;m為氣象再分析資料記錄的大氣層數(shù);Pi和ri分別為第i層大氣壓力和水汽混合比。其中,地面到對流層頂部的氣壓和水汽混合比可從NCEP數(shù)據(jù)中獲取。

3 實驗結果與分析

利用2種數(shù)據(jù)集對處理后的149幅干涉圖進行對流層延遲校正,篇幅有限,本文以2019-04-30～2019-03-25干涉對為例進行分析,校正效果如圖2所示。由圖2(a)可見,原始干涉圖的相位標準差為1.68,校正前的大氣延遲噪聲較多;圖2(b)和2(c)的相位標準差分別為1.37和1.33,大氣延遲相位均有所降低,說明2種數(shù)據(jù)集均能有效改善對流層延遲。但利用NCEP氣象再分析資料削弱大氣噪聲的效果差于實測氣象數(shù)據(jù),對實測氣象數(shù)據(jù)進行大氣校正后噪聲相位基本消失,效果最佳。

圖2 2019-04-30～2019-03-25干涉圖對流層延遲校正前后對比

由圖3可知,經(jīng)大氣校正后的大部分干涉相位標準差都有所降低,其中采用NCEP氣象再分析資料進行大氣校正后相位標準偏差降低的干涉圖共計100幅,占67.1%;采用地面氣象信息進行大氣校正后相位標準偏差降低的干涉圖共計106幅,占71.1%。由此說明,地面氣象信息具有更好的校正效果。

圖3 2種數(shù)據(jù)大氣校正后的相位標準差

3.1 精度評估

地面氣象站能夠對靠近地面的大氣層氣象要素進行實時觀測,具有較高的時空分辨率。為更好地利用氣象再分析資料計算天頂總延遲的精度,本文以實測氣象數(shù)據(jù)為真值,以再分析資料提供的數(shù)據(jù)為觀測值,計算各站點處再分析資料求解ZTD的平均偏差bias和均方根誤差RMSE,結果如表1(單位mm)所示。偏差bias和均方根誤差RMSE的計算公式為:

(10)

表1 再分析資料計算ZTD的精度

(11)

從表1可以看出,部分站點求解ZTD的bias和RMSE相對較大,最大偏差和均方根誤差分別為15.98 mm、16.05 mm。以最大誤差站點58569為例,對該氣象站點插值處的氣象數(shù)據(jù)進行分析。如圖4(a)和4(b)所示,2種氣象數(shù)據(jù)的一致性相對較好,證明利用反距離加權插值算法和三次樣條插值算法獲取的氣象數(shù)據(jù)具有可靠性。由圖4(c)可知,部分數(shù)據(jù)存在較大偏差。由于NCEP數(shù)據(jù)自身存在精度誤差,因此在解算對流層延遲時會出現(xiàn)誤差疊加,導致校正效果降低。

圖4 58569站地面氣象數(shù)據(jù)插值結果與實測數(shù)據(jù)比較

為進一步驗證InSAR監(jiān)測結果的準確性,將雷達視線向的地表形變投影至垂直方向,并引入2018～2020年外部水準數(shù)據(jù)進行對比驗證。選取分布在研究區(qū)的7個水準點數(shù)據(jù),分別與未加大氣校正、利用NCEP氣象再分析資料校正和利用地面氣象信息校正后的監(jiān)測數(shù)據(jù)進行對比,結果如表2(單位mm)所示。從表可以看出,未進行大氣校正的監(jiān)測結果精度較低,與水準測量的監(jiān)測結果差值較大,RMSE為5.62 mm;利用NCEP氣象再分析資料校正后的監(jiān)測結果精度有所提高,RMSE為3.86 mm;而經(jīng)地面氣象信息校正后的形變值與水準測量所得的監(jiān)測數(shù)據(jù)最為接近,RMSE為2.78 mm。綜上所述,經(jīng)過地面氣象信息和NCEP氣象再分析資料大氣校正后的監(jiān)測精度都有所改善,但利用地面氣象信息進行大氣校正的效果最好、可靠性最高。

表2 SBAS-InSAR監(jiān)測結果與水準數(shù)據(jù)對比

利用SBAS-InSAR技術對實測氣象數(shù)據(jù)大氣校正后的地區(qū)進行監(jiān)測分析,獲取鄞州區(qū)2018～2020年垂直方向的地表形變結果,如圖5所示(負值代表沉降,正值代表抬升)。由圖可見,大部分地區(qū)的形變都處于一個相對穩(wěn)定的狀態(tài),即未發(fā)生較大形變。但仍存在幾個明顯的沉降漏斗區(qū),本文選取5個沉降較為嚴重的區(qū)域展開分析。

圖5 鄞州區(qū)大氣校正后的地面沉降圖

3.2 沉降區(qū)域分析

A區(qū)域位于云龍鎮(zhèn)南部鄞城大道附近,該道路于2019-12正式通車,年平均形變速率約為-6.7～-55.8 mm/a,與甬新河交界處的累積沉降量最大達到165 mm。引起沉降的原因主要為監(jiān)測期間該段道路處于修建狀態(tài),大量工程材料的堆積和澆筑,施工期間對地面壓實不夠及修建完成后巨大車流量對道路的長期荷載等,而甬新河附近土質結構疏松,多為沙性土質且不易結成土塊,所以在河流附近的沉降量最大。

B區(qū)域位于首南街道西北部江山萬里小區(qū)附近,年平均形變速率約為-3.7～-48.2 mm/a,最大累積沉降量達142.3 mm。該區(qū)域形成沉降漏斗的主要原因是在形變監(jiān)測過程中,江山萬里小區(qū)六期和七期高層建筑物均處于施工階段,地基的開挖、地下水的大量抽取和重型載貨汽車的過往等均會對地面土質產(chǎn)生擾動,且施工區(qū)域臨近奉化江,土質相對較軟,導致地面的沉降現(xiàn)象更加顯著。

C區(qū)域位于東錢湖鎮(zhèn)西南部錢湖景苑小區(qū)附近,年平均形變速率約為-18.8～-54.1 mm/a,最大累積沉降量達147.1 mm。實地踏勘發(fā)現(xiàn),該小區(qū)南側寧波軌道交通正在建設施工,判斷該地區(qū)沉降主要原因為所在地基受軌道施工影響。

D區(qū)域位于邱隘鎮(zhèn)北部寧波藝術實驗學院教育集團明湖校區(qū)附近,年平均形變速率約為-2.6～-60.8 mm/a,最大累積沉降量達179.8 mm。該區(qū)域屬于新開發(fā)區(qū)域,整個監(jiān)測期間處于施工階段,大量工程建設活動對地表有較大影響,引發(fā)地面出現(xiàn)沉降漏斗。

E區(qū)域位于五香鎮(zhèn)南部四都裝村,附近有地鐵1號線,年平均形變速率約為-5.8～-56.7 mm/a,最大累積沉降量達164.7 mm。沉降原因主要有2點:1)該區(qū)域每天都有列車經(jīng)過,列車循環(huán)碾壓軌道的動荷載效應累積引發(fā)沉降;2)監(jiān)測期間正值寶涵路修建,壓實度不達標,加上車輛長期荷載對地面帶來的影響引起地面沉降。該區(qū)域四面環(huán)河,土質松軟,也加大了沉降風險。

4 結 語

1)利用NCEP氣象再分析資料校正后相位標準差降低的干涉對占67.1%,而經(jīng)過地面氣象信息校正后相位標準差降低的干涉對占71.1%,說明利用地面氣象信息能更好地削弱對流層延遲誤差的影響。

2)以地面氣象信息計算得到的ZTD為真值,NCEP氣象再分析資料計算得到的ZTD為觀測值。結果表明,最大偏差為15.98 mm,平均均方根誤差為4.29 mm。

3)未進行大氣校正的監(jiān)測結果最差,RMSE為5.62 mm;利用NCEP氣象再分析資料進行大氣校正后的監(jiān)測結果略有改善,RMSE為3.86 mm;而利用地面氣象信息進行大氣校正后的監(jiān)測結果能更真實地反映地面沉降情況,RMSE為2.78 mm。

4)2018～2020年鄞州區(qū)出現(xiàn)不規(guī)則地面沉降,選取5個沉降區(qū)域進行分析,結果表明,鄞州區(qū)地面沉降主要受地下水過度抽取、鐵路隧道開挖、高層建筑物施工及重型載貨汽車過往等因素影響。