一种融合不同分辨率SAR数据监测地表形变的方法与流程

本发明涉及一种利用合成孔径雷达监测地表形变的方法，具体是一种融合不同分辨率SAR数据监测地表形变的方法。

背景技术：

过度开采地下水、煤炭等资源引起的地表沉降已经成为我国主要的地质灾害之一。地表的沉降容易造成房屋开裂、道路塌陷，在海拔较低的沿海地区还容易引起洪涝。因此要全面地查明地面沉降的成因、现状以及发展趋势。

相比传统的观测技术手段，在合成孔径雷达(Synthetic Aperture Radar，SAR)基础上发展起来的合成孔径雷达差分干涉(Differential InSAR，DInSAR)技术，具有高分辨率、高精度、广覆盖的优势。尤其是在DInSAR技术上进一步发展的时序InSAR(Multi temporal InSAR，MT-InSAR)技术，如PS、SBAS等，通过对多景SAR数据的综合分析，能够有效抑制大气噪声以及时空失相干的影响，从而实现厘米甚至毫米级的地表沉降监测精度，并且能够展现地表沉降在时间序列上的演化过程。目前，DInSAR与时序InSAR技术已经广泛用于火山、地震、滑坡、城市地表沉降等形变监测中。

由于SAR数据中的相位是缠绕的，DInSAR与MT-InSAR技术需要通过相位解缠获取地表形变信息。为准确进行相位解缠，要求SAR数据的时间采样间隔τ满足Nyquist(奈奎斯特)采样定理，即：

式中，λ表示SAR数据的波长，Δv表示相邻像元或相邻相干点的形变速率之差。可以看出，Δv越大，要求采样间隔τ越短，即在一定的时间段内需要更多的SAR数据。

然而，我国的地面沉降情况越来越严重，部分地区的地表年沉降速率超过150毫米。但是，许多城市并没有积累足够多的高分辨率的SAR数据，实际应用中能够使用的高分辨率SAR数据也并不总是很理想。最为常见的问题就是：研究区域内可用的高分辨率SAR数据量少，并且时间分布不均匀，干涉图相干质量比较差，难以控制相邻像元及相邻相干目标点间的形变梯度，降低了形变监测精度。而同一地区的低分辨率SAR数据无法被利用，造成了数据的浪费；在地表形变剧烈的区域，该现象尤为严重。高分辨率SAR数据量的不足，严重制约了DInSAR与时序InSAR技术的实际应用能力。

技术实现要素：

针对上述现有技术存在的问题，本发明提供一种融合不同分辨率SAR数据监测地表形变的方法，可利用不同分辨率的SAR数据对地表形变进行监测，只需较少的高分辨率SAR数据就可具有较高的地表形变的监测精度，进而降低监测成本。

为了实现上述目的，本发明采用的技术方案是：一种融合不同分辨率SAR数据监测地表形变的方法，其具体步骤为：

第一步，利用低分、长波SAR数据(即低分辨率、长波的SAR数据)，生成差分干涉相位：

获取覆盖同一地区的低分、长波SAR数据，组成干涉对；然后利用该地区的DEM(即数字高程模型)数据模拟该地区的地形相位，并从每个干涉对中去除该相位，得到差分干涉相位；

第二步，提取低通形变分量Ⅰ：

根据低分、长波SAR数据量，利用DInSAR或时序InSAR方法从步骤一得到的差分干涉相位中提取地表形变相位；将该地表形变相位视为总体形变在时间与空间上的低通部分，即低通形变分量Ⅰ；

第三步，配准低分、长波SAR数据与高分、短波SAR数据(即高分辨率、短波的SAR数据)：

利用DEM数据将低分、长波SAR数据与高分、短波SAR数据从雷达坐标系地理编码到相同的WGS-84地心坐标系下，实现不同分辨率、波长SAR数据之间的匹配；

第四步，对低通形变分量Ⅰ进行空间插值及时间插值：

采用Kringing(克里金)与双线性插值法对低通形变分量Ⅰ进行空间插值；根据空间插值得到的形变信息采用二次多项式拟合时序形变曲线，通过时间插值最终得到高分、短波SAR成像时刻的低通形变分量Ⅱ；

第五步，获取高分、短波SAR成像几何空间下的低通形变分量Ⅲ：

将低通形变分量Ⅱ由低分、长波SAR数据的视线向转换到高分、短波SAR数据的视线向，得到高分、短波SAR数据成像几何空间下的低通形变分量Ⅲ；

第六步，获取高分、短波SAR数据的残余差分干涉相位：

首先将高分、短波SAR数据按照步骤一的过程生成差分干涉相位；然后减去步骤五得出的低通形变分量Ⅲ，得出高分、短波SAR数据的高分残余差分干涉相位；

第七步，提取高分、短波SAR数据的残余形变分量：

根据高分、短波SAR数据生成的高分残余差分干涉相位，利用DInSAR或时序InSAR方法得出高分残余形变分量，由于在上述步骤中大部分的形变相位已经去除，减缓了高分、短波SAR干涉相位中形变相位的梯度变化，因而可以保证相位解缠的准确性，提高了残余形变分量的解算精度；

第八步，获取地表总形变：

将解算后的高分残余形变分量与低通形变分量Ⅲ相加即可得到地表总的形变量。

上述步骤一具体为：首先获取覆盖同一地区的低分、长波SAR数据，若低分、长波SAR数据较少(比如少于12景)，则依时间顺序相邻两景低分、长波SAR数据依次组成干涉对；否则根据短基线干涉对选取原则，将所有满足条件的两景低分、长波SAR数据组成干涉对；然后利用该地区的DEM数据模拟该地区的地形相位，并从上述生成的每组干涉对中去除该相位，得到各个干涉对的差分干涉相位。

上述步骤二具体为：根据低分、长波SAR数据量，利用DInSAR或时序InSAR 方法从步骤一得到的差分干涉相位中提取地表形变相位；相对而言，由于采用的SAR数据空间分辨率低、波长较长，因此可以将该地表形变相位视为总体形变在时间与空间上的低通部分，即低通形变分量Ⅰ。

上述步骤三具体为：由于空间分辨率以及雷达波入射角等参数的不同，直接配准低分、长波SAR数据与高分、短波SAR数据是较为困难的。为此，通过DEM数据将两种SAR数据从雷达坐标系地理编码到相同的WGS-84地心坐标系下，实现不同SAR影像之间的配准；另外，在地理编码时，由于SAR卫星轨道的不精确性使得编码后的坐标可能存在一定程度的偏移，为此根据特定的地物目标对地理编码后的坐标进行校正。

上述步骤四具体为：由于两种SAR数据空间分辨率以及数据采样时间的不同，要将低分、长波SAR数据获取的低通形变分量Ⅰ融入到高分、短波SAR数据生成的差分干涉相位中，还需要进行空间与时间上的插值操作。空间插值为：首先采用Kringing插值法将低通形变分量规则格网化，然后采用双线性插值法获取高分、短波SAR数据像元处对应的形变分量；时间插值为：根据空间插值得到的相应形变分量，采用二次多项式拟合出每个像元的时序形变曲线，最后再插值得到高分、短波SAR数据采样时刻的形变信息，即低通形变分量Ⅱ。

上述步骤五具体为：通过已知的低分、长波SAR数据在几何框架下对地面目标的雷达波入射角和已知的高分、短波SAR数据在几何框架下对地面目标的雷达波入射角，可将低通形变分量Ⅱ由低分、长波SAR数据的视线向转换到高分、短波SAR数据的视线向，得到高分、短波SAR数据成像几何空间下的地表形变分量Ⅲ。

与现有技术相比，在高分、短波SAR数据较少的地区，本发明通过对不同分辨率的SAR数据进行融合，首先利用低分、长波SAR数据获取低通形变分量Ⅰ，进而经过空间插值与时间插值以及视向线的转换后得到高分、短波SAR数据成像几何空间下的低通形变分量Ⅲ；然后将该低通形变分量Ⅲ从高分、短波SAR的干涉相位中去除，得到高分残余干涉相位。由于大部分形变已经去除，因此可以利用高分残余干涉相位获取高精度的高分残余形变分量。最后，由低通形变分量Ⅲ与高分残余形变分量共同构成最终的地表形变量。因此本发明具有以下优点：1、由于联合不同分辨率、波长的SAR数据，可以克服单一SAR数据不足的缺陷，并且能够增强稳健性；2、由于可融合不同分辨率的SAR数据，因此够充分利用存档数据，有效降低对高分辨率SAR数据的数量要求，降低形变监测费用；3、只需利用少量高分辨率SAR数据，就可实现高精度、高空间分辨率的地表形变监测。

附图说明

图1是本发明的工作流程图；

图2是本发明的空间插值示意图，三维曲面由低通形变分量Ⅰ规则格网化得到，黑色点表示高分、短波SAR数据的像元；

图3是本发明中卫星的雷达波入射角示意图；

图4是本发明中去除低通形变分量Ⅲ前后的差分干涉对比图，其中(a)去除低通分量前的差分干涉相位图，(b)去除低通分量后的差分干涉相位图；

图5是采用本发明得到的地表形变与水准观测资料的对比图，其中(a)为仅使用11景高分TSX数据获取的形变速率与水准观测资料的对比，(b)为使用本发明，融合14景低分ASAR与11景高分TSX数据获取的形变速率与水准观测资料的对比。

具体实施方式

下面将对本发明作进一步说明。

实施例：实施例中的低分、长波SAR数据以Envisat ASAR(下述简称ASAR)为例，高分、短波SAR数据以TerraSAR-X(下述简称TSX)为例。

本发明的具体步骤为：

第一步：利用低分、长波ASAR数据，生成差分干涉相位

首先，获取覆盖同一地区的低分、长波ASAR数据。如果ASAR的数据量比较少的话(比如少于12景)，则依时间顺序相邻时间点的两景ASAR数据依次生成干涉对；否则，根据短基线干涉对选取原则，将所有满足条件的两景ASAR数据组成干涉对；然后利用DEM即数字高程模型(下述采用SRTM-DEM)数据模拟地形相位，并从每组干涉对中去除该相位，得到差分干涉图，差分干涉相位可记为：

上式中，上标S表示ASAR卫星，为地表变形引起的形变相位；为大气空间分布不均匀引起的大气相位；为轨道数据不精确造成的轨道误差相位；为DEM数据误差引起的地形误差相位；为失相干引起的噪声相位；

第二步：提取低通形变分量Ⅰ

根据ASAR的数据量，利用DInSAR或时序InSAR技术从差分干涉相位中提取地表形变相位相对于TSX数据而言，由于ASAR数据的空间分辨率低、波长较长，因此可将视为总体形变φ_Def在时间与空间上的低通部分，即低通形变分量Ⅰ。

第三步：配准ASAR与TSX两种影像

由于空间分辨率以及雷达波入射角等成像参数的不同，直接配准ASAR与TSX是较为困难的。为此，可以通过SRTM-DEM数据将ASAR与TSX从雷达坐标系地理编码到相同的WGS-84地心坐标系下，实现不同SAR数据之间的配准。另外，在地理编码时，由于SAR卫星轨道的不精确性使得编码后的坐标可能存在一定程度的偏移。为此，可以根据一些特定的地物目标(如道路交叉口、独立地物等)对地理编码后的坐标进行适当校正。

第四步：对低通形变分量Ⅰ进行空间插值及时间插值

可以分为2步：两者的PS点空间上不完全匹配，采样时间是不完全一致。

因为空间分辨率以及数据采样时间的不同，要将ASAR获取的低通形变分量Ⅰ融入到TSX数据生成的差分干涉相位中，还需要进行空间与时间上的插值操作。空间域插值方法为：首先采用Kringing(克里金)插值法将ASAR数据获取的低通形变Ⅰ规则格网化，然后利用双线性插值法获取TSX像元处的相应形变信息。时间域插值方法为：根据空间插值得到的相应形变信息，采用二次多项式拟合出每个像元的时序形变曲线，最后再插值得到TSX数据采样时刻的形变信息即低通形变分量Ⅱ：

第五步：转换低通形变分量Ⅱ的视线向

因为不同卫星采用的雷达波入射角不同，需要利用雷达波入射角，将利用低分、长波ASAR提取的低通形变分量Ⅱ(即)转换到高分、短波TSX数据所在的雷达视线向即低通形变分量Ⅲ；

上式中，表示低通形变分量转换到TSX框架下的形变相位。θ_S与θ_T分别表示地面目标在ASAR与TSX几何框架下的雷达波入射角。

第六步：获取高分残余差分干涉相位

首先利用高分、短波TSX数据生成差分干涉相位图，差分干涉相位可记为：

然后，将低通形变分量Ⅲ(即)从上式中去除，得到高分残余干涉相位

第七步：提取高分残余形变分量

根据高分残余干涉相位利用DInSAR或时序InSAR技术从公式(5)中获取地表残余形变分量因为上述步骤后，公式(4)中的大部分形变相位已经去除，减缓了高分TSX数据中形变相位的梯度变化，因而可以根据公式(5)有效地进行相位解缠，保证了残余形变分量的解算精度；

第八步：获取地表总形变

利用高分残余干涉相位解算出残余形变分量后，与低通形变分量Ⅲ相加即可得到地表总的形变相位φ_Def：

。