基于InSAR技術(shù)對陜西礦區(qū)塌陷地震早期預警的研究

吳 明，嚴 珊，王金碩，楊 敏，張晨蕾，王玥明

（1.陜西省地震局,西安710038；2.中國地震局第二監(jiān)測中心,西安710054）

0 引言

由于經(jīng)濟發(fā)展帶來的資源需求，陜西地區(qū)礦產(chǎn)資源開采持續(xù)增長，同時也帶來地面塌陷等伴生災害，根據(jù)陜西地震臺網(wǎng)記錄到的塌陷地震來看，其發(fā)生的區(qū)域基本與礦區(qū)位置吻合。隨著陜北能源基地的建設(shè)和煤炭資源的大量開采，采空區(qū)的塌陷災害越來越嚴重，伴生的地面形變、地下水流失、土地沙化、房屋開裂等問題對當?shù)鼐用竦纳a(chǎn)生活造成了困擾，自2004年以來陜西數(shù)字地震臺網(wǎng)記錄到的塌陷地震事件越來越多[1]。榆林地區(qū)位于鄂爾多斯地塊，構(gòu)造穩(wěn)定，地震活動較弱，自公元前780年至1960年記載的地震共9次，均小于6級，1990年至今未發(fā)生過除塌陷類型以外的地震，未見有地震形成的災害記錄[2]。而從中國地震臺網(wǎng)記錄的2004—2019年地震監(jiān)測數(shù)據(jù)來看，陜西省有記錄的塌陷地震多達324起，震級最大的是2004年10月14日發(fā)生在陜西省榆林市神木縣的ML4.2塌陷地震。因此采空區(qū)的塌陷地震早期預警研究具有重要意義。

塌陷監(jiān)測預警分為地下監(jiān)測和地表監(jiān)測，地下監(jiān)測如應力應變監(jiān)測[3]、土壤含水率監(jiān)測、爆破震動監(jiān)測等方法[4]，通常需要在地下埋設(shè)大量傳感器與設(shè)備，監(jiān)測預警成本高且應用范圍??；地表監(jiān)測包括全球定位系統(tǒng)（GPS）、熱紅外（TI）、合成孔徑雷達（InSAR）、大地水準測量等。傳統(tǒng)的地面形變監(jiān)測手段如大地水準測量、GPS等雖然方法成熟且單點精度較高，但是其仍存在成本高和空間分辨率低等問題。而InSAR技術(shù)因具有覆蓋范圍廣、全天時、全天候等特點被廣泛應用于區(qū)域形變監(jiān)測研究中[5-6]。1996年Carnec等人[7]首次利用D-InSAR技術(shù)監(jiān)測法國GARDANNE附近煤礦開采引起的地表形變，驗證了InSAR技術(shù)監(jiān)測礦區(qū)形變的可行性，2001年Berardino等[8]提出了SBAS-InSAR方法。InSAR用于礦區(qū)地面沉降監(jiān)測已有20多年的發(fā)展歷史，在滑坡監(jiān)測、同震形變監(jiān)測等領(lǐng)域獲得了很多成果[9-13]。隨著InSAR技術(shù)的發(fā)展，Ps-InSAR、SBAS-InSAR等基于時序的InSAR技術(shù)開始應用，在克服形變監(jiān)測中的大氣效應、時空失相干等方面效果更加顯著。為了更科學地研究塌陷地震預警，我們綜合應用InSAR以及沉降規(guī)律模型開展研究。

1 數(shù)據(jù)源及研究區(qū)域

地面沉降是一個時空演化的過程，一些發(fā)展比較快速的沉降可能會演化為塌陷，塌陷大多發(fā)生于沉降的發(fā)展期中，塌陷的發(fā)展也是一個時空漸變的演化過程，在這個過程中地表形變速率有其發(fā)展的演化規(guī)律，那么，在前期速率由慢到快的發(fā)展階段中，如果能夠使用InSAR技術(shù)手段監(jiān)測到塌陷發(fā)生前的速率特征，那么或許可以為塌陷地震的前期預警提供幫助。塌陷地震引發(fā)的地面垂直形變較大，塌陷區(qū)相位梯度過大導致失相干較嚴重，研究最終獲得的信息較少，目前基于InSAR的塌陷地震研究相對較少，對于塌陷地震早期預警的研究也存在困難。通過比較發(fā)現(xiàn)C波段SAR數(shù)據(jù)是用于觀測塌陷過程活躍階段沉降邊緣變形帶的較理想數(shù)據(jù)，用X波段數(shù)據(jù)觀測會因為相位梯度太小導致可觀測變形帶過窄，而使用L波段數(shù)據(jù)則因為波長過長，邊緣形變量不夠，導致無法在干涉圖上形成明顯的環(huán)狀干涉條帶。而如果使用時間基線過長的SAR數(shù)據(jù)做干涉，則會難以判斷沉降所處的階段，所以建議使用時間基線盡量短的數(shù)據(jù)來觀測，時間基線越短，可做出的判斷越準確。基于觀測效果以及時間成本綜合考慮，我們選擇C波段Sentinel-1A衛(wèi)星數(shù)據(jù)進行處理（數(shù)據(jù)源Sentinel-1A來自歐空局網(wǎng)站：https://scihub.esa.int/）。

根據(jù)中國地震臺網(wǎng)數(shù)據(jù)統(tǒng)計，2004—2019年陜西地區(qū)發(fā)生的324起塌陷地震，榆林地區(qū)達293起，占據(jù)了陜西省塌陷地震的絕大多數(shù)，同時考慮到本研究采用的哨兵衛(wèi)星數(shù)據(jù)從2014年開始升空監(jiān)測，我們從2014—2019年塌陷地震中選擇了震級最大的2次地震作為研究對象，分別為2017年9月5日陜西彬州ML3.2塌陷地震（以下簡稱彬州ML3.2塌陷地震）和2018年8月16日陜西府谷ML3.2塌陷地震（以下簡稱府谷ML3.2塌陷地震）（圖1）。ML3.2塌陷地震震中位于陜西省北部礦區(qū)，陜北地區(qū)地層松軟，煤礦埋藏較淺，上覆巖層普遍薄。陜北位于鄂爾多斯地塊內(nèi)部，這個塊體主要以緩慢的升降震蕩運動為主，具有整體性和穩(wěn)定性的特征，除了塊體周緣發(fā)育不同程度的褶皺與斷裂，塊體內(nèi)部完整且穩(wěn)定，構(gòu)造活動較少，沒有發(fā)生過6級以上地震[1]，在該區(qū)域記錄到的地震主要是采礦引發(fā)的塌陷地震。彬州ML3.2塌陷地震震中位于陜西省彬州市與咸陽市的交界處。

圖1 研究區(qū)域示意圖

2 沉降模型

對于采空區(qū)塌陷的理論模型也有很多學者做了相關(guān)研究，20世紀90年代錢鳴高院士首次提出了關(guān)鍵層理論用于揭示采空區(qū)上覆巖層力學機理，該理論認為采空區(qū)上覆巖層中較為堅硬的厚層狀巖層控制整個上覆巖層和地表的移動與變形，這種較厚層狀的堅硬巖層稱為關(guān)鍵層[14]。20世紀50年代波蘭學者Knothe對地下采礦引起的的地表塌陷進行了研究，基于地表下沉速度和該時刻下沉量與最終下沉量的差值成正比的假設(shè)前提，精確建立了下沉量時間序列模型。由于Knothe提出的模型與很多實際的監(jiān)測資料擬合的結(jié)果不符，李遠耀[15]在Knothe模型的基礎(chǔ)上，提出了修正的時間序列模型，在Knothe的時間函數(shù)上加了一個冪指數(shù)k，k為地面下沉系數(shù)，增加參數(shù)后，時間函數(shù)為：

李遠耀的模型模擬的沉降過程大致分加速度正負2個階段。郭增長[16]按照下沉速度對地表的影響程度不同，將地表下沉分為3個階段：①開始階段：下沉量達到10 mm的時刻為地表下沉開始階段，從下沉開始到沉降漏斗下沉速度達到1.67 mm/d為開始階段。②活躍階段：沉降漏斗下沉速度大于1.67 mm/d的階段，該階段的下沉量占總下沉量的85%～90%，地表的損壞主要發(fā)生在這個階段。③衰退階段：沉降漏斗下沉速度降到1.67 mm/d以下，并且半年內(nèi)各點的累計沉降量不超過30 mm。

通過研究分析了許多學者們用各種手段監(jiān)測得到的塌陷沉降時間序列曲線后，我們發(fā)現(xiàn)了一些有趣的規(guī)律，總結(jié)成如圖2的沉降規(guī)律模型，根據(jù)各時期的沉降量，將沉降分為4個階段，即開始階段（stage A）、前活躍階段（stage B）、后活躍階段（stage C）、穩(wěn)定階段（stage D）。圖中2條黑色曲線屬于開始階段（stage A），這個階段整個沉降漏斗沉降速率都低于1.67 mm/d；紫色曲線屬于前活躍階段（stage B），這一階段的特點是，整體下沉速度急劇加大，伴隨著沉降漏斗的邊緣水平擴展，這一時期漏斗邊緣基本成型，stage B的發(fā)展速度很快，持續(xù)時間通常也比較短；紅色曲線屬于后活躍階段（stage C），該階段的特點是，漏斗中央沉降速度依然很大，但是漏斗邊緣擴展的速度急劇下降，漏斗的中央逐漸發(fā)生超臨界沉降，隨著開采的推進，發(fā)展成為盆地，盆地中心的最大沉降值保持不變，沉降中心底部逐漸變平，stage C的持續(xù)時間通常受開采深度、開采方式、上覆巖層厚度及堅硬程度等影響；藍色曲線屬于穩(wěn)定階段（stage D），這個階段整個沉降漏斗速度都低于1.67 mm/d，由于在前一階段發(fā)生了塌陷，應力快速釋放并形成新的平衡，漏斗內(nèi)的活動速率迅速降低。塌陷絕大多數(shù)都發(fā)生在stage C，后活躍階段采空區(qū)上覆關(guān)鍵層全部斷裂，壓力平衡拱達到極限狀態(tài)，此時意外的外力因素極易導致應力失穩(wěn)形成塌陷。陜西的煤礦埋藏深度普遍較淺，上覆黃土和巖層，巖層一般較薄。

圖2 采空區(qū)沉降模型圖

如果我們假設(shè)的沉降模型合理，那么使用重訪周期12 d的Sentinel-1A衛(wèi)星C波段SAR數(shù)據(jù)對各個沉降階段的地面進行干涉處理，濾波和去平之后，可能得到的干涉的結(jié)果如下：開始階段由于整體下沉速率低，在12 d累計最大沉降不超過20 mm，干涉圖上應無明顯現(xiàn)象。前活躍階段，漏斗中央速度迅速增大，超越C波的相位梯度，而漏斗邊緣也在水平和垂直向有所運動，但是相位梯度相對中央平緩許多，在漏斗外邊緣和內(nèi)邊緣之間可能形成一個不超越C波段形變梯度的沉降環(huán)，所以這一階段得到的干涉圖可能會在沉降漏斗外邊緣附近形成若干干涉條紋，而條紋內(nèi)由于沉降速度超越形變梯度閾值，導致干涉條紋發(fā)生混疊。后活躍階段，漏斗中央沉降速度依然很快，但漏斗邊緣速度卻急劇下降，所以中央由于梯度過大無法形成干涉，而邊緣由于速度太低也無法形成干涉條紋。穩(wěn)定階段整體沉降速度都低于1.67 mm/d，也無法形成干涉條紋。也就是說，如果此沉降模型符合實際，那么使用重訪周期12 d的C波段SAR影像來不斷對整個沉降過程的地面做干涉的話，應該會在第2階段前活躍階段觀測到沉降漏斗邊緣的活躍，表現(xiàn)在干涉圖像上應該是在沉降中心的邊緣地帶形成若干環(huán)狀的干涉條紋。陜西的礦區(qū)大部分位于黃土覆蓋區(qū)，原濤[17]通過相似材料模擬實驗和FLAC數(shù)值模擬揭示了黃土區(qū)開采沉陷的基本規(guī)律，認為厚黃土層礦區(qū)地表移動具有發(fā)展快、穩(wěn)定快、活躍期短的特征，且黃土層厚度是影響地表動態(tài)移動規(guī)律的重要因素。因此在黃土覆蓋的薄巖層礦區(qū)第3階段的持續(xù)時間較短，觀測到邊緣的環(huán)狀干涉條帶距離塌陷發(fā)生時間也相應較短。

3 數(shù)據(jù)處理與驗證

在研究中我們選擇了最近幾年2次震級相對較大的塌陷地震事件，2次地震分別發(fā)生在府谷縣三道溝附近和彬州市新民鎮(zhèn)東務村附近，府谷和彬州都擁有豐富的煤礦資源，也長期受到采空區(qū)地面沉降的困擾。使用短周期時序InSAR對2次地震事件進行分析，驗證提出的模型是否符合預期。

3.1 2018年8月16日府谷ML3.2塌陷地震

據(jù)榆林地震臺網(wǎng)測定2018年8月16日13:09分在陜西榆林市府谷縣三道溝鎮(zhèn)（39.16°N，110.73°E）發(fā)生ML3.2地震（地震類型塌陷）（圖3）。我們嘗試用重訪周期12 d的C波段升軌哨兵衛(wèi)星數(shù)據(jù)，對發(fā)生塌陷的區(qū)域做干涉處理，影像從2018年1月1日到2018年8月29日，包含了共21景SAR數(shù)據(jù)及精密軌道數(shù)據(jù)，獲取以12 d為間隔的20對干涉圖像，干涉后的影像進行濾波和去平處理。

圖3 府谷塌陷地震位置示意圖（底圖來源于Google earth）

處理結(jié)果：

從處理后獲得的干涉序列（圖4）上來看，2018年1月1日—7月24日間的17幅干涉圖屬于沉降的stageA階段，2018年7月24日—8月5日的干涉圖出現(xiàn)了非常明顯的環(huán)狀干涉條紋，表示沉降進入stageB，沉降速度明顯加大且沉降范圍迅速擴大。塌陷發(fā)生于2018年8月16日，說明在該區(qū)域stageC時間跨度較短，長度大約為10～20 d，在塌陷發(fā)生后，沉降迅速進入stageD，沉降速率快速下降，在干涉圖上無明顯干涉現(xiàn)象。出現(xiàn)在stageB的干涉環(huán)最大半徑達到2 km。

圖4 府谷塌陷地震雷達干涉序列圖

SAR影像基本參數(shù)見表1，參考DEM為SRTM3 DEM數(shù)據(jù)，地面分辨率90 m。

表1 哨兵2A影像基本參數(shù)列表

3.2 2017年9月5日彬州市ML3.2級塌陷地震

2017年9月5日01:14:22，據(jù)中國地震臺網(wǎng)測定，陜西省彬州市（35.09°N，108.18°E）發(fā)生 ML3.2 地震（塌陷），震源深度0 km（圖5）。選擇從2017年4月2日—10月31日的15幅哨兵升軌數(shù)據(jù)，干涉處理后獲得14張間隔周期為12 d的干涉影像，濾波和去平處理后如圖5所示。塌陷發(fā)生于2017年9月5日，說明在彬州市的這次塌陷stageC時間跨度也較短，為大約5～15 d。

圖5 彬州塌陷地震位置示意圖（底圖來源于Google earth）

處理結(jié)果：

從干涉的結(jié)果（圖6）來看，2017年4月22日—9月25日之間的8幅干涉影像都屬于stage A，2017年8月20日—9月1日的干涉影像上出現(xiàn)明顯環(huán)狀干涉，表示沉降進入stage B，塌陷發(fā)生之后，迅速進入stage D，干涉圖上表現(xiàn)的沉降活動微弱。從stage B的干涉環(huán)測量得到的沉降區(qū)域最大半徑約為0.6 km，規(guī)模相對較小。

圖6 彬州塌陷地震干涉序列圖

SAR影像基本參數(shù)見表2，參考DEM為SRTM3 DEM數(shù)據(jù)，地面分辨率90 m。

表2 哨兵2A影像基本參數(shù)列表

4 討論與分析

從2次塌陷地震的干涉時間序列來看，發(fā)生于黃土覆蓋區(qū)的礦區(qū)塌陷地震，短周期的C波段雷達數(shù)據(jù)干涉均可以在塌陷地震發(fā)生前明顯地識別出前活躍階段（stage B）（圖7）。2次塌陷地震從stage B到地震發(fā)生之間的時間間隔較短，且地震發(fā)生后沉降活動迅速減弱，印證了原濤[17]關(guān)于厚黃土層礦區(qū)地表移動形變具有發(fā)展快、穩(wěn)定快、活躍期短的觀點。影響stage C時間長短的因素有很多，主要可以分為自然地質(zhì)因素和采礦技術(shù)因素[16]，自然地質(zhì)因素如煤層厚度、煤層傾角、埋藏深度、松散層厚度、斷層密度、地層不整合、巖石堅硬系數(shù)等，采礦技術(shù)因素如采空區(qū)位置、開采方法、工作面推進速度、頂板控制方法等。如果要準確地估算可能的塌陷地震與stage B的時間間隔，需要綜合前面的自然地質(zhì)因素與采礦技術(shù)因素來綜合考量。通過短周期的時序InSAR干涉測量判斷潛在塌陷地震區(qū)位置、所處階段，結(jié)合現(xiàn)場考察或資料搜集等獲得塌陷區(qū)的自然地質(zhì)因素以及采礦技術(shù)因素，可以為塌陷地震的早期預警提供有力的支撐。

圖7 采空區(qū)在階段B的三維剖面示意圖

2017年9月5日彬州ML3.2級塌陷地震在SAR的干涉影像上不如2018年8月16日府谷ML3.2級塌陷地震表現(xiàn)得明顯，主要是塌陷范圍小導致，從InSAR觀測結(jié)果上分析，彬州塌陷的最大直徑約600 m，而府谷塌陷的最大直徑為2 000 m。如果使用波長較短的X波段SAR數(shù)據(jù)來研究彬州塌陷地震效果會更加明顯。

需要注意的是，即使在短周期的干涉圖像中發(fā)現(xiàn)明顯的環(huán)狀干涉區(qū)域，該區(qū)域在未來一段時間是否會發(fā)生塌陷地震以及什么時候發(fā)生地震，這仍然需要結(jié)合區(qū)域?qū)嶋H的地質(zhì)地貌條件、采空區(qū)深度等多種因素來綜合判斷?？紤]到哨兵雷達影像的分辨率不高，如果塌陷的范圍過小，會影響到整體的監(jiān)測效果。小范圍的塌陷使用波長更短分辨率更高的雷達數(shù)據(jù)或許可以獲得更佳的監(jiān)測效果。

5 結(jié)論

結(jié)合前人對地表下沉速率與周期的研究以及大量現(xiàn)場監(jiān)測得到的沉降數(shù)據(jù)，總結(jié)沉降發(fā)育的時空規(guī)律經(jīng)驗模型，配合C波段SAR數(shù)據(jù)的特性獲得的短周期干涉序列，可以大范圍、低成本地精確篩選出處于stage B的潛在塌陷區(qū)域，結(jié)合野外實際調(diào)查以及資料收集獲得的采空區(qū)深度、上覆巖層厚度及堅硬系數(shù)、開采方法等信息，可以綜合判斷塌陷是否會發(fā)生以及可能發(fā)生在什么時間，從而為塌陷地震的早期預警提供支撐。

通過在黃土覆蓋礦區(qū)開展塌陷地震預警研究工作，可為未來在更廣泛的區(qū)域監(jiān)測塌陷地震、做塌陷風險評估、提前預警等開展更多研究提供經(jīng)驗。

致謝感謝袁志祥研究員在論文寫作過程中的悉心指導，給本文寫作提供了建設(shè)性的意見，在此表示感謝。