Nepal Earthquake Triggered Landslides Recognition and Deformation Analysis of Avalanches' Region
-
摘要: 2015-04-25尼泊尔Gorkha附近发生了Mw 7.8级强震，引发尼泊尔及中国西藏地区大量的山体滑坡和雪崩灾害，造成了严重的人员伤亡与经济损失。基于Landsat 8影像和亚像素相关性匹配技术识别形变失相关区域，并结合目视解译提取震害滑坡与雪崩沉积区域，统计分析了所提取的滑坡数目、面积等与地形因子的关系，同时对珠峰雪崩高发区的冰川形变进行了定量分析。实验结果表明，提取形变场中的形变失相关区域并结合目视解译验证是一种识别震害滑坡与雪崩沉积区域的有效方法；统计结果显示滑坡体随坡度分布呈现高斯特性，多发生在坡度大于30°的地形，易发于河谷两侧，且高坡度更容易触发单体面积大的滑坡体。该研究为利用Landsat 8影像数据进行震害滑坡与雪崩沉积区域提取、雪崩高发区的冰川流速监测预警工作提供了新视角。
-
关键词:
- 尼泊尔Gorkha地震 /
- 亚像素相关性匹配 /
- 滑坡 /
- 雪崩 /
- 冰川流速
Abstract: On April 25th 2015, an Mw 7.8 earthquake occurred near Gorkha, Nepal, which triggered substantial landslides and avalanches in Tibet and Nepal causing severe casualties and economic losses. In this study, based on Landsat 8 images and sub-pixel correlation technology, we acquire horizontal deformation fields of the earthquake region. Combined with the visual interpretation method, we detected the landslides and debris deposition areas of avalanches by extracting the corresponding decorrelation noises. We determined the correlation between the landslides and topography factors. Experimental results show that the extraction of decorrelation noises caused by geohazards from the deformation field, when combined with validation of visual interpretation method, is an effective method for locating the regions effected by landslides and avalanches. Furthermore, results show that landslides are more likely to occur at slope > 30° terrain, both sides of the valley, and high slope is more easily trigger larger single landslides. This study provides a new perspective on detection of earthquake-triggered landslides and debris deposition areas of avalanches, monitoring and early warning of glaciers' velocity based on the Landsat 8 satellite imagery.-
Keywords:
- Nepal Gorkha earthquake /
- sub-pixel correlation /
- landslides /
- avalanches /
- glacier velocity
-
数字高程模型(digital elevation model, DEM)是对地球表面地形地貌的离散化数字表达，其主要借助插值算法对不规则采样数据插值完成[1-2]。受仪器噪声、人为操作不当以及外界环境不理想等因素影响，DEM源数据不可避免地含有采样误差，即误差是空间采样数据的基本属性[3-6]。插值过程中，采样误差能够被传播和放大，进而显著影响DEM及其派生品的精度。为此，众多研究人员建议采用光滑插值算法代替准确插值用于曲面建模。例如，文献[7]使用克里金(Kriging)法值移除地球物理数据中的噪声，实现光滑地学曲面建模。文献[8]认为薄板样条算法具有较强光滑性，可以有效去除采样数据中误差。鉴于不规则三角网(triangulated irregular network, TIN)线性插值的理论缺陷，文献[9]借助二元二次多项式曲面加权平均实现了DEM光滑插值。考虑到航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission, SRTM)数据受噪声影响，文献[9]利用Kriging值完成SRTM重采样。为了抑制采样数据中噪声的影响，文献[10]发展了高精度曲面建模的光滑插值算法。理论上，光滑插值算法仅降低采样数据中的高程误差分量对曲面建模的影响，而忽略了水平误差分量。尽管仅移除高程误差分量可实现DEM光滑建模，但水平误差分量会导致DEM及其派生品(如地形和水系特征线)产生平面位置误差[11]，进而降低DEM在洪水风险分析[12]、山体滑坡[13]、遥感影像配准[14]等领域数字模拟结果的精度。因此，在水平位置精度要求较高的应用中，DEM建模时需同时考虑采样数据中的水平和高程误差。
总体最小二乘算法(total least squares, TLS)是近30年发展起来的一种能同时顾及系数矩阵和观测向量误差的数据处理方法[15]。在曲线(面)拟合中，TLS通过极小化所有采样点到待定曲线(面)的正交距离之和，来降低采样数据水平和高程误差分量对模拟结果的影响。大量实例分析表明，当采样点受水平和高程误差分量影响时，TLS模拟结果精度高于最小二乘算法[16-17]。实际应用中，TLS自身并不能作为插值算法直接用于DEM建模。传统插值算法中，多面函数(digital elevation model, MQ)以其高计算精度已成为DEM建模的主要技术手段[18-19]。因此，本文以MQ为基函数，借助TLS处理误差思路，发展总体最小二乘MQ算法(MQ-T)，以降低采样点水平和高程误差分量对模拟结果的影响，进而实现DEM高精度建模。
1 MQ-T原理
1.1 传统MQ
在曲面建模中，MQ基于数值逼近理论，假设任意平滑的规则或者不规则曲面均可以由简单曲面相互叠加进而逼近到任意精度等级[20]。MQ求算未知点(xi, yi)函数值的解析表达式为[21]：
$$ f\left( {{x_i},{y_i}} \right) = \sum\limits_{j = 1}^n {{w_j}\varphi \left( {{r_{ij}}} \right)} + \sum\limits_{j = 0}^m {{b_j}{p_j}\left( {{x_i},{y_i}} \right)} $$ (1) 式中，m表示多项式阶数，本文取m=1；n表示采样点数；pj(xi, yi)和bj分别表示多项式的项和对应权重系数；rij表示待求点到第j个采样点的距离；φ(r)和w分别表示映射函数和对应权重。因此，基于式(1)求算未知点函数值f(xi, yi)的前提是已知b和w。
因采样点总数小于未知参数个数(n < n+m)，直接利用最小二乘算法无法解算未知参数值。因此，可借助岭回归推导MQ基本方程组。岭回归实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法[22]。岭回归目标函数是在最小二乘方法目标函数基础上加一正则项。基于岭回归获得MQ目标函数为：
$$ \left\{ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}}} F\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}}} \right) = \frac{1}{2}{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}} + \frac{c}{2}{\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{e}}\\ {\rm{s}}.\;{\rm{t}}.\;{f_i} = {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} + {e_i},i = 1 \cdots n \end{array} \right. $$ (2) 式中，fi=f(xi, yi)；w=[w1…wn]T，b=[b0…bm]T；e=[e1…en]T，表示采样点高程误差；φi=[φ(ri1)]…φ(rin)]; pi=[p0(xi, yi)…pm(xi, yi)]；c表示光滑参数。
为了求式(2)极小值，由拉格朗日乘数法构建辅助函数为：
$$ \begin{array}{l} L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right) = \frac{1}{2}{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}} + \frac{c}{2}{\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{e}} + \\ \sum\limits_{j = 1}^n {{\alpha _j}\left( {{f_i} - {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} - {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} - {e_i}} \right)} \end{array} $$ (3) 式中，α表示拉格朗日算子。
对L(w, e, b, α)求关于w、e、b和α的导数并等于零得：
$$ \left\{ \begin{array}{l} \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{w}}}} = 0 \Rightarrow \mathit{\boldsymbol{w}} = \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{\alpha }}_i}\mathit{\boldsymbol{\varphi }}_i^{\rm{T}}} \\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\mathit{\boldsymbol{e}}_i}}} = 0 \Rightarrow c{\mathit{\boldsymbol{e}}_i} = {\mathit{\boldsymbol{\alpha }}_i},i = 1 \cdots n\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{b}}}} = 0 \Rightarrow \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{p}}_i}{\mathit{\boldsymbol{\alpha }}_i}} = 0\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\mathit{\boldsymbol{\alpha }}_i}}} = 0 \Rightarrow {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + \\ \;\;\;\;\;\;{\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} + {e_i} = {f_i},i = 1 \ldots n \end{array} \right. $$ (4) 式(4)中消去w和e得：
$$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}} + {\rm{diag}}\left( {\frac{1}{c}} \right)}&\mathit{\boldsymbol{P}}\\ {{\mathit{\boldsymbol{P}}^{\rm{T}}}}&\mathit{\boldsymbol{0}} \end{array}} \right]\left[ \begin{array}{l} \mathit{\boldsymbol{\alpha }}\\ \mathit{\boldsymbol{b}} \end{array} \right] = \left[ \begin{array}{l} \mathit{\boldsymbol{f}}\\ \mathit{\boldsymbol{0}} \end{array} \right] $$ (5) 式中，P=(pij)，pij=pj(xi, yi)；α=[α1…αn]T；f=[f1…fn]T；K=(kij)为MQ径向基函数，由核方法推导而来，即kij=k(rij)=φiφjT。文献[23]分析认为，三次曲面作为核函数其精度优于Hardy所采用的双曲面法，故本文采用k(r)=1+r3。
将$ \mathit{\boldsymbol{w}} = \sum\limits_{i = 1}^n {{\alpha _i}\mathit{\boldsymbol{\varphi }}_i^{\rm{T}}} $代入式(1)得：
$$ \begin{array}{l} f\left( {{x_i},{y_i}} \right) = \sum\limits_{j = 1}^m {{\alpha _j}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{\varphi }}_j^{\rm{T}}} + \sum\limits_{j = 0}^m {{b_j}{p_j}\left( {{x_j},{y_i}} \right)} \\ \;\;\;\;\; = \sum\limits_{j = 1}^m {{\alpha _j}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{\varphi }}_j^{\rm{T}}} + \sum\limits_{j = 0}^m {{b_j}{p_j}\left( {{x_j},{y_i}} \right)} \end{array} $$ (6) 其矩阵形式为：
$$ \mathit{\boldsymbol{f}} = \mathit{\boldsymbol{K\alpha + Pb}} $$ (7) 因此，基于MQ插值计算时必须首先解算方程组(5)获得α和b的值，然后基于式(6)或式(7)求算未知点函数值。由以上推导可见，MQ计算中仅考虑采样点高程误差对模拟结果的影响，而忽略了水平误差分量。当采样点同时含有水平和高程误差时，此方案势必使得MQ插值结果产生平面位置误差。
1.2 MQ-T
理论上，式(1)可看作为原始特征空间中的超平面方程，其法向量为(w －1)。因此，基于整体最小二乘算法思想，点的平面位置误差为点到超平面的垂直距离。对于MQ而言，采样点平面位置误差可表达为：
$$ {\varepsilon _i} = \frac{{{f_i} - {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} - {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}}}}{{\sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} }} $$ (8) 因此，MQ-T的目标函数表达为：
$$ \left\{ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}}} F\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{\varepsilon }}} \right) = \frac{1}{2}{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}} + \frac{c}{2}{\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}\mathit{\boldsymbol{\varepsilon }}\\ {\rm{s}}.\;{\rm{t}}.\;{f_i} = {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} + {\varepsilon _i}\sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} ,i = 1 \cdots n \end{array} \right. $$ (9) 类似MQ推导过程，对式(9)由拉格朗日乘数法构建辅助函数，求其关于w、ε、b、α的导数并等于零得：
$$ \left\{ \begin{array}{l} \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{w}}}} = 0 \Rightarrow \mathit{\boldsymbol{w}} = \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{\alpha }}_i}\mathit{\boldsymbol{\varphi }}_i^{\rm{T}}} \\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{\varepsilon }},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\varepsilon _i}}} = 0 \Rightarrow {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + \sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} {\varepsilon _i}\sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} ,i = 1 \cdots n\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{b}}}} = 0 \Rightarrow \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{p}}_i}{\mathit{\boldsymbol{\alpha }}_i}} = 0\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\mathit{\boldsymbol{\alpha }}_i}}} = 0 \Rightarrow {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + \sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} = {f_i},i = 1 \cdots n \end{array} \right. $$ (10) 对式(10)消去w和ε得：
$$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}} + {\rm{diag}}\left( {\frac{1}{c}} \right)}&\mathit{\boldsymbol{P}}\\ {{\mathit{\boldsymbol{P}}^{\rm{T}}}}&\mathit{\boldsymbol{0}} \end{array}} \right]\left[ \begin{array}{l} \mathit{\boldsymbol{\alpha }}\\ \mathit{\boldsymbol{b}} \end{array} \right] = \left[ \begin{array}{l} \mathit{\boldsymbol{f}}\\ \mathit{\boldsymbol{0}} \end{array} \right] $$ (11) 因wTw=αTKα，则式(11)迭代表达式为：
$$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}} + {\rm{diag}}\left( {\frac{{1 + {{\left( {{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}} \right)}^{\left( i \right)}}K{\alpha ^{\left( i \right)}}}}{c}} \right)}&\mathit{\boldsymbol{P}}\\ {{\mathit{\boldsymbol{P}}^{\rm{T}}}}&\mathit{\boldsymbol{0}} \end{array}} \right]\left[ \begin{array}{l} {\mathit{\boldsymbol{\alpha }}^{\left( {i + 1} \right)}}\\ {\mathit{\boldsymbol{b}}^{\left( {i + 1} \right)}} \end{array} \right] = \left[ \begin{array}{l} \mathit{\boldsymbol{f}}\\ \mathit{\boldsymbol{0}} \end{array} \right] $$ (12) 式中，i为迭代次数。
比较式(5)和式(12)得出，MQ为MQ-T特殊情况，即当αTKα=0时的情况。因此，式(12)计算时，应先借助式(5)求算α，然后基于式(12)迭代计算，直至α和b收敛。
2 DEM建模算法精度分析
本文分别以数值实验和实例分析验证MQ-T计算精度。数值实验中，以高斯合成曲面为研究对象，分别设计受不同误差分量污染的采样数据，基于MQ-T曲面建模，并与传统MQ方法计算结果比较，验证模型的去噪性。实例分析中，以全站仪获取的采样数据为研究对象，借助MQ-T构建测区DEM，并将其与传统插值算法(如IDW，Kriging和ANUDEM)计算结果比较，验证模型的实用性。
2.1 数值实验
高斯合成曲面表达式为：
$$ \begin{array}{l} z\left( {x,y} \right) = 3{\left( {1 - x} \right)^2}{{\rm{e}}^{ - x2 - {{\left( {y + 1} \right)}^2}}} - \\ 10\left( {\frac{x}{5} - {x^3} - {y^5}} \right){{\rm{e}}^{ - {x^2} - {y^2}}} - \frac{1}{3}{{\rm{e}}^{ - {{\left( {x + 1} \right)}^2} - {y^2}}} \end{array} $$ 计算区域为[－3, 3]×[－3, 3]，模拟区域网格数为101 × 101。从高斯合成曲面中随机采集31×31个采样点用于曲面建模，采样点分别受如下误差影响。
例1：仅受高程误差分量影响。其中，误差来源于正态分布N (0, 0.02)。
例2：仅受水平误差分量影响。其中，x和y坐标误差均来源于N (0, 0.01)。因此，例1中高程误差分量方差和例2中水平误差分量方差相等。此方案可用于比较水平和高程误差对模拟结果精度影响的程度。
例3：同时受水平和高程误差分量影响。其中，水平和高程误差分别来源于正态分布N (0, 0.01)和N (0, 0.02)。
实例中分别采用中误差(root mean squared error，RMSE)和平均误差(mean error, ME)作为精度指标比较MQ-T和MQ模拟结果精度。中误差和平均误差计算的公式为：
$$ {\rm{RMSE}} = \sqrt {\frac{{\sum\limits_{i = 1}^k {{{\left( {{z_i} - {{\hat z}_i}} \right)}^2}} }}{k}} $$ (13) $$ {\rm{ME}} = \frac{{\sum\limits_{i = 1}^k {\left( {{z_i} - {{\hat z}_i}} \right)} }}{k} $$ (14) 式中，zi和${{\hat z}_i} $分别表示第i个检核点的真实值和模拟值；k表示检核点个数。
由表 1可见，当采样点仅受高程误差分量影响时(例1)，两种方法模拟结果相当，表明MQ-T对高程误差具有较好的光滑性。当采样点仅受水平误差影响时(例2)，两种方法模拟结果精度均低于例1，且MQ-T误差小于MQ误差，表明MQ-T能较好消除采样点水平误差分量对曲面模拟的影响。当采样点水平和高程坐标分量均受误差影响时(例3)，两种方法模拟结果精度均低于例1和例2，表明采样点精度越低，模拟结果精度越低；MQ-T中误差明显小于MQ，但前者平均误差略大于后者。综上所述，当采样点含有水平误差分量时，以中误差为误差指标，MQ-T模拟结果精度高于MQ。
表 1 3种误差分布条件下，MQ和MQ-T模拟结果精度比较Table 1. Accuracy Comparison Between the Classical MQ and MQ-T Under Different Error Distributions算例 误差分布 方法 RMSE ME 水平分量 高程分量 例1 - N (0, 0.02) MQMQ-T 0.062 70.062 1 0.0070.007 例2 N (0, 0.01) - MQMQ-T 0.147 80.139 3 -0.003-0.002 例3 N (0, 0.01) N (0, 0.02) MQMQ-T 0.158 70.142 7 -0.014-0.015 2.2 实例分析
测区位于山东省济南市，测区面积为3.36 km2，最小和最大高程分别为109.6 m和364.5 m。利用全站仪获取该测区5 394个采样点，平均点间距(sample interval, SI)为25 m。基于全站仪测角和测距精度指标以及数据采集时各种观测条件，采样点平面和高程精度分别为5 cm和10 cm。为了验证采样点密度对DEM精度的影响，将原始采样数据依次抽稀，每次抽稀后数据点数和对应样点间距如表 2所示。除了MQ-T和MQ外，还采用传统插值算法如IDW、Kriging以及ANUDEM进行10 m分辨率DEM建模。为了验证DEM精度，借助实时动态定位软件RTK获取了该测区500个检核点(见图 1)。各种插值算法构建DEM精度如表 3所示。
表 2 采样点与采样点平均点间距对应关系Table 2. Relationship Between Number of Sample Points and Sample Interva采样点数 点间距/m 2 100 40 2 742 35 3 721 30 5 394 25 表 3 各种插值算法模拟结果中误差/mTable 3. RMSE of Different Interpolation Methods/m点间距 IDW Kriging ANUDEM MQ MQ-T 40 3.93 2.43 2.68 2.38 2.26 35 3.41 2.37 2.55 2.24 2.11 30 3.17 2.22 2.43 2.09 2.06 25 2.69 1.84 2.28 1.73 1.73 平均 3.30 2.22 2.49 2.11 2.04 由表 3可见，随着平均点间距增大，各种插值算法精度均有所下降。具体而言，IDW、Kriging、ANUDEM和MQ降幅分别为46%、32%、15%和38%；MQ-T精度降幅为31%。结果表明，ANUDEM模拟结果受采样点间距影响最小，而IDW影响最大。不管采用何种采样点间距，MQ-T计算精度均高于其他插值方法。MQ-T平均计算精度分别是IDW的1.61倍、Kriging法的1.09倍、ANUDEM法的1.22倍以及MQ的1.03倍。
各种插值方法基于30 m平均点间距采样点构建的10 m分辨率DEM山体阴影图如图 2所示。由图 2可见，IDW结果最粗糙，这主要是因为IDW为准确插值算法，无法光滑去噪。ANUDEM法结果最光滑，去噪效果明显，但导致很多微地形特征信息丢失(如图 2(c)中圆圈标注的区域)。Kriging法、MQ法和MQ-T法结果相似，但前两种方法模拟结果有零散坑洼和凸起(如图 2(b)和2(d)中方框标注区域)。相比较而言，MQ-T法模拟结果稍好于其他方法。但相比Kriging法，MQ-T法结果左上角的山脊有峰值削平现象(如图 2(e)中圆圈标记区域)。
![]()
图 2 各种插值算法构建DEM的山体阴影图Figure 2. Hillshades of IDW, Kriging, ANUDEM, MQ and MQ-T Algorithms3 结语
为了降低采样数据中的水平和高程误差分量对DEM建模精度的影响，本文以较高插值精度的多面函数(MQ)为基函数，借助整体最小二乘算法(TLS)思路，发展了整体最小二乘MQ算法(MQ-T)。数值实验分析表明：(1)当采样点仅受高程误差分量影响时，MQ-T法计算结果精度与MQ法相当；(2)当采样点受水平误差分量影响，以中误差为误差指标，MQ-T计算结果精度高于MQ法。实例分析中，分别基于MQ法和MQ-T法构建测区DEM，并将计算结果与传统插值算法(如IDW、Kriging、ANUDEM)比较。结果分析表明：(1)各种插值算法结果精度均随采样点密度的降低而逐步降低，其中，ANUDEM精度变化最小，IDW精度变化最大；(2)不管采用何种采样点间隔，MQ-T精度始终高于其他插值算法。
-
![]()
图 1 震区构造背景与Landsat 8影像覆盖图
Figure 1. Background Setting and Coverage of Landsat 8 for the 2015 Nepal Earthquake
![]()
图 5 吉隆、聂拉木地区不同坡度水平滑坡体数目、总面积与平均面积分布直方图
Figure 5. Landslides' Quantity, Gross and Average Area Distribution Histogram of Jilong and Nyalam in Different Slope Levels
![]()
图 6 郎塘形变场和雪崩沉积区提取图
Figure 6. Extraction of Deformation Fields and Avalanches' Debris Areas in Langtang
表 1 本文实验所采用的Landsat 8影像信息
Table 1 Information of Landsat 8 Images Used in This Study
影像对 时段 轨道号 分幅号 获取时间 云量/% 影像对1 震前 141 41 2014-03-26 2.4 震后 141 41 2015-04-30 30.2 影像对2 震前 141 40 2014-01-05 4.4 震后 141 40 2015-04-30 30.2 影像对3 震前 140 41 2015-03-22 11.8 震后 140 41 2015-05-25 25.6 影像对4 震前 140 41 2014-03-19 17.4 震后 140 41 2015-05-22 13.1 
下载: 导出CSV
表 2 吉隆山体滑坡编目表
Table 2 Landslides Catalog of Jilong
序号 滑坡位置 滑坡体面积/m2 技术监测手段 高程/m 坡向/(°) 坡度/(°) 东经/(°) 北纬/(°) SPC 目视判读 J01 85.363 28.360 163 646 √ √ 2 615 289.6 52.8 J02 85.359 28.354 29 052 √ √ 2 370 268.7 43.5 J03 85.337 28.376 37 580 √ √ 2 490 80.4 30.2 J04 85.324 28.389 49 045 √ √ 2 737 312.1 45.9 J05 85.305 28.401 5 599 √ √ 2 837 359.1 37.3 J06 85.304 28.406 33 891 √ √ 2 809 28.9 40.3 J07 85.350 28.364 2 603 √ √ 2 543 10.6 47.5 J08 85.345 28.369 11 552 √ √ 2 415 165.9 38.9 J09 85.337 28.374 5 967 √ √ 2 450 102.1 33.3 J10 85.340 28.381 7 431 线状滑坡× √ 2 839 208.5 53.9 J11 85.337 28.374 1 295 小滑坡× √ 2 603 24.7 27.1 J12 85.355 28.364 680 小滑坡× √ 2 472 120.0 45.7 J13 85.370 28.376 6 789 √ √ 2 925 212.4 44.8 J14 85.358 28.370 42 480 √ √ 2 439 196.8 52.6 J15 85.355 28.372 1 057 小滑坡× √ 2 550 172.0 41.4 J16 85.308 28.404 2 263 小滑坡× √ 2 804 191.8 33.3 J17 85.303 28.408 1 728 小滑坡× √ 2 805 49.1 42.9 J18 85.312 28.399 3 628 小滑坡× √ 2 741 172.1 25.6 J19 85.358 28.411 351 小滑坡× √ 2 984 135.1 37.0 J20 85.364 28.382 21 275 √ √ 3 152 256.7 48.2 J21 85.365 28.382 17 310 √ √ 3 152 256.7 42.7 注：√表示可以监测到对应滑坡体，×表示无法监测到对应滑坡体。
下载: 导出CSV
表 3 聂拉木山体滑坡编目表
Table 3 Landslides Catalog of Nyalam
序号 滑坡位置 滑坡体面积/m2 高程/m 坡向/(°) 坡度/(°) 东经/(°) 北纬/(°) N01 85.997 28.090 46 565 3 457 123.3 37.6 N02 85.998 28.087 7 153 3 264 177.7 31.7 N03 85.990 28.084 13 099 3 769 166.3 60.8 N04 85.997 28.085 6 360 3 238 270 26.9 N05 85.991 28.046 38 639 3 265 274.9 57.1 N06 85.994 28.039 53 518 3 049 171.8 32.8 N07 85.981 28.032 51 028 3 157 100.9 68.8 N08 85.980 28.021 12 191 2 602 106.2 45.4 N09 85.963 27.983 309 600 2 375 157.2 60.8 N10 86.039 27.975 381 165 3 462 106.5 42.2
下载: 导出CSV
-
[1] Wikipedia. April 2015 Nepal Earthquake[OL]. https://en.wikipedia.org/wiki/April 2015 Nepal earthquake, 2016
[2] Feng G C, Li Z W, Shan X J, et al. Geodetic Model of the 2015 April 25 Mw 7.8 Gorkha Nepal Earthquake and Mw 7.3 Aftershock Estimated from InSAR and GPS Data[J]. Geophysical Journal International, 2015, 203(2):896-900 doi: 10.1093/gji/ggv335
[3] 陶舒. 汶川地震滑坡遥感信息提取及灾害危险性评价研究[D]. 北京: 首都师范大学, 2009 Tao Shu. Information Extraction of Landslides Induced by the Wenchuan Earthquake and Susceptibi-lity Assessment Using Remotely Sensed Data[D]. Beijing: Capital Normal University, 2009
[4] 许冲, 徐锡伟, 沈玲玲, 等. 2014年鲁甸MS 6.5地震触发滑坡编录及其对一些地震参数的指示[J].地震地质, 2014, 36(4):1186-1203 https://www.researchgate.net/profile/Chong_Xu/publication/287556048_Inventory_of_landslides_triggered_by_the_2014_MS65_Ludian_earthquake_and_its_implications_on_several_earthquake_parameters/links/580ea47a08ae7525273d29d3.pdf?origin=publication_list Xu Chong, Xu Xiwei, Shen Lingling, et al. Inventory of Landslides Triggered by the 2014 MS 6.5 Ludian Earthquake and Its Implications on Several Earthquake Parameters[J].Seismology and Geology, 2014, 36(4):1186-1203 https://www.researchgate.net/profile/Chong_Xu/publication/287556048_Inventory_of_landslides_triggered_by_the_2014_MS65_Ludian_earthquake_and_its_implications_on_several_earthquake_parameters/links/580ea47a08ae7525273d29d3.pdf?origin=publication_list
[5] 武雪玲, 沈少青, 牛瑞卿. GIS支持下应用PSO-SVM模型预测滑坡易发性[J].武汉大学学报·信息科学版, 2016, 41(5):665-671 http://ch.whu.edu.cn/CN/abstract/abstract5447.shtml Wu Xueling, Shen Shaoqing, Niu Ruiqing. Landslide Susceptibility Prediction Using GIS and PSO-SVM[J]. Geomatics and Information Science of Wuhan University, 2016, 41(5):665-671 http://ch.whu.edu.cn/CN/abstract/abstract5447.shtml
[6] Lacroix P, Zavala B, Berthier E, et al. Supervised Method of Landslide Inventory Using Panchromatic SPOT5 Images and Application to the Earthquake-Triggered Landslides of Pisco (Peru, 2007, Mw8.0)[J]. Remote Sensing, 2013, 5(6):2590-2616 doi: 10.3390/rs5062590
[7] Leprince S, Barbot S, Ayoub F, et al. Automatic and Precise Orthorectification, Co-registration and Sub-pixel Correlation of Satellite Images, Application to Ground Deformation Measurements[J]. IEEE Transactions on Geoscience and Remote Sen-sing, 2007, 45(6):1529-1558 doi: 10.1109/TGRS.2006.888937
[8] 董晓燕, 丁晓利, 李志伟, 等.一种新的SAR像素偏移量估计流程及其在同震形变监测中的应用[J].武汉大学学报·信息科学版, 2011, 36(7):789-792 http://ch.whu.edu.cn/CN/abstract/abstract592.shtml Dong Xiaoyan, Ding Xiaoli, Li Zhiwei, et al. A New Chain for SAR Pixel Offset Eestimation and Its Application in Co-seismic Deformation Measurement[J].Geomatics and Information Science of Wuhan University, 2011, 36(7):789-792 http://ch.whu.edu.cn/CN/abstract/abstract592.shtml
[9] Ding C, Feng G C, Li Z W, et al. Spatio-Temporal Error Sources Analysis and Accuracy Improvement in Landsat 8 Image Ground Displacement Measurements[J].Remote Sensing, 2016, 8(11):937-958 doi: 10.3390/rs8110937
[10] Binet R, Bollinger L. Horizontal Coseismic Deformation of the 2003 Bam (Iran) Earthquake Measured from SPOT-5 THR Satellite Imagery[J]. Geophysi-cal Research Letters, 2005, 320(2):287-294 doi: 10.1029/2004GL021897
[11] Leung T K P, Admas J. The 2005, Mw 7.6 Kashmir Earthquake, Rupture Kinematics from Sub-pixel Correlation of ASTER Images and Seismic Waveforms Analysis[J]. Earth and Planetary Science Letters, 2006, 249(3):514-528 https://thesis.library.caltech.edu/2257/
[12] 冯光财, 许兵, 单新建, 等.基于Landsat 8光学影像的巴基斯坦Awaran Mw7.7地震形变监测及参数反演研究[J].地球物理学报, 2015, 58(5):1634-1644 doi: 10.6038/cjg20150515 Feng Guangcai, Xu Bing, Shan Xinjian, et al. Coseismic Deformation and Source Parameters of the 24 September 2013 Awaran, Pakistan Mw7.7 Earthquake Derived from Optical Landsat 8 Satellite Images[J]. Chinese Journal Geophysics, 2015, 58(5):1634-1644 doi: 10.6038/cjg20150515
[13] Avouac J P, Meng L, Wei S, et al. Lower Edge of Locked Main Himalayan Thrust Unzipped by the 2015 Gorkha Earthquake[J].Nature Geoscience, 2015, 8(9):1-6 http://palgrave.nature.com/ngeo/journal/v8/n9/abs/ngeo2518.html
[14] Kargel J S, Leonard G J, Shugar D H, et al. Geomorphic and Geologic Controls of Geohazards Induced by Nepal's 2015 Gorkha Earthquake[J]. Science, 2016, 351(6269):1-12 http://ecite.utas.edu.au/107060
[15] Jane Q. Killer Landslides: The Lasting Legacy of Nepal's Quake[OL]. http://www.nature.com/news/killer-landslides-the-lasting-legacy-of-nepal-s-quake-1.19803, 2016
[16] Scherler D, Leprince S, Strecker M R. Glacier-surface Velocities in Alpine Terrain from Optical Satellite Imagery-Accuracy Improvement and Quality Assessment[J]. Remote Sensing of Environment, 2008, 112(10):3806-3819 doi: 10.1016/j.rse.2008.05.018
-
期刊类型引用(10)
1. 曾广泉，马韬，张孟希，戴妍，陈凯文，丁继辉，俞双恩，王中文. 基于无人机多光谱影像的不同施氮量水稻LAI反演方法研究. 江苏农业科学. 2024(20): 41-48 . 
百度学术
2. 高钰琪，许桂玲，冯跃华，王晓珂，任红军，由晓璇，韩志丽，李家乐. 基于冠层高光谱植被指数的水稻产量预测模型研究. 中国稻米. 2023(05): 38-44 . 百度学术
3. 彭晓伟，张爱军，王楠，赵丽，杨晓楠. 高光谱技术在土壤及适种作物的研究进展. 遥感信息. 2022(01): 32-39 . 百度学术
4. 王晓珂，刘婷婷，许桂玲，冯跃华，彭金凤，李杰，罗强鑫，韩志丽，卢苇，PHONENASAY Somsana. 基于冠层高光谱遥感的杂交水稻植被指数氮素营养诊断模型. 中国稻米. 2021(03): 21-29 . 百度学术
5. 王浩淼，宋苗语，李翔，扈朝阳，鲁任翔，王翔，马会勤. 无人机高光谱遥感监测葡萄长势与缺株定位. 园艺学报. 2021(08): 1626-1634 . 百度学术
6. 刘雅婷，龚龑，段博，方圣辉，彭漪. 多时相NDVI与丰度综合分析的油菜无人机遥感长势监测. 武汉大学学报(信息科学版). 2020(02): 265-272 . 百度学术
7. 陈晓凯，李粉玲，王玉娜，史博太，侯玉昊，常庆瑞. 无人机高光谱遥感估算冬小麦叶面积指数. 农业工程学报. 2020(22): 40-49 . 百度学术
8. 落莉莉，常庆瑞，武旭梅，杨景，李粉玲，王琦. 夏玉米叶片光合色素含量高光谱估算. 干旱地区农业研究. 2019(04): 178-183 . 百度学术
9. 张良培，刘蓉，杜博. 使用量子优化算法进行高光谱遥感影像处理综述. 武汉大学学报(信息科学版). 2018(12): 1811-1818 . 百度学术
10. 李亚妮，鲁蕾，刘勇. 基于PROSAIL模型的水稻田缨帽三角-叶面积指数模型及其应用. 应用生态学报. 2017(12): 3976-3984 . 百度学术
其他类型引用(17)
计量
- 文章访问数: 1602
- HTML全文浏览量: 184
- PDF下载量: 332
- 被引次数: 27
