李志進(jìn) 章傳銀 王 偉 徐鵬飛 梁圣豪

1 山東科技大學(xué)測(cè)繪與空間信息學(xué)院,青島市前灣港路579號(hào),266590 2 中國(guó)測(cè)繪科學(xué)研究院,北京市蓮花池西路28號(hào),100036

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

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

1 研究區(qū)和實(shí)驗(yàn)數(shù)據(jù)

1.1 研究區(qū)概況

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

1.2 實(shí)驗(yàn)數(shù)據(jù)

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

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

2.1 時(shí)序 InSAR處理

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

(1)

假設(shè)在tA、tB(tA

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

(2)

Bv=δφ

(3)

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

DU=DLOS/cosθ

(4)

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

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

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

2.2 InSAR監(jiān)測(cè)量?jī)?yōu)化方法

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

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

(5)

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

2.3 融合CORS網(wǎng)與時(shí)序InSAR協(xié)同監(jiān)測(cè)

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

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

(6)

通過(guò)融合CORS網(wǎng)與時(shí)序InSAR數(shù)據(jù),可以將時(shí)序InSAR高空間分辨率與CORS站數(shù)據(jù)高精度的優(yōu)點(diǎn)結(jié)合起來(lái),構(gòu)建具有高分辨率、高精度的地面垂直方向形變量時(shí)間序列。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 InSAR結(jié)果分析

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

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

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

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

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

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

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

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

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

3.4 融合CORS網(wǎng)與時(shí)序InSAR監(jiān)測(cè)結(jié)果分析

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

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

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

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

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

4 結(jié) 語(yǔ)

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

1)通過(guò)去除InSAR監(jiān)測(cè)時(shí)間序列中非動(dòng)力學(xué)信號(hào)形變,可以得到更加可靠且穩(wěn)定的地面形變數(shù)據(jù)。

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

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