SARscape的DInSAR全流程

写在前面：

       鉴于隋大兴唐长安城遗址在面积上的尺度较大、现存遗产单体要素在空间上的跨度较为分散，与当今西安市主城区高度重合的实际情况。结合西安市地质活动长期发育，严重影响文化遗产安全及遗址本体保护的现实问题。兼顾开展科研攻关先遣工作：聚焦遗址本体、缩小调查范围而进行的前期数据沉淀，是支撑数据驱动决策、理论指导实践的关键。

       本文是与西安市隋唐长安城遗址保护中心（西安市世界遗产监测管理中心）联合开展科研攻关项目专栏的分支，详细记录从各项数据（SAR数据、DEM数据、精密轨道）下载到SARscape做DInSAR的科研工作全流程，感谢各位读者的阅读。

      目前，免费的且较为容易获取的SAR数据即Sentinel-1A（哨兵一号）数据，Sentinel-1原有2颗人造卫星，但Sentinel-1B出现故障，现已停止运行。S-1A现仍正常运转且不断更新，但美中不足的是Descending（降轨）数据在全球范围覆盖不全。因SARscape的DInSAR所需为L1 Single Look Complex（SLC）数据，能检索到研究区范围内最新的数据停更于02/20/2015,22:47:05。

       原则上，精确判断地表沉降变化，需结合Ascending（升轨）与Descending（降轨）两组数据的结果。因S-1A卫星对地观测为右摄，卫星沿轨道自南向北（S-N）飞行为升轨，对地观测即为自西向东；若沿轨道自北向南（N-S）飞行则为降轨，对地观测则变为自东向西。常规上，一个已知结果不能解释三个未知变量（即空间上X、Y、Z的变化），但好在西安为汾渭平原城市，单轨数据就可以很好的诠释相关内容。所以，即使研究区范围内缺失S-1A的Descending（降轨）数据，也并不影响通过Ascending（升轨）数据来进行沉降研究。

       当前，有诸多国内外平台支持下载Sentinel-1A（哨兵一号）数据，且有大量下载教程。本文在这里将几个主要的下载链接放在下方，对如何下载不再进行赘述，如有需求请自行搜索。

Sentinel-1A   数据下载（ASF）

Sentinel-1A   数据下载（欧空局）

Sentinel-1A   数据下载（中科院计算机网络信息中心）

Sentinel-1A（哨兵一号）数据下载（USGS）

       原则上，在SARscape的DInSAR中精密轨道数据不是必须的。一般SAR都会有自带轨道数据据，但精度较低，精密轨道一般是在S-1A数据发布后的一周至一个月之间进行公布。对应所下载的S-1A数据时段找到相应精密轨道文件进行下载。每个平台检索精密轨道文件的方式有所不同，本研究所用，以欧空局的POD Hub与Index of /aux_poeorb/两个平台检索的精密轨道文件为主。

两景Sentinel-1A（哨兵一号）数据：

对应的两组精密轨道文件： 

     本文将精密轨道文件下载链接放在下方，下载教程网上也有诸多详细讲解，如有需求请自行搜索。

S-1A  精密轨道文件下载 （欧空局）

S-1A  精密轨道文件下载 （Index of/aux_poeorb/）

      DEM，即数字高程模型（Digital Elevation Model）。目前，DEM数据产品陆地分辨率有1km、90m、30m、12.5m、5m等精度的数据，产品内容较为丰富；并不是产品分辨率越高越好，而要结合所需要的应用场景进行综合选择。因DInSAR初次需要DEM数据进行地理配准，依据经验，选择30m陆地分辨率的DEM数据是常规操作。

      针对SARscape所需要的DEM数据格式，为"_dem"结尾的文件而并非直接导入后输出的“.dat”格式。因为SARscape会识别“.”并结束文件读取，从而出现报错问题。所以，在ENVI Format导出时，请手动将后缀改为以"_dem"结尾的文件。

      本文将DEM数据下载链接，以及在ENVI 5.6中的导入与拼接工具路径放在下方，操作教程如有需求请自行搜索。

DEM 数据下载：SRTM Tile Grabber

       一景SAR数据很大，一般在4.5至7G之间。如果电脑性能不足处理起来会耗费很长时间，为提高运算处理效率、减少波段运算干扰，建议提前绘制好研究区范围的矢量文件，SARscape中可以导入.shp、.Kml、.Kmz格式的矢量作为裁剪范围。

       目前，主流的GIS软件都支持矢量文件的绘制，这一步也很简单在这里不进行细致描述；本文因研究对象是西安市隋大兴唐长安城遗址、西安市辖区境内的世界文化遗产以及隋唐时期不可移动文物要素，所以研究区矢量范围则按照行政区划直接进行选择。当然，研究区制定这一步不是必须的，如果没有研究范围对S-1A数据进行裁剪，系统会默认处理整个数据文件。

       本文将行政区划下载链接放在下方，下载文件请在QGIS中进行格式转换；如有特殊边界矢量需求，可自行在ArcGIS、QGIS等软件中进行绘制，操作教程如有需求请自行搜索。

全国行政区边界数据下载：DataV.GeoAtlas

      首先，选择适用于Sentinel-1的全套系统参数。路径：/SARscape/Preferences/Preferences specific，在左上的Load Preferences中选择Sentinel TOPSAR（IW-EW），在弹出的对话框面板选择“是”，参数面板设置点击“OK”。

      其次，是关于精密轨道文件的批量导入设置，DInSAR只需要两景数据，但SBAS-InSAR及PS则需要至少21景数据。若在导入S-1A数据时逐一选择匹配精密轨道文件不仅效率低下，也更难避免人工导入选择出错的情况，所以在运算初期设置好精密轨道文件的批量导入设置尤为关键。工具路径：/SARscape/Preferences/Preferences common，在选项面板的Directories and batch file name下Sentinel-1 auxiliary directory选项中设置orbit批量导入目录，这里文件路径不能出现中文，且存放精密轨道的文件夹名称必须是“orbit”。

      在进行DInSAR数据正式处理之前，首先要对S-1A数据进行导入，目的在于将SAR数据处理为软件可读取的相应格式，即将数据处理为pwr与slc，而不再是后缀为.xx,而是_xx的文件。工具路径：/SARscape/Import Data/SAR Spaceborne/Single Sensor/SENTINEL-1；在选择导入SAR数据时，应在解压文件下选择safe文件进行导入，精密轨道文件因前期做好了orbit设置，所以会随之自动导入，只需在导入后去文件夹检索slc_list.sml文件是否存在即可；在Import Sentinel-1面板Optional Files下的Area of Interest in Geographic Coordinates中导入之前下载或绘制好的研究区范围矢量数据；勾选自动命名，在Parameters下将Rename Output Using Parameters选择为True；最后，在Output Files选择设置输出文件夹地址即可完成数据导入操作。

工具路径：/SARscape/Interferometry/Interferometric Tools/Baseline Estimation

       原则上从本处开始，就要对Ascending（升轨）与Descending（降轨）分开计算，本文将使用升轨数据。

       首先，在做Baseline Estimation基线估算之前，需要定义主影像（master）与从影像（slave）常规的，早期数据为主（master）晚期数据为辅（slave）。因为，在相位转形变时需要*负号，若早期为slave数据，计算结果则会呈现相反态势，即原本应该为沉降的区域会变为抬升。

       本文在Baseline Estimation面板中Input Files下master为20230109T，slave则为20240209T。Optional Files和Parameters均为默认，点击Exec运行即可。

1.Baseline Estimation 基线估算报表

由基线估算报表可知：

2.InSAR

3.DInSAR

4.理论高程精度

5.理论形变精度 

路径：/SARscape/Interferometry/Phase Processing

       DInSAR操作流程比较简单，按照Phase Processing的操作顺序执行即可，但需注意些许参数面板的设置 。

1.1 - Interferogram Generation

       这一步即为生成干涉图，本环节耗时较久。在Interferogram Generation面板Input Files下参照Master为早期主影像，Slave为晚期辅影像的原则，输入两景按研究区范围裁剪过的SLC数据。产出结果为经过配准和多视的两景数据残差相位图以及主、从影像的强度图。其中，在DEM/Cartographic System中导入预设好的“_dem”文件，设置多视比及其它参数的同时请注意在Parameters中将Coregistration With DEM设施为True。其余参数选项默认，Exec执行即可。

查看干涉结果

       可以看出生成结果是镜像的，因为数据起初为SAR坐标，这里不必担心因为最后需要对结果通过地理配准进行矫正，转化为WGS_1984_Mercator；该图左侧为_cc,右侧为_dint。

2. 2 - Adaptive Filter and Coherence Generation

       这一步即自适应滤波及相干性生成。在面板的Input Files中，将第一步生成的干涉结果即_dint文件导入Interferogram File中，Master及Slave会自行匹配相对应的_pwr文件。值得注意的是Parameters中的参数设置，即也可以是默认内容，但也可以自行设置更改滤波窗口的大小。滤波窗口的大小就是滤波器的大小，以像元为基础单位，比如在Goldstein win Size中设置64，那么滤波器就是64*64像元。滤波窗口的设置并非越大越好，而是根据实际干涉图的多视比例进行不断调试，从而找到最佳效果。滤波后输出的结果后缀一般为fint。

       将滤波前的_dint与滤波后的_fint数据进行比对，可以发现整个滤波后的结果成像更加顺滑，噪点也随即少了很多。 

3.3 - Phase Unwrapping

      相位解缠主要有三个步骤，即输入相干性文件_cc和滤波后的干涉文件_fint、设置面板参数以及查看解缠图像upha。

       其目的在于解决2π模糊的问题。干涉相位只能以2π为模型，所以只要变化超过2π就会开始新的循环，相位解缠是对去平和滤波后的位相进行解缠处理。解缠后的干涉图显示研究区的总形变量的总位移，即S-1A卫星对地观测的前后距离，远离即为沉降，靠近即为抬升。解缠的算法有三种：区域增长法（Region Growing）、最小费用流算法（Minimum Cost Flow）、改进的最小费用流算法（Delaunay MCF），在参数面板即Phase Unwrapping中的Principal Parameters中设置一个不太大的相关性阈值（如1、0.3或0、0.15），避免相位解缠中不连续的区域产生“相位孤岛”输出结果为_upha文件。

       左侧为滤波后的_fint文件，右侧为对其解缠后的_upha。 

4.4 - Refinement and Re-flattening

       轨道精炼与重去平环节耗时较短，按照操作面板输入_CC文件，其余会自行匹配导入。值得提到的是在操作面板的Refinement GCPs（Ground Control Points）中手动输入制作轨道精炼的地面控制点，原则上打点选择的地方最好是地形较为平坦、形变较小且远离主形变区、没有阴影遮掩的区域。插入控制点完成后，需要结合_fint来检查GCPs控制点，若有在坡度较大与形变区的，就要移动或删除；在DEM/Cartographic System中导入DEM数据进行配准；在Parameters中进行参数调控，点击Exec执行运算，并查验结果。

       得到重去平后的输出结果_reflat_fint与_reflat_upha，并显示Refinement Results结果。由报告结果可知，均方根误差为Root Mean Square error (m)= 71.6498729947，小于官方公示的100满足最低要求，但越小越好，控制在10左右为最佳。

5.5A - Phase to Height Conversion and Geocoding

       相位转形变环节操作流程比较简单，即输入_CC与_upha文件与导入DEM数据，参收都是默认的，即可输出查看结果。

       查看 5A - Phase to Height Conversion and Geocoding计算后的形变结果，可以发现已将SAR坐标系下的影像转化为WGS-84坐标系，这样更便于后面在ArcGIS中生成专题图。

       通过在ENVI左侧的Layer Manager对_dispwf_disp文件进行右键，打开工具Quick Stats，对栅格数据的数值及特征进行简单查看。

       可以通过快速统计工具Quick Stats，查看形变量的最大值与最小值，即最大值约为0.123m，最小值约为-0.096146m。这也就是研究区范围内2023年01月09日至2024年02月09日的形变结果。

       研究结果回溯，是对隋大兴唐长安城遗址2023年01月09日至2024年02月09日，地质活动变化监测的阶段性总结，也是定量分析回归定性客观描述的阶段性成果。结合DInSAR生产的_disp文件，导入以下两组shp格式的数据内容：

陕西省文物局_文物保护单位名录

      对隋大兴唐长安城遗址整体、不可移动文物及世界遗产点，进行地质活动变化的时序对比观测。可以将_disp文件导入ArcGIS中进行分析，也可以在ENVI中对数据成果进行提取。

       本文以ENVI 5.6作为基础软件，对数据成果进行分批提取 。由呈现结果可知：

       首先，隋大兴唐长安城遗址当前有接近1/2的区域在观测周期内存在明显抬升，主要分布于遗址东南部与南部，北部大域相对稳定性较好；从整体情况综合判断，抬升与沉降的变化相对不大，但在空间尺度上的分布较广且具有连续性与周期性，应对该区域分布的文化遗产及不可移动文物进行持久关注。

       其次，对市区周围的32处不可移动文物即遗址本体进行对地观测，其中有世界文化遗产4处。成果梳理如下：

注：value/mm值，仅代表空间向量上Z轴的位移变化，并不解释具体位移形变方向。

       通过DInSAR对西安市隋大兴唐长安城遗址开展的年度周期性监测，发现成因是地质内部构造断裂带活动所导致的地面产生物理空间上的形变，而地面的沉降与抬升是直接影响城市发展与文化遗产保护最为常见的缓变型不可逆地质灾害。

       西安市地处汾渭断陷盆地中部，受临潼-长安断裂带影响，导致区域内构造活动较为强烈。由于早年西安市地表水资源匮乏，从上世纪50年代末开始大量开采地下水资源以满足城市人口不断扩张的需求，从而导致地下承压含水层水位急速下降，形成多个局部沉降中心（八里村、西工大、电子城等）。自1992年起，地面沉降的现象在西安市屡见不鲜，也自此引发了相关专家学者对该研究方向的持续关注。时至2010年，西安市城墙景区管委会综合执法大队统计西安市明城墙因地质构造活动引发的裂缝共有214处。众所周知的西安大雁塔倾斜问题，是由于其位置处于地质断裂带边缘且同时处于西北方向漏斗形下沉结构周边，其本质原因还是因为地下水资源开采所导致的不均匀沉降速率过快造成了塔身倾斜。由此可以看出，地表的抬升与沉降对文化遗产起着致命性打击与不可逆的隐形危害。

       回归到隋大兴唐长安城遗址，目前共计有11条（F2至F12，其中包含F8-1）相对均匀分布的断裂带，自东北向西南方向穿透郭城、皇城、宫城遗址的大部分区域。再加之自唐以后（包含唐朝在内），有史可查震源地在陕境内的地震共计60余起（参考：陕西地方志办公室—地震纪年），应在不同程度上对隋大兴唐长安城遗址造成了不可逆的迭进式形变损害。从技术层面来看，因本文是科研先遣的前期预估工作，目的在于缩小重点研究范围、聚焦遗产单体要素，所以采用DInSAR技术进行2023年至2024年的时段性监测。如果想要深化研究结果，从时空序列的周期性上提升数据的精准度，还需结合PS＋SBAS-InSAR来互相验证数据成果的精度是否可靠，从而通过小波变换以具体分析隋大兴唐长安城遗址地面沉降与抬升的时空演化特征与影响要素，为文化遗产保护与安全性评估提供更加健全、科学、严谨的技术支持。

       地质构造活动引发的地表沉降与抬升，是文化遗产预防性保护应当着重关注的一个重点。目前，文物部门对地质构造活动引发的物质运动机制了解相对不足，主要原因是缺乏跨学科理论知识与对遗址整体与单体要素的长期有效观测，特别是对西北地区夯土遗址的形变了解。后期，仍需不断加大科技创新力度与跨学科交叉集成研究，来筑牢西北地区大遗址与不可移动文物预防性保护的坚实防线，从而守护好跨越千年的“物化记忆”历史文化遗产。