杨云飞，李金平,2

(1.云南师范大学地理学部，昆明 650500；2.西部资源环境地理信息技术教育部工程研究中心，昆明 650500)

近年来城市爆炸时常发生，造成了巨大的人员伤亡和财产损失，如天津港爆炸[1]、响水爆炸[2-4]和贝鲁特爆炸[5]等。差分干涉测量(Differential InSAR，D-InSAR)技术是目前国内外在InSAR应用上最为成熟的技术，这项技术通过星载合成孔径雷达监测地球表面形变精度在厘米级甚至毫米级。Gabriel等[6]在1989年首次验证了D-InSAR技术具有监测地表微小形变的能力；Goldstein等[7]在1993年把D-InSAR技术应用于南极的冰川移动监测并获得成功，成果发表于Science。同年Massonnet等[8]利用D-InSAR技术获取了1992年加利福尼亚Landers地震同震形变场。此后，D-InSAR技术被广泛应用于由各种原因引发的地表形变的监测中，包括地震、滑坡、火山运动、冰川移动等自然现象[9-14]，还包括矿山开采、地下水抽取等人类活动[15-18]。

本文采用D-InSAR技术对爆炸产生地表形变场进行监测，同时对比天然地震产生形变场，得到爆炸产生的形变场特点，为爆炸后的重建提供形变信息支持；并基于形变值与萨道夫斯基公式构建城市爆炸形变模型。

1 研究区概况

天津港爆炸中心位于东经117.74°、北纬39.04°，从卫星影像(见图1)看，爆炸约554 m范围内设施严重损毁。响水爆炸位于东经119.76°、北纬34.34°，爆炸340 m范围内建筑损毁严重(见图2)。贝鲁特爆炸位于东经35.52°、北纬33.90°贝鲁特港口内，爆炸范围约553 m内港口设施摧毁殆尽，在约943 m范围内楼房有明显震碎物品坠落(见图3)。作为对比的天然地震选取2017年伊朗Mw7.3，其震中位于伊朗-伊拉克边界欧亚板块与阿拉伯板块碰撞交汇处的扎格罗斯褶皱冲断带，构造活动强烈地震多发(根据地球物理中板块运动模型NUVEL -1A计算结果显示，阿拉伯板块以30 mm/yr撞向欧亚板块[19])；从地貌特征来看，震中区域植被稀疏，利于合成孔径雷达进行干涉测量；以上两个因素也造就了此区域作为InSAR技术监测地震引发形变较为理想的试验场(见图4)。

图1 天津港爆炸卫星影像Fig.1 Satellite image of Tianjin Port explosion

图2 响水爆炸卫星影像Fig.2 Satellite image of Xiangshui explosion

图3 贝鲁特爆炸卫星影像Fig.3 Satellite image of Beirut explosion

2 数据及数据处理

Sentinel-1A数据来自于欧州航天局于2014-04-03发射的全新C波段SAR卫星， 它可以提供具有周期性且精确的大范围形变监测，性能比ERS-1/2、ENVISAT ASAR卫星获取的SAR影像更为完善，具有噪声小、影像干涉能力更强、监测精度更高、周期更短和免费申请获取等特点[20]。本文分别选取覆盖该爆炸、地震事件发生前后区域的一组 Sentinel-1A影像数据(SLC格式)，通过D-InSAR技术来获取爆炸产生形变场和地震产生形变场，影像参数如表1所示。

表1 研究区影像参数

D-InSAR技术其数据处理流程可分为：数据输入、干涉图生成、滤波和相干性计算、去平地效应及相位解缠、轨道精炼和重去平、相位转形变及地理编码[20]。鉴于轨道存在的系统误差、平地效应及不同地物具有不同相干性的实际特点，本次数据处理过程中通过引入外部数据和多次调整参数，获得了较为可靠结果，具体如下：数据输入中通过添加卫星精密轨道数据AUX_POEORB降低系统误差产生的影响提高数据处理精度；干涉噪声先采用距离向比方位向为1∶4的多视处理，之后选择高通滤波(Goldstein)去除相位噪声；利用STRM数据(30 m×30 m)外部数字高程模型(DEM)去除地形相位产生的影响；采用最小费用流法进行相位解缠，将缠绕的相位校正到高相干性参考点上，更好地还原地表形变引起的真实相位变化；在估算和去除残余的恒定相位后，轨道误差依然存在，由此选取GCP控制点，对干涉图进行重新估算，减小误差；经过以上处理获取的纯净相位信息转化为形变位移信息(包括生成雷达视线方向形变、垂直方向形变)并进行地理编码得到实际应用中对应的地表形变。

3 爆炸形变特征分析

天津港爆炸中心附近分别出现了一个较大沉降和隆升区域，其东向出现两个较大的沉降漏斗区域，整体监测区域垂直向形变范围在-3.7 ～ 4.8 cm(见图5)；响水爆炸中心附近出现了一个较小沉降区域但整体主要表现为隆升，其东向形变整体较西向剧烈，整体监测区域垂直向形变范围在-4.5～3.8 cm(见图6)；贝鲁特港爆炸中心附近主要表现为沉降但有小片区域表现为隆升，其东向主要表现为隆升西向主要表现为沉降，整体监测区域垂直向形变范围在-4.9～1.8 cm(见图7)。从天津港爆炸、响水爆炸及贝鲁特港爆炸雷达视线向形变场和垂直形变场综合来看：爆炸中心周边出现有密集且不连续(隆升与沉降交错)的形变，结合爆炸前后卫星影像地表隆升与沉降交错区域多发于建筑物密集且损毁严重区域，符合实际中建筑倒塌现象。由形变判断发生建筑损毁区域，天津港爆炸损毁范围半径约在2 034 m；响水爆炸损毁半径约在830 m；贝鲁特港爆炸损毁半径约在2 100 m。通过对比天然地震产生形变场(见图8)，地表爆炸所引发的地表形变更多的表现出一种杂乱的不连续波动。

图6 响水爆炸形变场Fig.6 Deformation field of Xiangshui explosion

图7 贝鲁特港口爆炸形变场Fig.7 Deformation field of Beirut Port explosion

图8 地震形变场Fig.8 Earthquake deformation field

4 城市爆炸形变模型构建

综合爆炸前后遥感影像目视解译及爆炸形变场(InSAR)获取信息，如表2所示。由表2发现从遥感影像目视解译能较为准确地获取爆炸源周边建筑物倒塌的范围(范围内损坏等级为完全破坏)；综合InSAR技术监测到的爆炸形变范围，可以获取形变边界至建筑物倒塌边界范围内的建筑物受损情况(范围内损坏等级为次严重破坏至基本无破坏)。

表2 遥感方式获取爆炸影响信息

由于爆炸发生在日常生产运输的地方(化工园区、港口)可认为爆炸是在普通地面发生，在此前提下爆炸产生的空气冲击波超压可以运用萨道夫斯基公式[21]计算:

(1)

式中：Δp为空气冲击波的超压值,kPa；m为TNT炸药的质量，kg；R为爆炸半径，m。

由公式(1)结合表2遥感手段获取的边界范围信息计算空气冲击波的超压值(见表3)。

表3 爆炸在遥感获取边界范围下超压值

在已知爆炸当量情况下，通过遥感手段可快速获取损毁范围和形变范围，即在较短的时间里获得损毁严重范围和次损毁严重范围内空气冲击波的超压值。

由于D-InSAR形变监测具有临界形变梯度[22]，由于临近爆炸中心，破坏等级在完全破坏范围内存在大量区域真实形变超过D-InSAR形变监测能力范围，因此本文在城市爆炸形变模型构建中不考虑完全破坏区域；通过获取的次严重区域至基本不受影响区域内空气冲击波的超压值及范围内形变值，以形变值作为自变量，超压值作为因变量，通过分析3次爆炸形变值、超压值曲线发现3次爆炸具有一致的曲线特征(见图9)；通过对3次爆炸形变曲线拟合得到基于3次爆炸的爆炸形变模型：y=-0.160x3+0.273x2-1.854x+13.150。

图9 次严重区域至基本不受影响区域形变—空气冲击波超压模型Fig. 9 Deformation from severe area to basically unaffected area-air shock wave overpressure model

5 结论

1)通过获取3次爆炸事件所引发的地表形变场，发现爆炸所引发的形变相较于天然地震具有不连续波动的特点。

2)结合爆炸前后卫星影像发现形变隆升区域为建筑密集区，以此为依据可通过形变隆升的区域评估爆炸带来的毁伤范围。

3)对于较大沉降区域，可能处于不稳定的地质区，灾害后重建时需重点评估。

4)通过总结3次爆炸形变与空气冲击波超压值曲线特点，构建了次严重区域至基本不受影响区域形变—空气冲击波超压模型。