多軌道時序InSAR擬穩基準平差方法研究

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

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

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

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

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

1.1 入射角引起的差異

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

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

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

1.2 參考點引起的差異

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

本文數據處理選擇多參考點,將綜合分析選取的擬穩點作為多參考點進行數據處理。為保證各條帶形變參考基準的統一,相鄰條帶間至少保證重疊2個擬穩點。

1.3 軌道殘余誤差

軌道殘余誤差主要由軌道估計的偏差所致,在數據處理結果上體現為趨勢性偏差,評估軌道殘余誤差的數學模型可表示為:

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

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

1.4 大尺度地形相關性誤差

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

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

2 擬穩基準平差方法

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

設不穩定點的未知數為X1,穩定點的未知數為X2,則誤差方程為[12]:

(1)

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

(2)

(3)

(4)

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

(5)

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

3 實驗數據與結果分析

3.1 數據源

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

表1 SAR數據Tab.1 SAR data

3.2 處理結果

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

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

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

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

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

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

依據InSAR監測結果,采用克里格方法對InSAR相干目標點平均沉降速率進行插值計算,表2為平均沉降速率分類統計結果。

表2 InSAR結果平均沉降速率分類統計Tab.2 Statistics of average subsidence rate of InSAR results

3.3 精度評定

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

分別對初始數據處理結果和擬穩基準平差結果按照上述方法統計重疊區域PS點速率差值的平均誤差和中誤差,結果如表3(單位mm/a)、表4(單位mm/a)所示。從表中可以看出,通過擬穩基準平差后內符合精度統計結果明顯優于初始數據處理結果。由圖5(a)可知,擬穩基準平差后平均誤差優于5 mm,中誤差優于9 mm。

表3 初始數據處理結果內符合精度統計Tab.3 Statistics of internal coincidence accuracy of initial data processing results

表4 擬穩基準平差后內符合精度統計Tab.4 Statistics of internal coincidence accuracy after quasi-stable datum adjustment

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

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

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

4 結 語

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

致謝：感謝中國地震局地殼形變臺網中心提供GNSS基準站位移結果。