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

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

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

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

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

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

1.1 研究區(qū)概況

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

1.2 實驗數(shù)據

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

2 數(shù)據預處理與融合方法

2.1 時序 InSAR處理

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

(1)

假設在tA、tB(tA

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

(2)

Bv=δφ

(3)

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

DU=DLOS/cosθ

(4)

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

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

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

2.2 InSAR監(jiān)測量優(yōu)化方法

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

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

(5)

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

2.3 融合CORS網與時序InSAR協(xié)同監(jiān)測

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

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

(6)

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

3 實驗結果與分析

3.1 InSAR結果分析

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

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

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

3.2 InSAR優(yōu)化方法結果

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

表1 InSAR監(jiān)測量處理前后數(shù)據統(tǒng)計Tab.1 Statistics of InSAR monitoring data before and after processing

3.3 CORS站數(shù)據處理結果

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

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

3.4 融合CORS網與時序InSAR監(jiān)測結果分析

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

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

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

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

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

4 結 語

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

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

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

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