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

1 自然資源部第一大地測量隊,西安市測繪路4號,710054

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

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

1 多軌道時序InSAR誤差源分析

1.1 入射角引起的差異

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

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

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

1.2 參考點引起的差異

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

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

1.3 軌道殘余誤差

軌道殘余誤差主要由軌道估計的偏差所致,在數(shù)據(jù)處理結(jié)果上體現(xiàn)為趨勢性偏差,評估軌道殘余誤差的數(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

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

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

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

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

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

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

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

(1)

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

(2)

(3)

(4)

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

(5)

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

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

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

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

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

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

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

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

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

3.3 精度評定

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

分別對初始數(shù)據(jù)處理結(jié)果和擬穩(wěn)基準(zhǔn)平差結(jié)果按照上述方法統(tǒng)計重疊區(qū)域PS點速率差值的平均誤差和中誤差,結(jié)果如表3(單位mm/a)、表4(單位mm/a)所示。從表中可以看出,通過擬穩(wěn)基準(zhǔn)平差后內(nèi)符合精度統(tǒng)計結(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)計Tab.3 Statistics of internal coincidence accuracy of initial data processing results

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

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

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

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

4 結(jié) 語

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

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