融合CORS網與時序InSAR研究北京地區地面形變

李志進 章傳銀 王 偉 徐鵬飛 梁圣豪

1 山東科技大學測繪與空間信息學院,青島市前灣港路579號,266590 2 中國測繪科學研究院,北京市蓮花池西路28號,100036

全球衛星導航系統(GNSS)是一種傳統的大地測量技術手段,具有全天候、高精度、實時服務能力強的特點,可直接獲取地面監測點的形變信息,但空間分辨率低,難以對整個區域進行監測。近年來新興的合成孔徑雷達干涉測量(interferometric synthetic aperture radar,InSAR)技術具有低成本、大范圍、高空間分辨率的特點,可以得到cm級甚至mm級地表形變監測結果。InSAR技術主要有合成孔徑雷達差分干涉測量(DInSAR)技術、永久散射體干涉測量(PS-InSAR)技術和短基線集干涉測量(SBAS-InSAR)技術。其中,SBAS-InSAR采用多主影像,時空基線為短基線,可有效解決時空失相干的問題,適用性更強。GNSS與InSAR兩種監測方法各有優缺點,將兩種方法相結合可以得到數據質量更優且更可靠的地面監測數據。楊國創等[1]對InSAR和GNSS技術所得的形變量進行分析,結果表明監測的沿海地區地表形變幅度和趨勢高度一致,說明兩者具備數據融合的條件。周文韜[2]提出利用方差分量估計方法融合兩者數據的三維形變監測并取得良好效果。但該方法缺少InSAR監測量的優化處理,由于InSAR技術監測的地表形變量包含自然或人為破壞、地表氣溫變化、降雨等氣象因素引起的誤差,無法直接用于與CORS數據融合處理。

本文對兩種技術所得的數據進行優化處理,首先通過時序InSAR處理獲取研究區地表形變量,將其分離出數米以下的巖土層形變量;然后以衛星導航定位基準站(CORS)數據為基準,對InSAR監測結果進行整體平差,修復其差分干涉尺度誤差,補償空間中長波誤差影響,控制InSAR監測量隨時間的累積誤差等,從而實現高分辨率、高精度的地面形變監測。

1 研究區和實驗數據

1.1 研究區概況

北京地處我國華北平原北部,自20世紀70年代以來,由于城市的開發建設,對地下水需求加大,導致地下水過度開采,同時城市高層建筑不斷增多,造成朝陽、通州等地出現嚴重的地面沉降[3]。對地面沉降進行大范圍、高效的時空監測不僅能為城市建設、經濟發展提供可靠的科學依據,也有助于做好地質災害防范工作。

1.2 實驗數據

時序InSAR處理使用31景Sentinel-1A衛星影像數據,起止時間為2018-01-03～2020-11-12。使用POD精密定軌星歷數據修正軌道誤差,采用空間分辨率為30 m的SRTM1 DEM數據共同參與本次處理。本文所用的CORS站數據為昌平站(CHPN)、牛欄山站(NLSH)、東三旗站(DSQI)、朝陽站(CHAO)、西集站(XIJI) 2018～2020年連續觀測數據。

2 數據預處理與融合方法

2.1 時序 InSAR處理

短基線集干涉測量(SBAS-InSAR)原理為:假設研究區內影像有N+1幅,選擇一幅影像為超級主影像,并且任意一幅影像至少可以與其他影像形成一個干涉對,共生成M幅差分干涉圖,M滿足以下條件:

(1)

假設在tA、tB(tA

δφj(x,r)=φ(tB,x,r)-φ(tA,x,r)≈

(2)

Bv=δφ

(3)

式中,B為系數矩陣,v為變化率,δ表示估計值,φ為相對差分相位。B在滿秩時可利用最小二乘法求解,秩虧時需用奇異值分解法求得最小范數解,進而得到研究區形變速率[4]。由式(4)可將LOS向形變量時間序列轉換為大地高方向形變量時間序列:

DU=DLOS/cosθ

(4)

式中,DU為大地高方向形變量,DLOS為雷達視線向形變量,θ為衛星雷達方向入射角。

本次實驗采用SARscape軟件,選擇SBAS-InSAR處理方法,時間基線閾值為120 d,共生成87對干涉對,時空基線分布如圖1所示。

圖1 時空基線分布Fig.1 Distribution of spatio-temporal baselines

2.2 InSAR監測量優化方法

由于InSAR監測量中存在野值、粗差和突變等非動力學形變信號,需要將其進行探測與分離。以InSAR監測量采樣歷元時刻為單元,由高斯低通濾波器構造低通監測量參考面,再以3倍標準差為閾值進行粗差探測。高斯低通濾波器[5]定義如下:

H(u,v)=e-D(u,v)2/2σ2

(5)

由于在每個采樣歷元中所分離的粗差點不同,為保證監測量在整個時序上的時空分布,依據動力學地面垂直形變量與動力源作用點距離或距離平方近似成反比的空間變化性質,以InSAR監測量采樣歷元時刻為單元,按高斯濾波算法計算粗差點的監測量,對粗差點進行修復,從而保證監測點的時空分布不變。通過對InSAR監測量進行優化處理,以抑制或削弱非地質動力學作用的極淺地表局部變化影響,從而得到地表數米以下的深巖土層垂直形變[6]。

2.3 融合CORS網與時序InSAR協同監測

采用GAMIT/GLOBK軟件進行CORS站連續觀測數據的基線解算,獲得單日解文件,計算CORS站大地高周變化時間序列。對時間序列進行粗差探測,并分離線性項和常數項,重構非線性項。采用連續Fourier和切比雪夫組合基函數,由不規則采樣時序構造低通濾波曲線參考基準,計算采樣值與參考值之差,并對殘差時序進行統計,將大于指定倍數殘差標準差的采樣記錄作為粗差并分離。剔除粗差后,利用連續切比雪夫基函數分離CORS站大地高周變化時間序列的常數項和線性項[7],并對非線性項進行低頻重構。將重構后的非線性項與線性項、常數項相加,得到新的CORS站垂直方向周變化時間序列。

去除InSAR監測量中極淺地表局部變化影響,得到地表數米以下的深巖土層垂直形變量,與CORS站數據進行融合。選取一個采樣歷元作為CORS數據與InSAR數據融合的參考歷元,將兩部分數據的參考歷元進行統一。由CORS站周邊InSAR監測量內插得到CORS站處的InSAR監測量,對CORS站大地高變化時序按照時間內插得到InSAR采樣歷元時刻的CORS站大地高變化。以CORS站為基準,根據每期InSAR散點間監測量的相對變化量,采用間接平差方法對全部InSAR監測量進行整體平差,從而實現CORS網與InSAR監測時序的高度統一與高精度傳遞。平差模型[8]為:

(6)

通過融合CORS網與時序InSAR數據,可以將時序InSAR高空間分辨率與CORS站數據高精度的優點結合起來,構建具有高分辨率、高精度的地面垂直方向形變量時間序列。

3 實驗結果與分析

3.1 InSAR結果分析

將北京地區SBAS-InSAR處理結果轉化為垂直方向形變量時間序列,繪制年平均形變速率,結果如圖2所示。由圖可知,沉降比較明顯的地區為朝陽區與通州區交界處、大興區南部、海淀區北部,沉降速率超過-60 mm/a,最大沉降速率約為-110 mm/a;抬升比較明顯的地區為昌平區中西部、門頭溝區和海淀區交界處以及房山區中北部,這些地區地表年平均形變速率約為10～20 mm/a。該結果與文獻[9]中結論大致相同。

圖2 InSAR監測量年平均形變速率Fig.2 Average annual deformation rate of InSAR monitoring data

將昌平站(CHPN)、牛欄山站(NLSH)、東三旗站(DSQI)、朝陽站(CHAO)、西集站(XIJI)5個CORS站在大地高方向的年平均變化率與SBAS-InSAR結果在大地高方向的年平均變化率進行比較,發現位于昌平區中心位置的CHPN站與位于順義區北部的NLSH站略微抬升,而位于昌平區東南部的DSQI站、朝陽區中心位置的CHAO站和通州區東部的XIJI站發生沉降,且CHAO站沉降最為明顯。CORS站數據與SBAS-InSAR結果相吻合,表明數據具有可靠性。

3.2 InSAR優化方法結果

對SBAS-InSAR結果進行優化處理后,選取采樣歷元進行對比。表1(單位mm)為2018-07-02監測量進行優化處理前后的數據統計,由表可知,處理后的標準差更小,平均值幾乎不變,并剔除了極值,表明InSAR數據去除極淺地表局部變化影響后獲得的數據質量更可靠。圖3為該歷元下監測量優化前后結果對比,由圖可知,處理后的誤差點得到明顯修復,進一步說明了優化方法的有效性。

表1 InSAR監測量處理前后數據統計Tab.1 Statistics of InSAR monitoring data before and after processing

3.3 CORS站數據處理結果

表2(單位mm/a)為2018～2020年5個CORS站在大地高方向的年平均變化率,其中CHPN和NLSH站表現為略微抬升,DSQI、CHAO和XIJI站呈現沉降趨勢,并且DSQI站沉降最為顯著,3 a間累積沉降量為87 mm。圖4為北京地區5個CORS站大地高周變化原始數據與其線性變化和非線性項重構結果。

表2 CORS站大地高年平均變化率Tab.2 Average annual change rate of geodetic height at CORS stations

3.4 融合CORS網與時序InSAR監測結果分析

以CORS數據為基準,對InSAR監測量進行間接平差,計算SAR影像覆蓋范圍內的散點在2018～2020年間年平均形變速率,結果如圖5所示。

圖5 數據融合后年平均形變速率Fig.5 Average annual deformation rate after data fusion

由圖5可知,朝陽區與通州區交界處地面沉降最嚴重,最大沉降速率達到-90 mm/a;昌平區中西部、海淀區西部、門頭溝區東部、石景山區、順義區北部地面抬升比較明顯,年平均形變速率約為10～20 mm/a;其他地區沉降變化不明顯,沉降速率約為-10～10 mm/a。

對比SBAS-InSAR結果可知,地面沉降與抬升情況一致。選取2018-03-04、2018-11-11和2020-09-13三個采樣歷元,利用距離反比方法獲取參考歷元為2019-04-21 00:00的SBAS-InSAR監測量以及數據融合后的監測量在昌平站等5個CORS站處的形變量,與CORS站大地高周變化時序相比并分別作差,結果見表3(單位mm)。由表可知,融合CORS網與時序InSAR的監測量更接近CORS站的監測量。

表3 CORS站形變量對比Tab.3 Comparison of monitoring data at CORS stations

4 結 語

本文基于SBAS-InSAR處理獲取的地表形變時間序列和CORS站大地高周變化時間序列,經過優化處理后,以CORS網為基準對InSAR數據進行間接平差,得到融合后的地面形變時間序列。對北京地區2018～2020年地面沉降進行分析,得到以下結論:

1)通過去除InSAR監測時間序列中非動力學信號形變,可以得到更加可靠且穩定的地面形變數據。

2)融合CORS網與時序InSAR進行監測整體可行,可將InSAR技術的高空間分辨率特點與CORS網的高精度特點結合起來,實現高分辨率、高精度的地面形變監測。融合CORS網與時序InSAR在監測方面具有很大優勢。

3)融合后的地面沉降監測結果顯示,朝陽區、通州區、大興區中南部地面沉降最為明顯,最大沉降速率達到-90 mm/a;昌平區中西部、海淀區西部、門頭溝區東部、石景山區、順義區北部存在較為明顯的抬升,年平均形變速率約為10～20 mm/a;其他地區地面形變相對較小,年平均形變速率約為-10～10 mm/a。