2014-2018年中巴經(jīng)濟(jì)走廊地表變形數(shù)據(jù)集

白艷萍，康建芳，李萌，韓守富，趙寶強(qiáng)，艾鳴浩，馬金輝＊

1. 蘭州大學(xué)資源環(huán)境學(xué)院，蘭州 730000

2. 中國(guó)科學(xué)院西北生態(tài)環(huán)境資源研究院，蘭州 730000

數(shù)據(jù)集（庫(kù)）基本信息簡(jiǎn)介表

引 言

中巴經(jīng)濟(jì)走廊區(qū)域地勢(shì)北高南低，區(qū)域構(gòu)造活躍、地勢(shì)險(xiǎn)峻、冰川活躍，地質(zhì)、水文等工程地質(zhì)條件極為惡劣，公路沿線雪崩、泥石流、坍方、滑坡、巖崩等地質(zhì)災(zāi)害頻發(fā)，嚴(yán)重影響公路的安全運(yùn)營(yíng)[1]。2005年克什米爾7.6級(jí)地震導(dǎo)致公路洪扎河（Huza River）以北段的巖體結(jié)構(gòu)嚴(yán)重破壞，這種歷時(shí)短、規(guī)模大的災(zāi)害已導(dǎo)致中巴公路的路面與路基嚴(yán)重?fù)p毀[2]。類似地質(zhì)災(zāi)害問題在中巴經(jīng)濟(jì)走廊區(qū)域普遍存在，而且隨著未來(lái)工程建設(shè)的進(jìn)行和氣候環(huán)境的變化，地質(zhì)災(zāi)害有可能進(jìn)一步頻發(fā)。

地質(zhì)災(zāi)害識(shí)別、活動(dòng)性調(diào)查與分析以及危險(xiǎn)性評(píng)價(jià)是走廊內(nèi)公路、鐵路、通信、油氣管道等重大工程建設(shè)、安全運(yùn)行與科學(xué)減災(zāi)的基礎(chǔ)。基于此，學(xué)術(shù)界和工程界采用多種手段開展了該區(qū)域地質(zhì)災(zāi)害研究。但由于該區(qū)域自然環(huán)境復(fù)雜、野外實(shí)際勘測(cè)難度很大，基于先進(jìn)的衛(wèi)星遙感技術(shù)和時(shí)序InSAR技術(shù)進(jìn)行地表變形的反演和地質(zhì)災(zāi)害的早期識(shí)別、監(jiān)測(cè)預(yù)警，已成為值得重視的研究手段。星載合成孔徑雷達(dá)干涉測(cè)量（Interferometric synthetic aperture radar，InSAR）因其高精度、高分辨率、全天候等優(yōu)點(diǎn)已迅速成為常用的大地測(cè)量技術(shù)之一[3-4]，旨在通過(guò)對(duì)同一區(qū)域兩次或多次過(guò)境的合成孔徑雷達(dá)（Synthetic Aperture Radar，SAR）影像的復(fù)共軛相乘，來(lái)提取地物目標(biāo)的地形或者形變信息[5]。如易靖松通過(guò)時(shí)序 InSAR技術(shù)進(jìn)行地表形變監(jiān)測(cè)，分析潛在滑坡的時(shí)空分布特征[6]。朱建軍提出 SAR數(shù)據(jù)的全球無(wú)縫覆蓋和高密度重返都會(huì)為潛在災(zāi)害識(shí)別和地質(zhì)監(jiān)測(cè)提供重要的數(shù)據(jù)源[5]。Hilleyge利用永久散射體星載合成孔徑雷達(dá)干涉測(cè)量（Persistent scatterer InSAR，PSInSAR）技術(shù)對(duì)美國(guó)Berkley地區(qū)的滑坡展開監(jiān)測(cè)，結(jié)果較好地反應(yīng)了滑坡形變的歷史演變進(jìn)程[7]。均可以說(shuō)明利用時(shí)序InSAR變形監(jiān)測(cè)技術(shù)，識(shí)別并監(jiān)測(cè)潛在地質(zhì)災(zāi)害的可能性。潛在災(zāi)害的識(shí)別和監(jiān)測(cè)是地質(zhì)災(zāi)害成災(zāi)機(jī)制研究風(fēng)險(xiǎn)評(píng)估和工程防治的基礎(chǔ)，也是科學(xué)預(yù)測(cè)自然災(zāi)害發(fā)展、演化趨勢(shì)，合理評(píng)估潛在自然災(zāi)害影響的關(guān)鍵科學(xué)問題。

本數(shù)據(jù)集為中巴經(jīng)濟(jì)走廊全部區(qū)域時(shí)序地表形變數(shù)據(jù)，通過(guò)PSInSAR技術(shù)對(duì)2014-2018年間覆蓋中巴經(jīng)濟(jì)走廊全部區(qū)域的Sentinel 1A數(shù)據(jù)進(jìn)行地表形變監(jiān)測(cè)而獲取，為中巴經(jīng)濟(jì)走廊的地質(zhì)災(zāi)害識(shí)別和經(jīng)濟(jì)建設(shè)提供參考數(shù)據(jù)。

1 數(shù)據(jù)采集和處理方法

1.1 原始數(shù)據(jù)說(shuō)明

1.1.1 DEM數(shù)據(jù)

DEM 數(shù)據(jù)是進(jìn)行地形因子計(jì)算的基礎(chǔ)，本數(shù)據(jù)集選擇 SRTM 90 m分辨率的 DEM 數(shù)據(jù)（https://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/），基礎(chǔ)影像共19景，將其進(jìn)行拼接、裁剪和空值填補(bǔ)后建立了覆蓋走廊范圍的完整 DEM，高程范圍為 986～7770 m，空間參考橢球體為WGS_1984。

1.1.2 InSAR數(shù)據(jù)

本研究根據(jù)SAR數(shù)據(jù)的分辨率、波長(zhǎng)、觀測(cè)角度、性價(jià)比和可獲得性等，針對(duì)研究區(qū)地形和變形特征選用歐空局合成孔徑雷達(dá)（C波段）Sentinel 1A數(shù)據(jù)（https://scihub.copernicus.eu/dhus/#/home。）。所涉及軌道為相對(duì)軌道107（下降軌道）、27（上升軌道）。覆蓋范圍如圖1所示。時(shí)間范圍從2014年10月到2018年6月。其中107軌道共400景、覆蓋74天，27軌道共340景、覆蓋85天，數(shù)據(jù)基本特征如表1所示。空間覆蓋范圍如圖1所示，時(shí)間覆蓋范圍如圖2所示。107軌道和27軌道影像的采樣時(shí)間分布如圖2。雖然107軌道缺少20161014、20161026、20161107、20161119、20161201、20161213、20161225等日期的數(shù)據(jù)，但是采樣時(shí)間分布相對(duì)均勻，可以滿足圖像高精度配準(zhǔn)要求。

圖1 數(shù)據(jù)集覆蓋范圍

圖2 影像數(shù)據(jù)時(shí)間分布圖

表1 Sentinel 1A 數(shù)據(jù)基本信息

1.2 數(shù)據(jù)集處理方法及過(guò)程

1.2.1 采用原理和技術(shù)路線

本文采用PSInSAR技術(shù)獲得地表形變量，其基本原理是利用多景同一地區(qū)的SAR影像，通過(guò)統(tǒng)計(jì)分析時(shí)間序列上幅度和相位信息的穩(wěn)定性，探測(cè)不受時(shí)間、空間基線去相關(guān)影響的穩(wěn)定點(diǎn)目標(biāo)[8]。PSInSAR技術(shù)通過(guò)對(duì)主輔影像進(jìn)行干涉處理，將獲取的強(qiáng)度與相位信息進(jìn)行統(tǒng)計(jì)分析，選取在時(shí)間上散射特性相對(duì)穩(wěn)定、回波信號(hào)較強(qiáng)的永久散射體（Persistent Scatterer，PS）目標(biāo)點(diǎn)作為觀測(cè)對(duì)象。這些目標(biāo)點(diǎn)可能是人工建筑物、裸露的巖石、人工布設(shè)的角反射器等。由于它們?cè)跁r(shí)間序列SAR影像中幾乎不受斑點(diǎn)噪聲的影響，經(jīng)過(guò)長(zhǎng)時(shí)間間隔仍然保持穩(wěn)定的散射特性，所以被稱作PS。構(gòu)建PS點(diǎn)網(wǎng)絡(luò)，網(wǎng)絡(luò)中任意相鄰PS點(diǎn)間經(jīng)過(guò)差分處理的方式分離出形變信息，進(jìn)而求出相鄰點(diǎn)間的高程差值和形變速度差值[9]。以高程差和形變速度差作為原始數(shù)據(jù)，利用相位解纏算法求解出所有PS點(diǎn)的形變、高程參數(shù)，再通過(guò)最小二乘法等求得研究區(qū)內(nèi)PS點(diǎn)的形變量和形變速度[10]。

本文采用SARProZ軟件——SAR影像處理系統(tǒng)，主要面向雷達(dá)影像數(shù)據(jù)PSInSAR處理技術(shù)及相關(guān)應(yīng)用，基于并行多線程處理技術(shù)，能夠快速處理分析海量雷達(dá)影像數(shù)據(jù)，可識(shí)別區(qū)分建筑物和地面上的PS點(diǎn)，區(qū)分建筑物形變和地表形變；精確反演影像采集時(shí)間的大氣相位圖，剔除大氣相位殘值；對(duì)于城市區(qū)域或強(qiáng)反射硬目標(biāo)形變速率監(jiān)測(cè)精度可達(dá)3～5 mm。其數(shù)據(jù)處理流程如圖3。

1.2.2 中巴經(jīng)濟(jì)走廊PSInSAR的數(shù)據(jù)處理過(guò)程

（1）中巴經(jīng)濟(jì)走廊計(jì)算范圍設(shè)計(jì)

計(jì)算范圍選取107軌道和27軌道覆蓋的中國(guó)境內(nèi)的經(jīng)濟(jì)走廊區(qū)域。具體計(jì)算時(shí)，采用分幅計(jì)算的方式，各分幅計(jì)算范圍如圖4所示。本文僅展示境內(nèi)重點(diǎn)區(qū)域的計(jì)算結(jié)果。

境內(nèi)重點(diǎn)區(qū)域的覆蓋范圍為圖4中紫色區(qū)域所示，涵蓋了慕士塔格山和塔什庫(kù)爾干塔吉克自治縣兩個(gè)典型沉降區(qū)，其四角點(diǎn)經(jīng)緯度分別為：（38°42′12″N，75°03′41″E）、（38°34′16″N，75°48′54″E）、（37°26′35″N，75°29′34″E）、（37°31′22″N，74°46′01″E）。面積約為 8000 km2。

圖3 PSInSAR數(shù)據(jù)處理流程圖

圖4 樣本數(shù)據(jù)覆蓋范圍

（2）InSAR主影像選擇和影像連接拓?fù)鋱D構(gòu)建

針對(duì)研究區(qū)Sentinel 1A數(shù)據(jù)集，獲取所有影像對(duì)的時(shí)間、空間基線，將時(shí)間、空間基線之和為最小的影像作為主影像，構(gòu)建時(shí)間和空間影像連接拓?fù)鋱D（圖5），其中主影像日期為2016年12月31日。

圖5 PSInSAR影像連接拓?fù)鋱D

（3）影像配準(zhǔn)

配準(zhǔn)過(guò)程中先采用精確軌道參數(shù)進(jìn)行粗配準(zhǔn)，為精配準(zhǔn)提供初始偏移量，以加快精配準(zhǔn)速度。然后參考90 m分辨率的SRTM數(shù)據(jù)，選擇整景影像上均勻分布的同名像素點(diǎn)建立多項(xiàng)式方程進(jìn)行精配準(zhǔn)，配準(zhǔn)精度小于1/8個(gè)像元。

（4）干涉相位計(jì)算

采用Goodstein mode濾波器，濾波窗口為15×15。對(duì)已經(jīng)過(guò)濾波的主輔影像像元對(duì)生成干涉相位值，逐像元計(jì)算干涉相位，且采用DEM去除地形相位，生成時(shí)間序列干涉圖集。干涉圖用于評(píng)價(jià)InSAR的計(jì)算結(jié)果，干涉環(huán)越明顯，形變效果越好。根據(jù)圖6所示干涉圖（共85幅干涉圖，本文選擇部分干涉圖進(jìn)行展示），可以看出殘留地形相位不明顯，且存在明顯干涉環(huán)，表明各時(shí)期影像與主影像相干性較好，干涉結(jié)果可以用于后繼PSInSAR計(jì)算。

圖6 干涉處理結(jié)果

（5）干涉相位解纏

由于干涉圖整體相干性較高時(shí)，采用枝切法進(jìn)行相位解纏，解纏效果良好。

（6）候選PS點(diǎn)選擇、空間網(wǎng)絡(luò)構(gòu)建、大氣效應(yīng)剔除

在進(jìn)行候選PS點(diǎn)計(jì)算時(shí)，采用Amplitude Stable. index進(jìn)行PS點(diǎn)選擇，閾值設(shè)置為0.65，線性趨勢(shì)估值設(shè)置為-20和20，高程的估值參數(shù)設(shè)置為-50和50，多項(xiàng)式項(xiàng)數(shù)設(shè)置為1，計(jì)算過(guò)程中將相干性（Coherence）作為權(quán)重，以此估算大氣效應(yīng)，并反演得到地表形變速度和各時(shí)期形變量，最終PS點(diǎn)的相干性大于0.7，符合國(guó)際PS點(diǎn)選擇的慣例。

2 數(shù)據(jù)樣本描述

本數(shù)據(jù)集采用PSInSAR技術(shù)獲取中巴經(jīng)濟(jì)走廊地表形變數(shù)據(jù)，本文選擇兩個(gè)典型的境內(nèi)重點(diǎn)區(qū)域?qū)?shù)據(jù)集樣本進(jìn)行描述。第一個(gè)典型重點(diǎn)區(qū)域位于塔什庫(kù)爾干塔吉克自治縣，以東經(jīng) 75°17′，北緯37°39′為區(qū)域中心，沿軌道方向長(zhǎng)為98公里、寬42公里的長(zhǎng)方形區(qū)域。圖7所示為27軌道數(shù)據(jù)計(jì)算結(jié)果，圖8所示為107軌道數(shù)據(jù)計(jì)算結(jié)果。對(duì)比兩條軌道數(shù)據(jù)的計(jì)算結(jié)果可以看出，庫(kù)孜滾村的東南方向2.2公里的位置均發(fā)生明顯的下沉，下沉位置完全吻合。其余地方?jīng)]有明顯變形。

圖7 PSInSAR計(jì)算結(jié)果樣本數(shù)據(jù)（27軌道）

圖8 PSInSAR計(jì)算結(jié)果樣本數(shù)據(jù)（107軌道）

第二個(gè)典型重點(diǎn)區(qū)域位于克孜勒蘇柯爾克孜自治州阿克陶縣，有著世界冰山之父之名的慕士塔格山，以東經(jīng)75°08′，北緯38°13′為區(qū)域中心。本文選用27軌道的計(jì)算結(jié)果，圖9所示的慕士塔格山的冰川下緣，沿著冰川觸水徑流方向有明顯的下沉現(xiàn)象，其在山腳下西南方向也有局部明顯沉降。

3 數(shù)據(jù)質(zhì)量控制

本數(shù)據(jù)集基礎(chǔ)數(shù)據(jù)選用歐空局合成孔徑雷達(dá)（C波段）Sentinel 1A數(shù)據(jù)，它具有良好的基線控制和短暫的重返周期，能夠在短暫時(shí)間積累大量影像數(shù)據(jù)，有利于InSAR時(shí)序分析，并且提升雷達(dá)圖像參數(shù)反演的精度。在PSInSAR數(shù)據(jù)處理過(guò)程中，通過(guò)影像配準(zhǔn)、干涉成像、相位解纏、高程解算和地理編碼等關(guān)鍵步驟，有效地去除地平、高程相位以及大氣效應(yīng)，獲得了理論精度毫米級(jí)的形變信息。并且文中采用了107軌道和27軌道兩個(gè)軌道，分別對(duì)研究區(qū)進(jìn)行計(jì)算，不同軌道計(jì)算結(jié)果吻合較好。

4 數(shù)據(jù)使用方法和建議

使用數(shù)據(jù)時(shí)，可利用ArcGIS地理信息系統(tǒng)軟件對(duì)shape文件進(jìn)行查詢、編輯等操作，其shape中的屬性表可作為分析與統(tǒng)計(jì)地表形變的基本數(shù)據(jù)。用Google Earth可對(duì)KML文件進(jìn)行查詢、分析等操作。本數(shù)據(jù)集可靠性高，代表性較強(qiáng)，可以結(jié)合地質(zhì)巖性結(jié)構(gòu)數(shù)據(jù)，用于中巴經(jīng)濟(jì)走廊地質(zhì)變化研究，滿足用戶對(duì)中巴經(jīng)濟(jì)走廊自然災(zāi)害科學(xué)數(shù)據(jù)的需求。