高建東,王 勇,3,安江雷,姜俊狄
(1. 江蘇省測繪工程院,江蘇 南京 210013; 2. 自然資源部國土衛星遙感應用重點實驗室,江蘇 南京 210013; 3. 武漢大學衛星導航技術研究中心,湖北 武漢 430079; 4. 河海大學地球科學與工程學院,江蘇 南京 211100)
地面沉降是指在自然或人為因素作用下,地表土層緩慢緊縮而導致區域性或局部地面下沉的地質現象,嚴重時可造成不可逆轉的自然資源和生態環境損失[1]。測繪上常規的地面沉降監測手段主要有精密水準監測、GNSS及基巖標分層標監測等方法,大多通過監測點測量結果的變化值反演地面沉降量。雖然這些技術手段單點監測精度高,但是需要根據先驗的區域地面沉降信息,在合適位置布設監測點。而且外業觀測工作量大、周期長、成本高,監測點在空間上一般都呈現出離散且稀疏的分布特點,只能反映監測點區域的地面沉降信息,無法獲取大范圍沉降細節,降低了數據的可用性[2]。
InSAR技術的引入彌補了精密水準和GNSS等傳統測手段在監測范圍上的不足。但是InSAR變形監測也面臨諸多挑戰,比如,其只能監測地面形變在雷達視線方向(line-of-sight, LOS)上的一維投影;對于低相干區,很難捕捉到有效形變信息;InSAR的形變監測結果基本都依賴于實地測量數據(如水準、GNSS等)進行標定[3]。
因此,本文綜合利用精密水準測量、GNSS技術、InSAR技術3種方法的優勢進行聯合監測,以期獲取更為精確全面的區域地表形變信息[4-5]。
(1)時間基準統一。一般而言,精密水準、GNSS和時序InSAR 3種監測手段的監測周期不會完全一致,因此在進一步數據處理前,需要對精密水準、GNSS和時序InSAR監測數據進行時間基準的統一,即將3種監測手段獲取的監測數據統歸算到相同的時間間隔上[6]。監測所得的形變值可作為時間的一元函數,一般首先采用插值法、多項式法進行擬合逼近,然后截取相同時間段的數據進行分析[7],統一時間基準。
(2)空間基準統一。由于精密水準進行地面沉降監測的基準面是似大地水準面,采用正常高系統,基于GNSS技術的監測是基于橢球面,監測成果屬于大地高系統,而時序InSAR監測處理的結果一般為雷達視線方向(LOS)的變化值[8]。要進行3種監測技術的沉降監測數據融合,必須將其統一到相同的空間基準上。
采用統計學中常用的Bland-Altman(B-A)法[9]進行監測數據的一致性評價與分析。該方法是定量與定性方法的結合,其原理是對于兩種評定間的差異進行隨機效應分析,來解釋說明一致性問題。該方法假定測量誤差不會影響變量間的相關性,但會影響一致性,其實質是對兩種評定之間差異的一種統計分析,即利用統計學原理定量進行兩種評定間的系統差異,進而依此進行一致性評價。
(1)重合點位數據融合方法。重合點即兼有精密水準、GNSS、InSAR監測結果的點位。經監測數據時空基準統一和一致性分析后,篩選出一致性較好的數據,直接采用最小二乘法進行數據融合。
(2)非重合點位數據融合方法。考慮監測成本、點位環境要求等因素,精密水準和GNSS監測點數量一般較少。因此,兼具3種監測方法成果的點位不會太多。考慮精密水準和GNSS監測手段精度高,經過重合點融合計算后能提高時序InSAR監測結果的精度和可靠性。參考控制測量思路,以重合點位改正值為控制框架,對其余非重合點位進行插值改正。
利用m種監測方法在某一監測首期內對同一研究區域進行地面沉降監測,根據前述分析,假設第i個重合點位分別獲取了地面沉降量di1,di2,…,dim,融合計算后的地面沉降量為di,則融合后的地面沉降量與各種監測手段的差值分別為Δdi1,Δdi2,…,Δdim,容易得:Δdi1=di-di1,Δdi2=di-di2,…,Δdim=di-dim。
基于重合點位在同一坐標框架下的坐標值和監測值差值,采用多面函數擬合的方法即可得到研究區域整體的分布情況。公式為
(1)
式中,Δd為擬合某種監測手段監測值與融合后的地面沉降量之間的差值;aj為多面函數的系數;F(x,y,xi,yi)為x、y的二次函數,其核心點位于(xi,yi)處,常用的二次核函數為
(2)
選擇m個已知點為核心點,令
Qij=F(xi,yi,xj,yj)
(3)
列出誤差方程為
v=Qa-D
(4)
式中,

(5)
在vTPv=min的前提下,可求得系數a為

(6)
將a代入原多面函數中,并利用這些點位的坐標,計算出區域內某種監測手段任意點位處的監測值相對于融合后數據的差值Δdi。這些差值加上原監測值即為融合后的沉降量hi,即
hi=di+Δdi
(7)
通過融合精密水準、GNSS和時序InSAR等手段獲取的研究區域沉降值,可得到研究區域高精度的面狀地面沉降信息。
某沉降區位于江蘇省西北部,包含徐州市沛縣和豐縣區域,如圖1所示。該區域北靠山東省濟寧市微山縣和魚臺縣,西側與安徽省碭山縣相接,東側毗鄰微山湖和邵陽湖,南側為徐州市銅山區。該區域富含煤礦資源,煤田總面積約160 km2[10]。

圖1 研究區位置
2.2.1 精密水準數據
使用數據為研究區域精密水準沉降監測網(如圖2所示)兩期觀測數據,按照國家一、二等水準測量規范施測,屬于較高精度的精密水準測量數據[11]。該地面沉降監測網第1期觀測時間為2019年7月23日至30日,第2期觀測時間為2020年10月10日至10月17日,兩期觀測采用相同的觀人員、儀器、路線。平差計算精度統計見表1。

表1 地面沉降監測水準網平差計算精度統計 mm

圖2 研究區域水準點及水準路線分布
2.2.2 GNSS數據
GNSS監測數據來源于徐州市GNSS地面沉降監測網(如圖3所示)數據,按照B級GNSS控制網精度施測,第1期靜態GNSS監測網觀測時間從2019年7月26日至8月3日,第2期觀測時間從9月24日至10月1日。基線解算采用GAMIT軟件進行,平差計算采用武漢大學研制的PowerNET軟件。計算結果精度統計見表2。

表2 GNSS監測網平差計算精度統計 mm

圖3 GNSS監測網點分布
由于精密水準和GNSS均為點狀監測,選取既有精密水準監測結果又有GNSS監測結果的16個監測點,將兩期監測點的垂直方向變化速率作為地面沉降速率進行插值繪圖,如圖4所示。

圖4 基于精密水準測量和GNSS測量的地面沉降速率分布
2.2.3 時序SAR數據處理
采用的數據為Sentinel-1A衛星Level-1級產品:干涉寬幅(IW)模式的單視復數(SLC)的影像,每景影像包含3個Subswath子帶,每個Subswath子帶由9~10個Burst組成,覆蓋范圍在南北方向約180 km,東西方向約250 km,重訪周期為12 d,分辨率為5 m×20 m,時間跨度為2019年7月21日至2020年10月8日。
采用SBAS-InSAR技術[12]進行時序InSAR數據處理。該方法是一種基于多個主影像的時序InSAR方法,克服了單一主影像導致干涉結果較差的缺點,提高了干涉圖的相干性和運算效率[13],所需的SAR影像數據較少,在植被覆蓋地區的處理效果更好,處理結果如圖5所示。

圖5 時序InSAR地面沉降速率分布
2.2.4 融合處理
以時序InSAR監測數據時間為基準,取多源地面沉降監測手段數據采集時間的交集,將精密水準和GNSS方式獲取的沉降值歸算到2019年7月21日至2020年8月8日這一時間周期上3種監測手段得到的監測數值。選取了16個兼有精密水準、GNSS和時序InSAR監測結果的點位進行融合分析。
對精密水準監測結果和GNSS方法監測結果進行一致性檢驗,兩種監測方法獲取的累計沉降量如圖6所示。

圖6 精密水準和GNSS方法獲取的累計沉降量
在時序InSAR監測結果中選取和水準點位置一致的同名點進行分析,若水準點處沒有時序InSAR監測結果,則采取鄰近點加權平均的方法進行同名點提取。兩種監測方法獲取的累計沉降量如圖7所示。

圖7 精密水準和時序InSAR獲取的累計沉降量
采用B-A一致性分析方法中的差值比較法進行分析計算。精密水準測量監測數據與GNSS監測數據的一致性為93.8%,精密水準測量監測數據與時序InSAR監測數據的一致性為95.2%。
采用最小二乘法進行數據融合。按照1.1.2節所述,對16個重合點位融合后的地面沉降值與單一監測手段獲取的地面沉降值求差。基于重合點位在同一坐標框架下的坐標值和監測值差值,采用多面函數擬合的方法即可得到研究區域整體的分布情況,進行插值繪圖,如圖8所示。可以看出,相較于精密水準和GNSS檢測,融合后的結果能反映大面積區域的沉降信息;相較于單一時序InSAR檢測,附加精密水準和GNSS監測數據改正后,監測結果的精度和可靠性得到了提升,融合后的地面沉降速率圖比單一數據表現更加豐富,細節更加清晰。
2.3.1 研究區域地面沉降時空特征分析
將離散監測點通過克里金差值后疊加行政村界,如圖9所示。結果表明:研究區域主要沉降區域位于沛縣安國鎮、龍固鎮、朱寨鎮、大屯街道及豐縣分鳳城街道,在監測周期內的最大沉降速率為29.8 mm/a,位于沛縣安國鎮冠英村。監測周期內,地面沉降速率超過10 mm/a的沉降區面積為2 018.8 km2,占研究區域總面積61.9%。
2.3.2 采礦區域地面沉降分析
徐州煤炭開采已有130多年的歷史,為我國經濟發展提供了大量煤炭資源,而沛縣北部礦區作為徐州煤炭資源的主體,包括姚橋煤礦在內的一大批高產量煤田,該區域探明的煤炭儲量超過20億t[14]。除此之外,沛縣還有姜梨園和新莊兩處鐵礦開采。而豐縣區域有李堂煤礦、瑞豐鹽礦和金馬鹽礦等3處采礦區域。如圖10所示。

圖10 多源監測數據融合后地面沉降速率分布(疊加采礦區)
由圖10可以看出,除閉礦較早(2013年)的沛城煤礦區域外,其他礦區均存在較嚴重的大面積地面沉降現象,地面沉降現象的發生與礦區的分布呈明顯的正相關關系。長期的煤礦和其他礦產開采造成了嚴重的生態惡化與環境污染。一方面,隨著煤炭資源的持續開采,地面標高降低、耕地積水、建構筑物損害、土地鹽漬化等生態環境問題,直接影響到礦群關系及煤礦的正常生產。另一方面,礦井關閉后,隨著地下水位上升,采動覆巖結構受地下水作用,將產生二次移動變形,從而導致地表移動變形、土地塌陷、建構筑損害等災害。近年來,隨著國家對自然資源和生態環境保護的重視,對礦產開采和地下水抽取進行了整頓治理[15],地面沉降現象得以一定程度遏制。
本文采用源地面沉降監測數據融合方法,綜合利用精密水準測量、GNSS技術、InSAR技術3種方法的優勢進行聯合監測,獲取了更為精確全面的區域地面沉降信息。但也存在計算模型較為簡單的不足之處,對范圍較大、情況復雜區域適用性還有待進一步提升,需探索研究更加科學嚴密的插值模型,以獲取更為準確的地沉降監測數據融合結果。
隨著地面沉降監測研究的深入,通過單一化的技術手段和測量方法往往難以準確把握區域地面沉降全貌,需建立天地一體化地面沉降動態監測體系,融合精密水準監測、GNSS監測、時序InSAR監測、基巖標分層標志監測、地下水變化監測、陣列式微震監測和光纖監測等多種技術手段建立自動化監測系統,開展全周期、多視角、高精度、系統性的監測工作,揭示地面沉降現象全貌,準確探究地面沉降分布和發育規律。