多軌道時(shí)序InSAR擬穩(wěn)基準(zhǔn)平差方法研究

杜釗鋒 周曉敏 程小凱 鄭文科 朱新杰

1 自然資源部第一大地測(cè)量隊(duì),西安市測(cè)繪路4號(hào),710054

在利用時(shí)序InSAR(interferometric synthetic aperture radar)技術(shù)進(jìn)行大范圍地面沉降監(jiān)測(cè)時(shí),通常涉及多個(gè)軌道的數(shù)據(jù)融合,其中最主要的問題是不同軌道形變場(chǎng)參考基準(zhǔn)的統(tǒng)一[1]。許多學(xué)者[2-3]針對(duì)多軌道時(shí)序InSAR技術(shù)展開研究,Shirzaei[4]和Pepe等[5]對(duì)多軌道多時(shí)相InSAR拼接算法進(jìn)行分析,提出對(duì)公共區(qū)域目標(biāo)點(diǎn)在空間范圍和時(shí)間序列的擬合進(jìn)行多軌道數(shù)據(jù)基準(zhǔn)統(tǒng)一的方法。部分學(xué)者[1,6]利用Envisat、Radarsat、Sentinel等衛(wèi)星數(shù)據(jù)對(duì)我國(guó)主要沉降區(qū)開展大范圍InSAR監(jiān)測(cè),通過(guò)計(jì)算多軌沉降結(jié)果的重疊區(qū)偏移值,將同名點(diǎn)的形變速率取算術(shù)平均,以實(shí)現(xiàn)大區(qū)域形變速率的鑲嵌。文獻(xiàn)[7-8]對(duì)大范圍多軌道基準(zhǔn)轉(zhuǎn)換與數(shù)據(jù)拼接問題作進(jìn)一步研究,分析引起相鄰軌道沉降速率差異的因素,提出采用區(qū)塊法、插值法等方法解決異軌永久散射體(persistent scatterer, PS)沉降量偏差。

然而,上述研究大都基于重疊區(qū)評(píng)估恒定的偏差量,通過(guò)消除偏差量來(lái)實(shí)現(xiàn)基準(zhǔn)統(tǒng)一,但該類方法容易將不同軌道參考基準(zhǔn)的差異當(dāng)作誤差組分,尤其對(duì)形變量較大的區(qū)域其拼接后的形變速率容易產(chǎn)生扭曲。因此,本文通過(guò)對(duì)多軌道時(shí)序InSAR誤差源進(jìn)行分析與處理,引入擬穩(wěn)基準(zhǔn)平差方法,消除多軌道形變監(jiān)測(cè)結(jié)果的偏差,達(dá)到不同軌道形變場(chǎng)參考基準(zhǔn)的統(tǒng)一,增強(qiáng)多軌道時(shí)序InSAR技術(shù)監(jiān)測(cè)大范圍地面沉降的適用性。

1 多軌道時(shí)序InSAR誤差源分析

1.1 入射角引起的差異

入射角引起的差異主要由SAR斜視特征所決定,InSAR數(shù)據(jù)處理得到的直接結(jié)果為視線向變化量。如圖1所示,S1和S2表示相鄰軌道2次衛(wèi)星過(guò)境,假設(shè)某點(diǎn)真實(shí)地面沉降量為L(zhǎng),SAR數(shù)據(jù)近地點(diǎn)和遠(yuǎn)地點(diǎn)入射角分別為α、α+β,雷達(dá)視線向變化量分別為L(zhǎng)cosα、Lcos(α+β),差值為ΔL=Lcosα-Lcos(α+β)。以哨兵數(shù)據(jù)為例,當(dāng)真實(shí)地面沉降速率為50 mm/a時(shí),雷達(dá)視線向近地點(diǎn)和遠(yuǎn)地點(diǎn)最大沉降偏差為12.3 mm/a;當(dāng)真實(shí)地面沉降速率為100 mm/a時(shí),雷達(dá)視線向近地點(diǎn)和遠(yuǎn)地點(diǎn)最大沉降偏差為24.5 mm/a。由此可以推斷,地面沉降越嚴(yán)重的區(qū)域由入射角引起的差異越明顯。

圖1 雷達(dá)入射角差異示意圖Fig.1 Diagram of radar incident angle difference

本文研究主要考慮的是垂直方向地面沉降量,水平方向變化量由地殼水平運(yùn)動(dòng)速度場(chǎng)決定,屬于板塊運(yùn)動(dòng)尺度,在本文研究區(qū)內(nèi)可認(rèn)為是均質(zhì)的。因此,數(shù)據(jù)處理中將視線向形變量轉(zhuǎn)換到垂直方向,不考慮水平方向位移。

1.2 參考點(diǎn)引起的差異

參考點(diǎn)引起的差異與相干目標(biāo)點(diǎn)距參考點(diǎn)的距離有關(guān),一般來(lái)說(shuō),參考點(diǎn)本身也存在一定程度的年際和季節(jié)性變化,選取參考點(diǎn)首先應(yīng)考慮穩(wěn)定性,即參考點(diǎn)在監(jiān)測(cè)周期內(nèi)處于相對(duì)穩(wěn)定的狀態(tài),但隨著監(jiān)測(cè)范圍的逐步擴(kuò)大,參考點(diǎn)距離增加,由單個(gè)參考點(diǎn)得到的形變量逐漸失真。朱武等[9]對(duì)多參考點(diǎn)PS-InSAR數(shù)據(jù)處理方法進(jìn)行研究,對(duì)比單參考點(diǎn)和多參考點(diǎn)下的形變場(chǎng),認(rèn)為多參考點(diǎn)結(jié)果比單參考點(diǎn)結(jié)果更接近地面觀測(cè)資料。前人研究范圍較小,而當(dāng)研究范圍不斷擴(kuò)大時(shí),參考點(diǎn)引起的差異將更加明顯。

本文數(shù)據(jù)處理選擇多參考點(diǎn),將綜合分析選取的擬穩(wěn)點(diǎn)作為多參考點(diǎn)進(jìn)行數(shù)據(jù)處理。為保證各條帶形變參考基準(zhǔn)的統(tǒng)一,相鄰條帶間至少保證重疊2個(gè)擬穩(wěn)點(diǎn)。

1.3 軌道殘余誤差

軌道殘余誤差主要由軌道估計(jì)的偏差所致,在數(shù)據(jù)處理結(jié)果上體現(xiàn)為趨勢(shì)性偏差,評(píng)估軌道殘余誤差的數(shù)學(xué)模型可表示為:

1階:z=ax+by+c

1.5階:z=ax+by+cxy+d

2階:z=ax2+by2+cxy+d

3階:z=ax3+by3+cx2y+dy2x+ex2+

fy2+gxy+h

本文采用精密軌道星歷文件計(jì)算基線,考慮到按上述模型進(jìn)行軌道誤差評(píng)估時(shí)可能導(dǎo)致形變分量被誤評(píng)估為軌道誤差,因此數(shù)據(jù)處理中暫不進(jìn)行軌道殘余誤差的評(píng)估。

1.4 大尺度地形相關(guān)性誤差

InSAR干涉測(cè)量中通常將大氣延遲分解為大氣垂直分層延遲和水平紊亂延遲[10],其中,大氣垂直分層延遲主要由靜力學(xué)延遲和小部分天頂濕延遲組成,該部分誤差與地形具有較強(qiáng)的相關(guān)性;水平紊亂延遲則主要由天頂濕延遲組成,該部分誤差在時(shí)間上具有隨機(jī)性特征。

本文主要研究大區(qū)域垂直形變場(chǎng)基準(zhǔn)的統(tǒng)一,大氣延遲以低頻信號(hào)為主,與地形具有較強(qiáng)的相關(guān)性,因此采用GACOS模型[11]對(duì)大氣延遲進(jìn)行改正,再計(jì)算平均沉降速率和形變時(shí)間序列。評(píng)估完上述誤差組分后,其余誤差主要由大氣誤差組成,以高頻信號(hào)為主,通過(guò)濾波方法即可去除,濾波時(shí)需考慮空間窗口和時(shí)間維度,保留一定程度的季節(jié)性變形幅度,最后得到相干目標(biāo)點(diǎn)的平均沉降速率和形變時(shí)間序列結(jié)果。

2 擬穩(wěn)基準(zhǔn)平差方法

擬穩(wěn)平差方法是針對(duì)變形監(jiān)測(cè)網(wǎng)中存在部分穩(wěn)定點(diǎn)和部分不穩(wěn)定點(diǎn)的情況,不同于監(jiān)測(cè)網(wǎng)的秩虧自由網(wǎng)平差計(jì)算,擬穩(wěn)平差中將所有點(diǎn)都當(dāng)作未知數(shù)進(jìn)行解算,顧及到變形監(jiān)測(cè)網(wǎng)中存在部分相對(duì)穩(wěn)定的點(diǎn)(稱為“擬穩(wěn)點(diǎn)”),對(duì)該類點(diǎn)附加相關(guān)的未知量范數(shù)極小的條件,此前的秩虧自由網(wǎng)便能求得唯一解。由擬穩(wěn)點(diǎn)構(gòu)成的參考基準(zhǔn)則稱為監(jiān)測(cè)網(wǎng)的擬穩(wěn)基準(zhǔn)。擬穩(wěn)平差主要應(yīng)用于沉降監(jiān)測(cè)分析,如水準(zhǔn)網(wǎng)沉降分析,本文將其應(yīng)用于InSAR監(jiān)測(cè)的地面沉降平差計(jì)算與分析。

設(shè)不穩(wěn)定點(diǎn)的未知數(shù)為X1,穩(wěn)定點(diǎn)的未知數(shù)為X2,則誤差方程為[12]:

(1)

式中,n表示InSAR監(jiān)測(cè)結(jié)果中PS點(diǎn)由不同參考點(diǎn)所組成的多余觀測(cè)量,t1和t2分別表示不穩(wěn)定PS點(diǎn)和穩(wěn)定PS點(diǎn)數(shù)量,t1+t2=t,R(A)=r,秩虧d=t-r,A2維數(shù)大于秩虧數(shù),即m>d,在InSAR監(jiān)測(cè)的地面沉降網(wǎng)中m>1,必須有2個(gè)或2個(gè)以上的擬穩(wěn)點(diǎn)。R(A1)=t-m,A1為列滿秩矩陣。觀測(cè)值l的權(quán)重P按距離定權(quán),類似于水準(zhǔn)監(jiān)測(cè)網(wǎng),P滿秩,按照最小二乘原則VTPV=min,構(gòu)成法方程組為NX=ATPl,可寫成:

(2)

(3)

(4)

將式(3)、(4)代入式(1)得到改正數(shù)的計(jì)算公式:

(5)

由以上公式可以看出,擬穩(wěn)平差是在最小二乘原則和最小范數(shù)條件下得到的,與自由網(wǎng)不同之處在于最小范數(shù)條件是僅相對(duì)擬穩(wěn)點(diǎn)而言,在使監(jiān)測(cè)網(wǎng)形不扭曲的同時(shí),得到監(jiān)測(cè)網(wǎng)擬穩(wěn)基準(zhǔn)的變化量。

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

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

根據(jù)研究區(qū)范圍獲取SAR影像,本文選擇Sentinel-1衛(wèi)星干涉寬幅成像模式升軌數(shù)據(jù),影像空間分辨率為5 m×20 m,覆蓋研究區(qū)共10景,每景超100期數(shù)據(jù),共計(jì)1 120期數(shù)據(jù),每期數(shù)據(jù)覆蓋范圍約為250 km×150 km,SAR數(shù)據(jù)如表1所示。

表1 SAR數(shù)據(jù)Tab.1 SAR data

3.2 處理結(jié)果

本文采用StaMPS方法[13]進(jìn)行時(shí)序InSAR數(shù)據(jù)處理,基于幅度離散和相位空間相關(guān)特性選點(diǎn),該方法能夠在城區(qū)和郊區(qū)識(shí)別出較多的PS點(diǎn)。數(shù)據(jù)處理中依據(jù)空間基線、時(shí)間基線、相干性等因素選取主圖像,40條帶、142條帶、69條帶主圖像分別為2018-11-28、2018-03-04、2018-12-24,干涉對(duì)空間基線在[-150, 150]之間,依次進(jìn)行圖像配準(zhǔn)、干涉處理、去除平地效應(yīng)、計(jì)算候選點(diǎn)、去除噪聲點(diǎn)、相位解纏、評(píng)估大氣影響相位、軌道誤差相位、地形殘差相位等[14],將形變相位轉(zhuǎn)換為年均沉降速率和時(shí)間序列變化量,即為初始數(shù)據(jù)處理結(jié)果。

本文研究區(qū)主要以平原為主,北部和西部有一定范圍的山區(qū),一般認(rèn)為山區(qū)較少發(fā)生大范圍的地表形變,選定為InSAR參考點(diǎn)。各條帶初始數(shù)據(jù)處理結(jié)果均是以InSAR參考點(diǎn)為基準(zhǔn),類似于水準(zhǔn)監(jiān)測(cè)網(wǎng)中給定固定點(diǎn)的間接平差計(jì)算結(jié)果。

在擬穩(wěn)點(diǎn)選取時(shí),充分考慮InSAR初始數(shù)據(jù)處理結(jié)果中變形量較小區(qū)域的空間分布,同時(shí)綜合已收集資料中CORS站地表形變結(jié)果、國(guó)家高程基準(zhǔn)網(wǎng)基巖點(diǎn)分布、大地控制點(diǎn)垂向變化結(jié)果等因素,首先選取一定數(shù)量的擬穩(wěn)點(diǎn)作為備選點(diǎn),通過(guò)一致性檢驗(yàn)剔除相對(duì)不穩(wěn)定的點(diǎn),最終確定參與平差的擬穩(wěn)點(diǎn)。如圖2所示,40條帶、142條帶、69條帶擬穩(wěn)點(diǎn)數(shù)量分別為5、5、4。通過(guò)編寫MATLAB代碼實(shí)現(xiàn)多軌道時(shí)序InSAR擬穩(wěn)基準(zhǔn)平差數(shù)據(jù)處理,該方法可以有效降低單參考點(diǎn)或區(qū)域造成的系統(tǒng)偏差,尤其是可以抑制由于參考距離較遠(yuǎn)造成誤差傳播遞增的現(xiàn)象,使測(cè)量平差計(jì)算結(jié)果更加符合客觀規(guī)律,提升地面沉降監(jiān)測(cè)網(wǎng)的可靠性。

經(jīng)過(guò)擬穩(wěn)基準(zhǔn)平差處理,得到華北平原40、142、69條帶InSAR數(shù)據(jù)處理結(jié)果,將年均沉降速率和年際沉降速率進(jìn)行疊加,結(jié)果如圖3所示。從圖中可以看出,各條帶之間沉降速率過(guò)渡自然,無(wú)明顯分界,說(shuō)明InSAR數(shù)據(jù)處理整體情況良好。

通過(guò)年際沉降速率(圖4)可以發(fā)現(xiàn),沉降漏斗演變和發(fā)展各有規(guī)律,華北平原主要沉降漏斗在2019年最為嚴(yán)重,漏斗逐年變化形態(tài)也有所差異,北京-天津地區(qū)沉降漏斗呈逐年減少趨勢(shì),河北地區(qū)(巨鹿縣-南宮市、臨漳縣-廣平縣-曲周縣、景縣-阜城縣、饒陽(yáng)縣-肅寧縣-高陽(yáng)縣)主要沉降漏斗呈持續(xù)高速發(fā)展趨勢(shì),沉降速率為-150～100 mm/a。沉降漏斗變化發(fā)展與年度降水量、地下水開采量、南水北調(diào)補(bǔ)水等存在一定關(guān)系,待后續(xù)進(jìn)一步研究。

圖4 年際沉降速率Fig.4 Interannual subsidence rate

依據(jù)InSAR監(jiān)測(cè)結(jié)果,采用克里格方法對(duì)InSAR相干目標(biāo)點(diǎn)平均沉降速率進(jìn)行插值計(jì)算,表2為平均沉降速率分類統(tǒng)計(jì)結(jié)果。

表2 InSAR結(jié)果平均沉降速率分類統(tǒng)計(jì)Tab.2 Statistics of average subsidence rate of InSAR results

3.3 精度評(píng)定

精度評(píng)定以內(nèi)符合精度為主,對(duì)40-142重疊區(qū)域、142-69重疊區(qū)域內(nèi)PS點(diǎn)進(jìn)行分析,分別提取PS點(diǎn)各年度平均沉降速率,通過(guò)建立一一對(duì)應(yīng)關(guān)系,將各年度平均沉降速率進(jìn)行差值,計(jì)算其平均誤差和中誤差,40-142重疊區(qū)域建立PS點(diǎn)樣本數(shù)目542 184個(gè),142-69重疊區(qū)域建立PS點(diǎn)樣本數(shù)目1 389 857個(gè)。

分別對(duì)初始數(shù)據(jù)處理結(jié)果和擬穩(wěn)基準(zhǔn)平差結(jié)果按照上述方法統(tǒng)計(jì)重疊區(qū)域PS點(diǎn)速率差值的平均誤差和中誤差,結(jié)果如表3(單位mm/a)、表4(單位mm/a)所示。從表中可以看出,通過(guò)擬穩(wěn)基準(zhǔn)平差后內(nèi)符合精度統(tǒng)計(jì)結(jié)果明顯優(yōu)于初始數(shù)據(jù)處理結(jié)果。由圖5(a)可知,擬穩(wěn)基準(zhǔn)平差后平均誤差優(yōu)于5 mm,中誤差優(yōu)于9 mm。

表3 初始數(shù)據(jù)處理結(jié)果內(nèi)符合精度統(tǒng)計(jì)Tab.3 Statistics of internal coincidence accuracy of initial data processing results

表4 擬穩(wěn)基準(zhǔn)平差后內(nèi)符合精度統(tǒng)計(jì)Tab.4 Statistics of internal coincidence accuracy after quasi-stable datum adjustment

圖5 InSAR結(jié)果與CORS站結(jié)果對(duì)比Fig.5 Comparison of InSAR and CORS results

本文通過(guò)研究區(qū)內(nèi)10個(gè)CORS站1999～2019年位移計(jì)算結(jié)果(中國(guó)地震局地殼形變臺(tái)網(wǎng)中心提供),與擬穩(wěn)基準(zhǔn)平差后的時(shí)序InSAR進(jìn)行對(duì)比分析。將CORS站與InSAR結(jié)果歸算到地表垂向變化時(shí)間序列上進(jìn)行數(shù)據(jù)對(duì)比,結(jié)果如圖5(b)～(h)所示。由圖可知,兩者吻合程度較好,BJFS、HELY等站點(diǎn)變化幅度與InSAR結(jié)果較為一致;BJSH、BJYQ、SXYC等站點(diǎn)變化幅度相比InSAR結(jié)果略小;HECX、TJBH、TJWQ等站點(diǎn)連續(xù)多年的地面沉降趨勢(shì)與InSAR結(jié)果整體一致,兩者存在一定程度的差異主要是由于空間位置無(wú)法完全重合所致,CORS站點(diǎn)為特定點(diǎn)位,而InSAR結(jié)果的相干目標(biāo)點(diǎn)代表一定區(qū)域的面散射體特征[15]?？傮w來(lái)講,CORS站位移計(jì)算結(jié)果與InSAR結(jié)果在刻畫季節(jié)性變化的細(xì)節(jié)形變信息上較為一致,從外符合精度角度進(jìn)一步說(shuō)明InSAR結(jié)果具有可靠性。

CORS站優(yōu)勢(shì)主要體現(xiàn)在單站點(diǎn)數(shù)據(jù)采樣足夠豐富,計(jì)算結(jié)果更利于表達(dá)年際變化規(guī)律和季節(jié)性變化規(guī)律,但站點(diǎn)分布數(shù)量少。InSAR監(jiān)測(cè)優(yōu)勢(shì)主要體現(xiàn)在面狀地面沉降信息的提取,大數(shù)據(jù)量更有利于刻畫空間分布的細(xì)節(jié)信息,但在單點(diǎn)時(shí)間序列擬合上數(shù)量相對(duì)較少。

4 結(jié) 語(yǔ)

本文以我國(guó)重點(diǎn)沉降區(qū)(華北平原)為研究區(qū),分析多軌道時(shí)序InSAR誤差源,綜合考慮水準(zhǔn)基巖點(diǎn)、InSAR形變參考點(diǎn)、CORS站等因素選取擬穩(wěn)點(diǎn),提出采用擬穩(wěn)基準(zhǔn)平差方法進(jìn)行多軌道時(shí)序InSAR平差計(jì)算,實(shí)現(xiàn)大范圍多軌道InSAR垂直形變場(chǎng)基準(zhǔn)的統(tǒng)一。研究結(jié)果表明,經(jīng)過(guò)擬穩(wěn)基準(zhǔn)平差計(jì)算后,研究區(qū)各條帶之間沉降速率過(guò)渡自然,無(wú)明顯分界,擬穩(wěn)基準(zhǔn)平差后內(nèi)符合精度統(tǒng)計(jì)結(jié)果明顯優(yōu)于初始數(shù)據(jù)處理結(jié)果,擬穩(wěn)基準(zhǔn)平差后平均誤差優(yōu)于5 mm,中誤差優(yōu)于9 mm,說(shuō)明時(shí)序InSAR擬穩(wěn)基準(zhǔn)平差結(jié)果內(nèi)符合精度良好。對(duì)比分析研究區(qū)CORS站垂直方向變化量與InSAR結(jié)果發(fā)現(xiàn),兩者吻合程度較好,從外符合精度角度說(shuō)明時(shí)序InSAR擬穩(wěn)基準(zhǔn)平差結(jié)果具有可靠性。

致謝：感謝中國(guó)地震局地殼形變臺(tái)網(wǎng)中心提供GNSS基準(zhǔn)站位移結(jié)果。