融合分布式散射體時序InSAR技術(shù)在礦區(qū)形變調(diào)查中的應(yīng)用

賈會會, 張海清, 李克達，張小朋

1.華北地質(zhì)勘查局五一四地質(zhì)大隊，河北 承德 067000 2.河北省地質(zhì)災(zāi)害監(jiān)測預(yù)警技術(shù)創(chuàng)新中心，河北 承德 067000 3.河北大學建筑工程學院，河北 保定 071002

0 引言

河北省承德市灤平縣張百灣鎮(zhèn)周臺子村鐵礦區(qū)由于多年礦山開采導致歷史遺留大量的淺埋采空區(qū)，其中部分采空區(qū)已發(fā)生明顯的地面塌陷。其中：20170814，灤平縣啟星礦業(yè)集團啟泰采區(qū)華兵選礦廠發(fā)生采空塌陷事件，形成深約40 m的塌陷坑，致使華兵選礦廠全部陷入地下，造成巨大的經(jīng)濟損失；201904—201905，距離周臺子村約90 m的三采區(qū)連續(xù)發(fā)生兩次采空區(qū)地面塌陷，塌陷坑深度約40 m，直徑約50 m，造成礦路損毀斷交。潛在的塌陷安全隱患會對周邊建筑和居民區(qū)造成巨大威脅，因此對礦區(qū)進行形變調(diào)查和監(jiān)測，對采空區(qū)塌陷進行災(zāi)害預(yù)警十分必要。

傳統(tǒng)礦區(qū)形變監(jiān)測的方法有水準測量、GPS測量、三維激光掃描[1-3]等。這些以“點”為基礎(chǔ)的形變監(jiān)測技術(shù)耗費大量人力、物力、財力，且監(jiān)測結(jié)果空間和時間分辨率較低，因此無法進行大范圍、高密度的礦區(qū)地表形變監(jiān)測[4]。近幾十年發(fā)展起來的InSAR(合成孔徑雷達干涉測量技術(shù))[5]具有監(jiān)測精度高、全天時全天候、監(jiān)測結(jié)果時間分辨率和空間分辨率較高等優(yōu)勢，其中監(jiān)測精度可達到厘米甚至毫米級，為礦區(qū)形變監(jiān)測和范圍調(diào)查提供了新的有效監(jiān)測手段[6-12]。傳統(tǒng)時序InSAR技術(shù)采用永久散射體來獲取地表形變信息，但這種方法的缺陷是自然地表測量點較少、空間采樣率低，且很難確定出鐵礦采空區(qū)的沉降范圍和形變趨勢。針對礦區(qū)內(nèi)監(jiān)測點數(shù)量不足的問題，本文采用DS(distributed scatterer，分布式散射體)和高相干點時間序列InSAR技術(shù)來對礦區(qū)進行形變解算。DS點是指分辨單元內(nèi)具有相同散射特性的散射體，例如沙漠、草甸等。DS點處理流程包括同質(zhì)像素點的提取和最優(yōu)相位估計。提取同質(zhì)像素點的現(xiàn)有方法是利用配準后的時間序列SAR(合成孔徑雷達)振幅影像進行t檢驗、KS(Kolmogorov-Smirnov)檢驗、AD(Anderson-Darling)檢驗等[13]。這些方法雖然在效率和魯棒性方面存在不足，但大多數(shù)都能有效地解決同質(zhì)像素點的識別問題?，F(xiàn)有的最優(yōu)相位估計方法主要有MLE(極大似然估計)、PL(phase linking, 相位鏈接)算法和PT(phase triangulation,相位三角)算法[14-16]等。然而，這些方法不能在一個分辨單元中分解多個散射機制，否則可能會影響同質(zhì)像素點相位的估計精度。由于中低分辨率SAR數(shù)據(jù)的分辨單元包括很多不同類型的地面目標；因此，可采用樣本協(xié)方差矩陣或者相干矩陣EVD(特征值分解)方法進行單個分辨率單元中多個散射機制的區(qū)分，并進一步提高相位優(yōu)化精度。此外，可基于分布式散射體的時間序列InSAR技術(shù)獲得高密度監(jiān)測點，提取鐵礦開采區(qū)引起的地表形變不同時空特征，進而對礦區(qū)形變進行精確識別和監(jiān)測。

為查明灤平縣周臺子村鐵礦采空區(qū)空間分布，本文采取融合分布式散射體和相干散射體(distribute scatterer-coherent scatterer InSAR,DS-CSInSAR)的方法對研究區(qū)地表形變進行監(jiān)測。即先采用AD檢驗的方法提取同質(zhì)像素點，再采用基于協(xié)方差矩陣特征值分解算法來對DS點的最優(yōu)相位進行估計，然后通過融合相干點來對礦區(qū)地表形變進行反演，并確定礦區(qū)形變的空間分布；以期為采空區(qū)塌陷災(zāi)害預(yù)警及治理提供技術(shù)支撐。

1 研究區(qū)和數(shù)據(jù)集

1.1 研究區(qū)介紹

研究區(qū)位于河北省承德市灤平縣周臺子村，鐵礦采區(qū)位于灤河北岸。研究區(qū)的地理位置見圖1a，范圍為117.447°E—117.678°E，40.867°N—41.082°N，紅色邊框為研究區(qū)覆蓋范圍。該區(qū)地貌類型為低山侵蝕構(gòu)造地貌(圖1b)，高程范圍為269～1 650 m；氣候?qū)俅箨懠撅L氣候，冬長寒冷，夏短而炎熱，多年平均氣溫9.1 ℃，最熱月(7月)平均氣溫24.4 ℃，最冷月(1月)平均氣溫-9.4 ℃，極端最高氣溫41.5 ℃，極端最低氣溫-24.2 ℃，最大日溫差23.8 ℃；歷年最大降水量835.9 mm，最小降水量326.7 mm，平均降水量557.9 mm。研究區(qū)GPS監(jiān)測點分別位于華兵礦業(yè)、JC33廣場、JC34二采、JC35路邊(圖1c，d)，A—J鐵礦區(qū)的位置見圖1d。

a.研究區(qū)地理位置；b.研究區(qū)地質(zhì)地貌圖；c.研究區(qū)數(shù)字高程模型(DEM)圖；d.鐵礦和GPS監(jiān)測點位置圖。a圖中，藍色邊框為Sentinel-1數(shù)據(jù)的覆蓋范圍，紅色邊框為研究區(qū)；b圖中，紅色邊框為研究區(qū)；c圖中，紅色圓圈為該地區(qū)GPS監(jiān)測點；d圖中，A—J為鐵礦區(qū)位置。

1.2 數(shù)據(jù)集

本文采用C波段48景250 km幅寬的Sentinel-1數(shù)據(jù)來對調(diào)查區(qū)進行形變反演。Sentinel-1數(shù)據(jù)集時間跨度為20170316—20201002，IW(interfe-rometric wide swath)TOPS(terrain observation with progressive scans，漸進掃描地形觀測)模式，升軌軌道號為69，距離向和方位向分辨率分別為2.3 m和13.9 m。48景SAR圖像的Sentinel-1數(shù)據(jù)時空基線圖見圖2，設(shè)置時間基線和空間基線的閾值分別為120 m和80 d，同時剔除一些受大氣影響嚴重的干涉對，最終得到88個干涉對，其中最大的垂直基線和最大的時間基線分別為115.435 m和60 d。

圖2 Sentinel-1數(shù)據(jù)時空基線圖

2 技術(shù)方法

為了更精確地識別和監(jiān)測調(diào)查區(qū)，首先提取出DS點并與相干點進行融合，得到研究區(qū)更精確的地表形變信息；然后使用沉降信息進行礦區(qū)的時空分析及沉降影響因素分析；最后確定出礦區(qū)的沉降區(qū)域。具體技術(shù)流程見圖3，主要包括Sentinel-1 TOPS模式預(yù)處理、相干點和DS點處理、時序形變解算等。

γ.平均相干系數(shù)；γt.時間相干性；γth.平均相干系數(shù)的固定閾值。

2.1 Sentinel-1 TOPS模式預(yù)處理

進行預(yù)處理時，首先采用30 m分辨率的SRTM(shuttle radar topography mission，航天飛機雷達地形測繪使命)數(shù)據(jù)進行DEM(digital elevation model，數(shù)字高程模型)配準；其次，由于TOPS成像模式中多普勒中心頻率的變化，因此需要采用ESD(enhanced spectral diversity，增強譜分集)配準方法來減少Burst間的相位跳變，并進一步優(yōu)化方位向的偏移；然后，在將輔圖像重新采樣到主圖像的坐標系中之前，需要對輔圖像進行降噪去斜，同時對每個Burst都進行去斜、插值和反去斜的流程[17]，從而使SAR圖像方位向上的配準精度優(yōu)于1/1000像素[18]；之后，將每個條帶的Burst進行拼接，因為研究區(qū)位于第2個條帶，因此不需要進行條帶拼接；最后按照研究區(qū)的范圍進行ROI(感興趣區(qū))的提取。

2.2 相干點和DS點處理

本文利用研究區(qū)時間序列的相干系數(shù)圖得到平均相干系數(shù)圖，利用平均相干系數(shù)閾值法，即γ≥γth，提取平均相干系數(shù)γ在固定閾值(γth=0.7)以上的相干點作為高相干點。為了提取研究區(qū)的DS點，利用AD檢驗提取同質(zhì)像素點，并利用特征值分解算法對同質(zhì)像素點進行最優(yōu)相位估計。

2.2.1 AD檢驗

相對于KS檢驗對尾部差異不敏感的缺點，AD檢驗在尾部增加了權(quán)因子，以有效減少第二類錯誤[19]。對配準后時間序列的SAR振幅圖像，其中任意兩個像素a和b的AD檢驗統(tǒng)計參數(shù)A2可表示為

(1)

式中：N為SAR圖像個數(shù)；Fa(x)和Fb(x)分別為兩個像素a和b的振幅的CDF(經(jīng)驗累積分布函數(shù))；Fa,b(x)為a和b兩像素點的CDF。當A2滿足特定的閾值，則a和b為SHP(statistically homogeneous pixels，同質(zhì)像素點)。

為了提取同質(zhì)像素點，取13×13窗口作為本文的搜索窗口，使用上述的AD檢驗去比較該窗口內(nèi)的中心像素點和其他像素點的檢驗統(tǒng)計量。本文中顯著性水平取0.95，同質(zhì)像素點閾值取20。當初始的同質(zhì)像素點提取后，先判斷其與中心像素點是否連接，舍棄與中心像素點沒有連接的像素點，則剩下的像素點即為初始的DS點集。

2.2.2 同質(zhì)像素點的最優(yōu)相位估計

最優(yōu)相位估計是SHP識別的關(guān)鍵步驟，本文基于協(xié)方差矩陣特征值分解算法對最優(yōu)相位進行估計，該方法優(yōu)勢是不用進行相干矩陣求逆操作，避免了迭代耗時且計算效率高，且該方法適用于中低分辨單元內(nèi)的多種散射機制[20]。

上述提取DS點的協(xié)方差矩陣C見下式：

(2)

式中：E為數(shù)學期望；p為時間序列像素x的復(fù)數(shù)向量，即p=[p1,p2,p3,…,pN]?；?為厄米特共軛轉(zhuǎn)置；k為同質(zhì)像素點集合的第k個點；Ω為Nx個同質(zhì)像素點的集合。

因此基于協(xié)方差矩陣的特征分解可見下式：

(3)

(4)

由于分辨單元內(nèi)主導散射機制數(shù)量n沒有建立先驗知識且不同特征值對應(yīng)的特征向量相互正交；因此，求解最大特征值λi對應(yīng)的特征向量vi的相位分量作為優(yōu)化相位的估計值。

2.2.3 DS點與相干點融合

根據(jù)時間相干性γt對DS點進一步篩選，先將平均相干系數(shù)γ大于γt的DS點提取出來，即

(5)

式中：M為最優(yōu)干涉相位值數(shù)目；φjk為原始干涉相位；θ為時間域濾波后求得的最優(yōu)相位，j和k分別表示第j個和第k個干涉對;i為虛數(shù)單位。

然后，將上述提取的高相干點與DS點進行融合，同時將DS點最優(yōu)相位估計與相干點的原始相位進行融合，進行時序形變解算。

2.3 時序形變解算

將DS點和相干點融合后，對這些點集構(gòu)建Delanuary三角網(wǎng)。步驟如下：首先，構(gòu)建線性形變模型進行鄰近測量點的干涉相位差分，對三角網(wǎng)中每條邊進行相對形變速率和相對高程誤差解算，采用空間搜索方法對三角網(wǎng)中每條邊進行求解，邊的閾值取0.6，小于該閾值的邊則被剔除，不參與形變解算，得到三角網(wǎng)每條邊的相對形變速率和高程誤差。然后，選擇一個參考點，建立絕對形變速率與相對形變速率的系數(shù)矩陣，采用最小二乘算法對絕對形變速率和絕對高程誤差進行求解。再將求出的形變速率和高程誤差從干涉相位中減去，即可得到每個點的殘余相位。最后，對殘余相位進行MCF(minimum cost flow，最小費用流)相位解纏，得到其真實相位。殘余相位包括大氣相位、非線性形變相位和噪聲相位，利用大氣相位在時域上表現(xiàn)為高頻信號、在空域上表現(xiàn)為低頻信號的特性，采用時間域的高頻和空間域的低頻濾波對大氣相位進行去除并估計出大氣相位，將上述求得線性形變與非線性形變相加即可得到研究區(qū)的時序形變量結(jié)果。

3 結(jié)果與分析

3.1 相干點和DS點預(yù)處理結(jié)果

同質(zhì)像素點的提取結(jié)果能夠反映出同質(zhì)區(qū)域同質(zhì)點的分布情況，采用AD檢驗識別的同質(zhì)像素點集見圖4a。該區(qū)域大部分同質(zhì)像素點的數(shù)目在100～255之間。圖4b為利用相干性系數(shù)大于0.7篩選出來的高相干點，可見研究區(qū)由于位于自然地表區(qū)域，相干點非常少，僅為3 514個。圖4c為利用時間相干性閾值為0.5篩選的最終的DS點，數(shù)量為57 766，圖4d為DS點和相干點融合的結(jié)果?？梢娊?jīng)過DS點加密，測量點的數(shù)量得到極大的提升，提高了時序形變解算的精度，且減小了相位解纏的誤差。

a.同質(zhì)像素點集分布圖；b.高相干點目標分布圖；c.DS點分布圖；d.DS點與相干點融合圖。b、c、d圖中，紅色點代表選出的點。

利用DS點所有干涉對的相位信息，采用基于協(xié)方差矩陣特征值分解算法對單一主圖像干涉對的最優(yōu)相位進行估計，以20181013—20181130和20181130—20181224兩個干涉對為例，結(jié)果見圖5。從圖5a、5c可看出原始干涉相位去相干噪聲較嚴重；對分布式散射體進行相位優(yōu)化后，相位質(zhì)量和整體的相干性得到明顯的提高，在低相干地區(qū)，干涉相位的相干性也得到提高(圖5b、5d)。

a.20181013—20181130原始相位；b.20181013—20181130最優(yōu)相位；c.20181130—20181224原始相位；d.20181130—20181224最優(yōu)相位。

3.2 鐵礦區(qū)形變結(jié)果

通過本文的融合分布式和相干散射體時序InSAR處理后，得到鐵礦區(qū)年平均形變速率(圖6)，可知研究區(qū)的形變速率區(qū)間為-34.50～24.50 mm/a，形變趨勢較大的地區(qū)主要集中在圖6紅色矩形框內(nèi)。提取圖6中紅色邊框的形變速率數(shù)據(jù)，結(jié)果見圖7a。由該地區(qū)A—J鐵礦區(qū)域紅色曲線范圍(圖7b、c)可發(fā)現(xiàn)，這些地區(qū)的形變速率較大，且均出現(xiàn)不同程度的沉降。提取圖7b和圖7c中A—J鐵礦區(qū)域中點目標(白色圓圈)進一步分析時序形變位移。圖8為A—J形變區(qū)域提取的點目標時序形變位移曲線圖，可發(fā)現(xiàn)：A—D區(qū)域在每年6—9月出現(xiàn)不同程度的沉降，沉降幅度在-33.40～32.70 mm，其中負值表示下沉，正值表示抬升；E區(qū)域在20200827沉降量最大，為50.76 mm；F—I區(qū)域在每年6—9月形變量變化區(qū)間為-35.02～26.03 mm，其中周臺子村附近礦區(qū)F、G區(qū)域的最大沉降量分別為32.00 mm和34.00 mm，離路邊最近的礦區(qū)H區(qū)域最大沉降量為33.00 mm；在J區(qū)域窯嶺溝礦區(qū)出現(xiàn)明顯的抬升現(xiàn)象，最大抬升量為24.03 mm。

紅色矩形框內(nèi)代表形變趨勢較大的地區(qū)。

a.年平均形變速率圖；b.鐵礦區(qū)A—G形變區(qū)域；c.鐵礦區(qū)H—J形變區(qū)域。

3.3 精度驗證

為驗證本文時序InSAR結(jié)果的準確性和精度，本研究獲取灤平縣礦區(qū)華兵礦業(yè)、JC33廣場、JC34二采、JC35路邊4個GPS監(jiān)測點20200715—20201002期間24 h間隔的GPS實測數(shù)據(jù)和圖8中GPS點對應(yīng)的InSAR時序測量值進行對比驗證。二者InSAR時序形變量與GPS的結(jié)果對比見圖9，進而可計算出4個監(jiān)測點與InSAR測量值之間的絕對誤差分別為-1.70～3.70，-0.80～4.20，-3.42～3.77，-2.86～4.79 mm。InSAR與GPS監(jiān)測相關(guān)性對比結(jié)果(圖10)表明：JC33廣場和JC35路邊的相關(guān)系數(shù)分別為0.93和0.81，這是由于這2個GPS監(jiān)測點位于建筑物附近，形變趨勢穩(wěn)定，且InSAR干涉圖的相干性較好，在進行形變反演時受到的大氣、地形誤差影響較小，InSAR形變精度較高，因此與GPS測量值有著較高的相關(guān)性；華兵礦業(yè)和JC34二采相關(guān)系數(shù)分別為0.51和0.65，這是由于這2個GPS監(jiān)測點位于礦區(qū)附近，屬于自然地表，在InSAR處理過程中，易受到去相干的影響，形變解算精度沒有建筑區(qū)的高。此外InSAR測量值與入射角和LOS(line-of-sight，雷達視線方向)、殘余大氣相位也有關(guān)系，因此與GPS測量值有著較低的相關(guān)性，但是它們的絕對誤差范圍較小，從而證明了InSAR監(jiān)測灤平縣礦區(qū)的地表形變結(jié)果與GPS野外觀測結(jié)果一致性較好，結(jié)果可信度較高。

a.A—E區(qū)域；b.F—J區(qū)域。

a.華兵礦業(yè)；b.JC33廣場；c.JC34二采；d.JC35路邊。

a.華兵礦業(yè)；b.JC33廣場；c.JC34二采；d.JC35路邊。

為了對比分析DS-CSInSAR 技術(shù)結(jié)果的準確性，我們采用CSInSAR技術(shù)對該鐵礦區(qū)進行了PSI(永久散射體雷達干涉)分析、Delaunay三角網(wǎng)構(gòu)建、相對形變和絕對形變解算等步驟得到研究區(qū)的年平均形變速率，結(jié)果見圖11。圖11表明：兩種方法得到的形變趨勢基本一致，但與 CSInSAR 技術(shù)對比，本文提出的DS-CSInSAR方法極大地提高了點的密度；通過 DS 的空域濾波和最優(yōu)相位估計，DS-CSInSAR方法提高了低相干鐵礦區(qū)干涉圖相位的質(zhì)量和相干性，從而間接證明了該算法的可靠性。

圖11 CSInSAR技術(shù)(a)與DS-CSInSAR技術(shù)(b)結(jié)果對比圖

4 結(jié)論與建議

1)為彌補礦區(qū)監(jiān)測點目標不足的情況，本文提出了DS-CSInSAR技術(shù)的鐵礦區(qū)識別與監(jiān)測的方法，利用雙樣本AD檢驗和協(xié)方差矩陣特征值分解方法進行相位優(yōu)化，大大提高了礦區(qū)的測量點數(shù)量和形變解算的精度。與CSInSAR技術(shù)方法相比，DS-CSInSAR技術(shù)方法在監(jiān)測點密度和干涉圖的質(zhì)量上均有明顯改善。

2)通過對研究區(qū)地表形變的時空分析發(fā)現(xiàn)，研究區(qū)年平均形變速率為-34.50～24.50 mm/a。不同礦區(qū)在每年的6—9月出現(xiàn)不同程度的沉降，其中周臺子村和路邊附近F、G、H區(qū)域的最大沉降量分別為32.00、34.00、33.00 mm；窯嶺溝礦區(qū)出現(xiàn)明顯抬升現(xiàn)象，最大抬升量為24.03 mm。InSAR 監(jiān)測結(jié)果與4個GPS站點結(jié)果一致，相關(guān)性較好，絕對誤差為-3.42～4.79 mm。

3)本文利用高密度的DS點確定了鐵礦區(qū)形變趨勢，相對于其他技術(shù)具有空間分辨率高、監(jiān)測范圍大的優(yōu)勢。該方法較好地解決了在地形復(fù)雜鐵礦區(qū)進行形變監(jiān)測的問題，進一步提高了利用InSAR形變結(jié)果識別和監(jiān)測鐵礦區(qū)的能力；同時，根據(jù)沉降結(jié)果進行的分析也為預(yù)測鐵礦區(qū)形變強度和位移趨勢提供了重要信息。