基于多閾值目標(biāo)提取的時(shí)序InSAR礦區(qū)地表沉降監(jiān)測研究

劉澤洲，盧才武，章 賽，李 萌，和鄭翔

(1.西安建筑科技大學(xué)資源工程學(xué)院，陜西 西安 710055；2.西安建筑科技大學(xué)建筑學(xué)院，陜西 西安 710055)

礦區(qū)資源長期地下開采會導(dǎo)致地表沉降和形變，從而引發(fā)地表沉陷、裂縫以及滑坡等地質(zhì)災(zāi)害，不但帶來嚴(yán)重的經(jīng)濟(jì)損失,而且還會危害居民正常生活及人身安全。為了有效控制礦區(qū)災(zāi)害的發(fā)生，必須通過監(jiān)測手段持續(xù)準(zhǔn)確地對礦區(qū)地表變化情況進(jìn)行監(jiān)測，全面并及時(shí)地獲取地表沉降情況。傳統(tǒng)煤礦形變的監(jiān)測方法雖然能精確獲取測量點(diǎn)的形變信息[1]，但這種僅依賴有限離散觀測數(shù)據(jù)的方法無法真實(shí)地反映地表形變的動態(tài)進(jìn)程。 近年來，合成孔徑雷達(dá)差分干涉測量技術(shù)(differential interferometric synthetic aperture radar，D-InSAR)以其突出優(yōu)勢發(fā)展迅猛[2-4]，尤其是多種時(shí)間序列InSAR方法的提出[5-7]，有效彌補(bǔ)了D-InSAR方法監(jiān)測精度和可靠性受時(shí)間、空間去相干以及大氣效應(yīng)等因素影響的問題[8-10]。

受自身特殊條件影響，礦區(qū)地表沉降呈現(xiàn)出地形條件復(fù)雜、梯度大、沉降速率大、形變范圍小等特點(diǎn)，使得監(jiān)測難度較大[11-12]。短基線集(small baseline subsets，SBAS)技術(shù)能夠有效解決礦區(qū)地表緩慢形變監(jiān)測問題[13]。在SBAS技術(shù)處理過程中，一般會通過人工選取目標(biāo)點(diǎn)進(jìn)行軌道精煉與沉降反演[14]，這種方式雖然可以有效去除殘留相位與平地效應(yīng)，但若沒有提前對研究區(qū)進(jìn)行考查，獲取地面目標(biāo)點(diǎn)信息就盲目地選擇，反而會引入誤差[15]。由于礦區(qū)開采極易發(fā)生大量級地表沉陷現(xiàn)象，使得地表像元相干性較差，無法輕易準(zhǔn)確獲取穩(wěn)定像元點(diǎn)，因此目標(biāo)點(diǎn)的選取是當(dāng)下礦區(qū)InSAR監(jiān)測中的重點(diǎn)研究問題。

針對以上問題，本文提出一種基于多閾值目標(biāo)提取的SBAS礦區(qū)地表監(jiān)測方法。通過設(shè)定離差閾值參數(shù)、區(qū)域窗口閾值參數(shù)與相干性閾值參數(shù)提取最為穩(wěn)定的目標(biāo)點(diǎn)，并將其應(yīng)用于SBAS-InSAR處理中。以蘭州窯街礦區(qū)為研究區(qū)，探究多閾值目標(biāo)提取的時(shí)序InSAR方法在礦區(qū)地表沉降監(jiān)測中的應(yīng)用效果。

1 方法原理

1.1 基于多閾值目標(biāo)點(diǎn)的提取方法

地物目標(biāo)的干涉相位中一般含有大量噪聲相位，通過衡量噪聲貢獻(xiàn)程度即可選取干涉相位較為純凈的地物目標(biāo)。SAR影像是由多個(gè)像元組成，假設(shè)任意像元的高斯噪聲標(biāo)準(zhǔn)差為σn，相應(yīng)的振幅?信息為Rice分布，可得到式(1)。

(1)

式中：A為某一像元；I為Bessel函數(shù)；g為地面目標(biāo)的正實(shí)數(shù)反射能量。

若該分布近似接近于高斯分布，則需要g/σn>4。此時(shí)，振幅離差指數(shù)可以近似為式(2)。

(2)

式中：σ?為地物目標(biāo)所對應(yīng)的振幅標(biāo)準(zhǔn)差；σnI為虛部標(biāo)準(zhǔn)差；m?為地物目標(biāo)所對應(yīng)的時(shí)序振幅均值；D?為時(shí)序振幅離差指數(shù)。

相干系數(shù)數(shù)學(xué)表達(dá)式見式(3)。

(3)

式中：S1、S2為一組干涉對影像；M、N為局部窗口大小；i、j為提取像元的像素坐標(biāo)。

本文提出的方法首先通過設(shè)定合理的離差指數(shù)閾值參數(shù)對地物目標(biāo)進(jìn)行首次篩選，提取干涉相位較為純凈的地物目標(biāo)。其次，由于相干性能直接反應(yīng)出數(shù)據(jù)處理生成的干涉相位質(zhì)量的高低程度，因此設(shè)定相干性閾值參數(shù)對地物目標(biāo)進(jìn)行第二次篩選，消除大部分相干性差的像素點(diǎn)，最終獲得探測結(jié)果。

在計(jì)算相干系數(shù)時(shí)，局部移動窗口大小的設(shè)定是保證相干系數(shù)精度的重要因素之一。一般情況下，窗口越大，可靠性越高，但在一定程度上降低了分辨率，易將附近的不穩(wěn)定點(diǎn)誤判為穩(wěn)定目標(biāo)點(diǎn)，且易遺漏離散分布單獨(dú)存在的個(gè)別目標(biāo)點(diǎn)。因此，為了平衡圖像分辨率和相干系數(shù)估計(jì)精度這兩個(gè)互為影響的因素，在進(jìn)行第二次篩選之前對研究區(qū)域進(jìn)行窗口劃分。本文采用通過設(shè)定局部移動窗口閾值參數(shù)將整個(gè)區(qū)域劃分成若干個(gè)子區(qū)域，再對每個(gè)子區(qū)域的地物目標(biāo)進(jìn)行相干性篩選，提取大于相干性閾值的目標(biāo)點(diǎn)為參考點(diǎn)，最后對所有子區(qū)域進(jìn)行鑲嵌。

1.2 SBAS-InSAR方法原理

假設(shè)SAR影像在tA和tB(tA>tB)兩個(gè)時(shí)刻干涉得到第j幅干涉圖，在方位向和距離向的坐標(biāo)為(x，r)，則在該點(diǎn)的相位可以表示為式(4)。

δφj(x)=φ(tB,x)-φ(tA,x)≈

(4)

式中：λ為雷達(dá)波長；d(tB，x)和d(tA，x)為雷達(dá)視線方向(LOS)的累計(jì)形變，d(t0,x)=0。

在(x，r)這一像元點(diǎn)，假設(shè)主影像IE、輔影像IS對應(yīng)生成M期干涉圖可表示為式(5)。

(5)

第j(j=1，2，…，M)幅干涉圖的主影像獲取時(shí)間較輔影像晚，即IEj>ISj，則其相位組成的觀測方程見式(6)。

δφj=δφ(tIEj)-δφ(tISj)

(6)

式(6)表示含有N個(gè)未知干涉相位的M個(gè)觀測方程，其矩陣形式見式(7)。

Aφ=δφ

(7)

A是一個(gè)M×N的近似關(guān)聯(lián)矩陣，數(shù)據(jù)中生成的干涉圖可以決定矩陣A。數(shù)據(jù)都在一個(gè)小基線子集里，當(dāng)M接近或者小于N時(shí)，采用奇異值分解的方法可以求得該方程最小范數(shù)意義上的最小二乘解，最終獲取時(shí)序形變量。

2 研究區(qū)概況與數(shù)據(jù)

2.1 研究區(qū)概況

窯街煤礦地處甘肅青海兩省交界處，經(jīng)度范圍為102°50′～102°56′，緯度范圍為36°20′～36°26′，行政區(qū)劃屬蘭州市紅古區(qū)，該區(qū)域地處黃土高原西部，海拔1 800～2 400 m，地形較復(fù)雜[15]。

該礦區(qū)現(xiàn)有窯街三礦、金河煤礦、海石灣煤礦等三個(gè)開采礦區(qū)。由于多年來不斷的地下開采，礦區(qū)內(nèi)地表沉陷、滑坡等地質(zhì)災(zāi)害頻發(fā)，嚴(yán)重影響礦區(qū)內(nèi)各類建設(shè)設(shè)施以及采礦生產(chǎn)安全[16]，研究區(qū)信息如圖1所示。

圖1 研究區(qū)影像圖Fig.1 Image map of the study area

2.2 實(shí)驗(yàn)數(shù)據(jù)準(zhǔn)備

本研究選用具有較短時(shí)間及空間基線的哨兵一號衛(wèi)星(Sentinel-1 A)干涉寬幅模式(IW)下的單視復(fù)數(shù)數(shù)據(jù)，該衛(wèi)星搭載5.6 cm的C波段合成孔徑雷達(dá)，重訪周期為12 d，影像分辨率為5 m×20 m。本文選取的實(shí)驗(yàn)影像時(shí)間跨度為2018年1月—2020年4月，影像雷達(dá)入射角為39.06°，方位向和距離向分辨率為13.91 m與3.70 m。此外，為提高影像配準(zhǔn)的精度，加入精密軌道星歷數(shù)據(jù)降低軌道誤差。參考DEM數(shù)據(jù)采用90 m分辨率SRTM-DEM數(shù)據(jù)。

3 實(shí)驗(yàn)數(shù)據(jù)處理

3.1 多閾值目標(biāo)點(diǎn)提取

傳統(tǒng)SBAS-InSAR方法采用人工選取目標(biāo)點(diǎn)，但在未掌握地面目標(biāo)點(diǎn)信息的情況下進(jìn)行選擇很有可能會增大誤差。因此，本文提出一種基于多閾值的目標(biāo)點(diǎn)提取方法。對原始單視復(fù)數(shù)數(shù)據(jù)進(jìn)行多視配準(zhǔn)之后，得到研究區(qū)強(qiáng)度信息圖，采用本文提出的方法對所有影像數(shù)據(jù)進(jìn)行分析，提取地面目標(biāo)點(diǎn)。首先，為了最大程度提取最優(yōu)點(diǎn)作目標(biāo)點(diǎn)，設(shè)定離差閾值參數(shù)為3.2，進(jìn)行第一次篩選，結(jié)果如圖2(a)所示，實(shí)驗(yàn)共篩選出目標(biāo)點(diǎn)29 453個(gè)。其次，由于研究區(qū)窗口較大，計(jì)算像元相干性時(shí)，在較高相干性區(qū)域易產(chǎn)生簇成團(tuán)的控制點(diǎn)，在較低相干性區(qū)域易將附近的不穩(wěn)定點(diǎn)誤判為控制點(diǎn)，且易遺漏離散分布單獨(dú)存在的個(gè)別控制點(diǎn)。因此，為了平衡高圖像分辨率與高相干系數(shù)精度，設(shè)定局部窗口閾值參數(shù)將整個(gè)區(qū)域分解成若干個(gè)子區(qū)域，再進(jìn)行相干性計(jì)算，經(jīng)過反復(fù)實(shí)驗(yàn)，設(shè)定區(qū)域窗口閾值為16 sqkm，將整個(gè)區(qū)域劃分為132個(gè)子區(qū)間，在每個(gè)子區(qū)間內(nèi)進(jìn)行篩選。最后，設(shè)定相干性閾值參數(shù)為0.9，在第一次篩選結(jié)果的基礎(chǔ)上進(jìn)行第二次篩選，篩選后的地面控制點(diǎn)如圖2(b)所示，實(shí)驗(yàn)最終總共篩選地面控制點(diǎn)103個(gè)，將最終篩選結(jié)果作為后續(xù)軌道精煉的目標(biāo)點(diǎn)用以進(jìn)行地面沉降反演。

圖2 篩選的目標(biāo)點(diǎn)Fig.2 Target point for screening

3.2 改進(jìn)SBAS-InSAR處理

對研究區(qū)數(shù)據(jù)進(jìn)行SBAS-InSAR處理，空間基線閾值設(shè)定臨界極限的5%，時(shí)間基線閾值設(shè)定為180 d，得到的時(shí)間和空間基線連接圖如圖3和圖4所示。

圖3 空間基線連接圖Fig.3 Diagram of spatial baseline connection

圖4 時(shí)間基線連接圖Fig.4 Diagram of time baseline connection

在干涉處理結(jié)束后，檢查所得到的所有干涉像對，去除相干性差的像對，避免形變提取時(shí)產(chǎn)生誤差。而后將基于多閾值方法提取的目標(biāo)點(diǎn)代替人工選點(diǎn)方法，用作SBAS-InSAR軌道精煉與重去平處理。最后，利用奇異值分解法估算形變速率與殘余地形，并通過空間域低通濾波算法與時(shí)間域高通濾波算法去除大氣相位與殘余地形相位，獲取時(shí)序形變速率。

4 結(jié)果分析

4.1 沉降速率分析

經(jīng)過上述處理后，得到研究區(qū)2018年1月—2020年4月期間的雷達(dá)視線方向(LOS)年平均速率圖。由式(8)可以得到研究區(qū)年平均垂直沉降速率情況，如圖5所示。

Δh=Δr/cosθ

(8)

式中：θ為雷達(dá)入射角；Δr為雷達(dá)視線方向沉降量；Δh為垂直方向沉降量。

圖5 研究區(qū)年平均垂直沉降速率圖Fig.5 Annual average vertical settlement rate map

由圖5可知，研究區(qū)內(nèi)大部分地區(qū)地表沉降較為穩(wěn)定，年平均沉降速率主要集中在-24～13 mm/a。A、B、C三處區(qū)域存在明顯地表沉降，具體表現(xiàn)為由邊緣區(qū)域向中心逐漸遞增的沉降漏斗狀，沉降速率為-25～-156 mm/a沉陷和14～87 mm/a抬升，這與文獻(xiàn)[16]得出的結(jié)果基本一致。結(jié)合光學(xué)影像資料對比，圖5中三處沉降漏斗與實(shí)際該煤礦開采工作面基本吻合。圖6為沿HD方向剖面?zhèn)纫晥D，結(jié)合礦區(qū)地下開采資料進(jìn)行對比分析可知，礦區(qū)地下開采是導(dǎo)致三個(gè)沉降區(qū)發(fā)生沉陷的主要原因。

圖6 HD方向剖面?zhèn)纫晥DFig.6 Side view of section in HD direction

4.2 時(shí)序沉降分析

為分析礦區(qū)地表沉降時(shí)序演變規(guī)律，對整個(gè)研究區(qū)進(jìn)行時(shí)間序列分析計(jì)算后，獲得研究期內(nèi)各個(gè)時(shí)間點(diǎn)相對于第一期影像(2018年1月2日)的沉降累計(jì)量，取其中較為重要的11個(gè)時(shí)間節(jié)點(diǎn)沉降累計(jì)圖(圖7)。由圖7可知，2018年5月之前研究區(qū)沒有明顯沉降，沉降量集中在-5～-15 mm之間。從2018年5月開始，研究區(qū)出現(xiàn)較為明顯沉降，并且隨著時(shí)間推移，沉降值和沉降范圍逐漸增大，到2020年4月沉降值達(dá)到-376 mm。其中，2018年9月、2019年1月、2020年4月沉降程度相對較大，沉降值較前分別有58 mm、54 mm和45 mm的增長。

圖7 研究區(qū)時(shí)序沉降量Fig.7 Time series subsidence in the study area

為了進(jìn)一步分析開采沉陷區(qū)沉降發(fā)育過程，在研究區(qū)內(nèi)選取三處(A區(qū)、B區(qū)、C區(qū))沉降最為明顯的區(qū)域，對其做時(shí)序沉降趨勢分析，分析結(jié)果如圖8所示。與研究區(qū)煤礦地下開采資料對比分析，三處區(qū)域分別位于窯街煤礦區(qū)域中的窯街三礦、金河煤礦、海石灣煤礦三個(gè)開采礦區(qū)。三處礦區(qū)在整個(gè)研究時(shí)段內(nèi)大體呈現(xiàn)線性下沉趨勢，其中A區(qū)和B區(qū)整體下沉速率較為穩(wěn)定，C區(qū)呈階段性變化趨勢，在2018年3月之前下沉趨勢較為嚴(yán)重，沉陷值可達(dá)-78 mm，2018年4月之后下沉趨勢較為穩(wěn)定。

圖8 時(shí)序沉降趨勢圖Fig.8 Time series subsidence trend

4.3 監(jiān)測結(jié)果可靠性分析

為驗(yàn)證多閾值目標(biāo)提取的SBAS方法的有效性，使用相同的原始數(shù)據(jù)，利用傳統(tǒng)SBAS-InSAR方法獲取A區(qū)、B區(qū)、C區(qū)的時(shí)序地表沉降值，將兩方法監(jiān)測結(jié)果進(jìn)行對比分析，對比結(jié)果見圖9。

圖9 兩種方法監(jiān)測結(jié)果對比Fig.9 Comparison of monitoring results of two methods

由圖9可知，A區(qū)兩種方法監(jiān)測結(jié)果中差值最大為9.30 mm，最小為0.12 mm，差值的平均值為0.76 mm。對24期數(shù)據(jù)的差值離散程度進(jìn)行計(jì)算，得到標(biāo)準(zhǔn)差范圍在0.06～4.66 mm之間，標(biāo)準(zhǔn)差平均值為2.06 mm。B區(qū)兩種方法監(jiān)測結(jié)果中差值最大為18.80 mm，最小為0.39 mm，差值的平均值為6.80 mm。對24期數(shù)據(jù)的差值離散程度進(jìn)行計(jì)算，得到標(biāo)準(zhǔn)差范圍在0～9.40 mm之間，標(biāo)準(zhǔn)差平均值為3.96 mm。C區(qū)兩種方法監(jiān)測結(jié)果中差值最大為42.20 mm，最小為0.80 mm，差值的平均值為28.70 mm。對24期數(shù)據(jù)的差值離散程度進(jìn)行計(jì)算，得到標(biāo)準(zhǔn)差范圍在0.40～21.10 mm之間，標(biāo)準(zhǔn)差平均值為12.30 mm。

綜上所述，基于多閾值目標(biāo)提取的SBAS方法與傳統(tǒng)SBAS-InSAR得到的監(jiān)測結(jié)果大致相吻合，基于以往SBAS-InSAR方法的監(jiān)測精度，可認(rèn)為本文提出方法得出的結(jié)果是可靠的。

5 結(jié) 論

由于礦區(qū)地表環(huán)境復(fù)雜，很難有效提取穩(wěn)定的目標(biāo)點(diǎn)用以SBAS處理。因此，本文提出了一種基于多閾值目標(biāo)提取的時(shí)序InSAR方法，利用此方法反演出研究區(qū)的年平均沉降速率和時(shí)序累計(jì)沉降值等時(shí)間序列沉降信息，并與傳統(tǒng)SBAS方法監(jiān)測結(jié)果進(jìn)行實(shí)例對比。研究得出結(jié)論如下所述。

1) 兩種方法監(jiān)測結(jié)果對比分析，三處礦區(qū)沉降差值平均值分別為0.76 mm、6.80 mm、28.70 mm，標(biāo)準(zhǔn)差平均值分別為2.06 mm、3.96 mm、12.30 mm，兩種方法監(jiān)測結(jié)果較為吻合，證明基于多閾值目標(biāo)提取的SBAS方法與傳統(tǒng)SBAS-InSAR方法在礦區(qū)地表沉降提取精度均可達(dá)到亞厘米級。

2) 傳統(tǒng)SBAS-InSAR方法為保證提取到更多地表相位較為穩(wěn)定的目標(biāo)點(diǎn)，在選擇研究區(qū)方面表現(xiàn)較為苛刻，這使得其在實(shí)際應(yīng)用過程中存在許多局限性。而本文提出方法可在保證精度的同時(shí)，克服傳統(tǒng)SBAS-InSAR方法選取目標(biāo)點(diǎn)在實(shí)際應(yīng)用中存在的局限性，尤其在礦區(qū)發(fā)生大量級地表沉陷現(xiàn)象，導(dǎo)致地表像元較差的情況下，本文提出方法具有更好的適用性，可為后續(xù)礦區(qū)沉降監(jiān)測分析提供新的監(jiān)測手段。

3) 監(jiān)測結(jié)果表明，研究區(qū)范圍內(nèi)有三處開采沉陷盆地，分別位于三礦(A區(qū))、金河煤礦(B區(qū))、海石灣煤礦(C區(qū))。最大沉陷值為-376 mm，最大年平均沉降速率為-156 mm/a，且沉降下沉趨勢在逐步增加，因此需對該礦區(qū)進(jìn)行長時(shí)間序列的監(jiān)測，用以提供災(zāi)害預(yù)警。