吳 明，嚴(yán) 珊，王金碩，楊 敏，張晨蕾，王玥明

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

0 引言

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

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

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

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

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

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

2 沉降模型

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

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

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

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

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

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

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

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

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

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

處理結(jié)果：

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

圖4 府谷塌陷地震雷達(dá)干涉序列圖

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

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

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

2017年9月5日01:14:22，據(jù)中國(guó)地震臺(tái)網(wǎng)測(cè)定，陜西省彬州市（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時(shí)間跨度也較短，為大約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)狀干涉，表示沉降進(jìn)入stage B，塌陷發(fā)生之后，迅速進(jìn)入stage D，干涉圖上表現(xiàn)的沉降活動(dòng)微弱。從stage B的干涉環(huán)測(cè)量得到的沉降區(qū)域最大半徑約為0.6 km，規(guī)模相對(duì)較小。

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

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

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

4 討論與分析

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

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

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

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

5 結(jié)論

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

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

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