-
时序合成孔径雷达InSAR干涉测量技术主要通过长时间序列的SAR数据进行分析,以去除(或削弱)干涉相位中的轨道、大气、数字高程模型(digital elevation model,DEM)误差以及低相干性等因素的影响,从而提取更精准的地表形变。自20世纪90年代中后期起,经过近十余年的发展,时序InSAR 技术获得很大的发展,并产生了众多分支,主要有相位叠加技术[1]、永久散射点(persistent scatterers or permanent scatterers,PS)InSAR[2]、最小二乘方法(least square,LS)[3]、小基线集方法(small baseline subset algorithm,SBAS)[4]、相干目标方法(coherence target/coherence point target,CT)[5]和临时性相干点(temporal coherence point,TCP)InSAR[6]等。时序InSAR技术的发展显著提高了地表形变测量的精度和可靠性[7, 8],也使InSAR技术与固体地球物理学、冰川学、水文学等学科交叉融合逐渐深入,在地球科学研究中发挥的作用越来越大。
在时序InSAR技术发展中,通过建立InSAR观测值的函数模型进行参数的估计,虽然可以降低大气、轨道、DEM等误差对观测结果的影响,但是由于误差影响因素的复杂性和不确定性,时序InSAR技术的随机模型发展一直以来都很缓慢,难以对InSAR观测值的参数精度进行有效评估。在已有的时序InSAR技术中,未知参数(主要是视线向LOS(line of sight)位移)的精度估计主要依赖于稳定区域的均方根、误差的相干性或两者的联合[9-12],或者在形变参数估计过程中不考虑各观测值的精度差异,而假定各观测值的精度相当并进行等权参数估计。在对时序InSAR技术的随机模型研究中,Kampes等[13]首先利用了方差分量估计(variance component estimation,VCE)对PSInSAR技术中的干涉图质量进行评估,来估计参数的精度;González等[7]基于InSAR干涉图中误差的时空相关特性,建立了Gauss-Markov数学模型评估SBAS方法的观测误差,并建立观测值的协方差矩阵;Zhang[6]基于TCP模型方法,融合VCE方法来建立随机模型评估待求参数的精度。已有的时序InSAR技术中采用先验权或等权方法进行参数估计,虽然已经取得了许多成功的应用,也是迄今时序InSAR技术中较为通用的处理方法,但是这种定权方法是否合理,对结果影响如何,一直以来没有定论,使得时序InSAR技术的精度和可靠性受到一定质疑。时序InSAR误差模型的构建及精度评定,对时序InSAR的数学模型完善有着重要的作用,研究如何建立一个适合的随机模型进行参数的估计和精度评估有着重要的意义。
本文结合已有的时序InSAR函数模型,拓展了Kampes[15]的随机函数模型,建立了时序InSAR的一般数学模型,完善了时序InSAR的数学理论体系。利用模拟数据,分别基于单主影像和多主影像时序InSAR的加权和等权参数解算精度进行了模拟试验。
-
测量平差问题的数学模型包括函数模型和随机模型,其主要作用是依据函数模型中给出的观测值与未知量之间的函数关系,顾及观测量的先验方差和协方差,确定观测值的协因数阵或权阵,按最小二乘原理作出未知量的最佳估值。
函数模型是描述观测量与未知量的数学函数关系的模型,是确定客观实际的本质或特征的模型。函数模型不同,其未知参数的估计平差方法相应不同。函数模型可分为线性和非线性模型。对于InSAR观测值与待求参数的关系,根据Gauss-Markoff 模型[14]表示如下:
(1) (2) 其中,φ表示干涉相位观测值;A是设计矩阵;Cφ是方差协方差阵;Qφ是协因数阵;σ2是单位权先验方差因子(σ02=1)。
根据InSAR相位φ的组成部分,具体函数模型可以写为:
(3) 其中,φtopo是DEM高程误差导致的相位变化,可以表示为转换系数β和高程误差Δh的乘积;φdefo是形变相位,可以表示为速率v和时间t的乘积;λ为SAR传感器波长;ε表示不相关随机噪声,E{ε}=0;待估参数为Δh和v。
由于干涉图是基于SAR影像的组合,对干涉图的误差研究可以从单视复型(single look complex,SLC)影像的基础上进行分析。假设在第k幅SLC影像中,原始相位观测量(在H点处)矢量为:
(4) 其中,k=0,…,K。相位观测值的方差为[15]:
(5) 其中,Qnoisek表示热噪声;Qatmok为大气的方差协方差阵;la,b(a和b为一幅干涉图中按顺序排列的像素点,分别可取0,1,2,…,H)为像素点a和b的实际空间距离;大气相位特征利用协方差函数描述[7, 16]。
对于所有的SLC影像,其相位观测值为(K+1)H×1矢量:
(6) 可推导N幅干涉图的相位:
(7) 其中,Λ为系数阵。
对于时序InSAR处理,作为一种相对空间测量技术,在同一影像中存在参考基准点问题,为简便可以假定第一个点为参考点,对干涉图中选定的时序数据集与参考点作差,可以得到平差网观测方程为:
(8) 其中,Ω为系数阵。
为了解算方便,在实际数据处理时对式(8)按干涉图顺序排列乘以一个变换矩阵P,得到按观测点顺序排列的矩阵:
(9) 根据方差传播定律可得到干涉图中相位差分量的方差-协方差阵[15]:
(10) 从而可得出式(3)中的待估参数的方差-协方差阵如下[15]:
(11) 其中,
,待估参数的最小二乘解为[15]:(12) -
基于单主影像的时序InSAR技术,是指在干涉过程中,其中一景SAR影像为主影像,其余影像为辅影像的时序InSAR技术。代表算法为PSInSAR技术[2, 17],这种算法从干涉图中估计形变参数,不用考虑临界基线的限制。假设有(K+1)幅SLC影像,选定k=0的SAR影像为主影像,则可形成K幅干涉图,式(7)中的系数阵Λ可表示为:
(13) 将式(13)代入式(9)即可以对干涉相位的随机模型进行估计。
-
基于多主影像的时序InSAR技术,主要是选取干涉像对中干涉质量高的像对进行处理。代表算法有SBAS[4]、CT[3]和TCP-InSAR[6]。多主影像算法为了选取高相干性的像对,一般选取时空基线短的像对进行处理。设(K+1)幅SLC影像,形成了N幅干涉图,其中
,式(7)中的系数阵Λ可表示为:(14) 将式(14)代入式(9)即可以对干涉相位的随机模型进行估计。
-
迄今为止,对时序InSAR的随机模型研究虽然有一定的时展,但是在实际处理中应用并不多。在实际过程中对时序InSAR利用等权估计还是随机模型定权,值得进一步探讨,下面分别就等权与随机模型定权的时序InSAR 参数估计方法进行模拟研究。
本文以ERS和ENVISAT卫星参数为基础,模拟具有相似时空基线分布的SAR影像,比较等权以及利用随机模型定权条件下的参数估计。以没有误差的二维匀速位移场开始,利用Mogi模型,模拟区域大小为10 km×10 km,深度为7.5 km(X为6.5 km、 Y为6.5 km),体积变化速率为-0.25×10-3km3/a的形变。对每景SAR影像加入随机噪声,噪声均值为15°,标准差为10°,大气噪声范围为-2~2 cm,DEM误差范围为-10~10 m。
-
利用PSInSAR技术进行单主影像的数据模拟研究,随机得到了1992~2001年间30个历元的降轨SAR数据,其中PS点数量为1 000个。数据的时空基线分布如图 1(a) 所示,选取1996-01-17的影像为主影像,利用30景数据可组成29景模拟的干涉图。根据PSInSAR方法的原理[2],对PS点组成平差网,分别利用等权方法(即不定权)和随机模型定权方法进行参数估计,得到估计的DEM误差以及形变速率结果如图 2。
图 2 单主影像时序InSAR等权(a~f)与随机模型加权(g~l)参数估计模拟
Figure 2. Parameters Estimated by Equal-Weight Model and Weighted Model for Single-Master Image TS-InSAR
等权条件下的DEM及速率估计误差均方差分别为0.3 m、0.16 mm/a;加权条件下的DEM及速率估计误差均方差分别为0.29 m、0.15 mm/a。比较等权与加权模型下的单主影像参数估计误差结果,加权条件下的均方差较小,结果有所改善,但是计算效率较等权条件低。
-
利用TCP方法进行多主影像时序InSAR的数据模拟研究,随机得到了1999~2006年间30个历元的降轨SAR数据,其中高相干点数量为1 000个。利用30景数据形变干涉像对,选取垂直基线200 m的短基线像对,共获得58景模拟的干涉图,其时空基线分布如图 3。根据TCP方法[6],对高相干点组成三解网,利用等权方法(即不定权)选取可靠的弧段进行平差解算,得到估计的DEM误差以及形变速率如图 4。
图 3 多主影像时序InSAR模拟干涉像对时空基线分布
Figure 3. Distribution of Temporal Spatial Baseline for Multi-master Image TS-InSAR
图 4 多主影像时序InSAR等权(a)~(f)与随机模型加权(g)~(l)参数估计模拟
Figure 4. Parameters Estimated by Equal-Weight Model and Weighted Model for Multi-master Image TS-InSAR
等权条件下的DEM及速率估计误差均方差分别为0.7 m、0.23 mm/a。而按给出的随机模型对观测量进行加权,再进行参数估计,得到估计的DEM误差以及形变速率如表 1(测试用电脑配置中的CPU为4核16个CPU,主频为2.12 GHz,内存为32 GB)。加权条件下的DEM及速率估计误差均方差分别为0.7 m、0.22 mm/a。比较等权与加权模型下的多主影像参数估计误差结果,加权条件下的均方差较小,结果有所改善,但计算效率较等权条件低。
表 1 等权与加权模型的参数估计误差
Table 1. Parameters Estimation Error by Equal-Weight Model and Weighted Model
最大值 最小值 平均值 均方差 计算耗时/s 单主影像 等权 DEM/m -1.1 0.5 -0.2 0.3 213 形变速率/(mm\5a-1) 0.55 -0.36 0.1 0.16 加权 DEM/m -0.83 0.80 0.02 0.29 731 形变速率/(mm\5a-1) 0.54 -0.33 0.14 0.15 多主影像 等权 DEM/m -3.5 1.5 -0.7 0.7 329 形变速率/(mm\5a-1) -1.1 0.54 -0.17 0.23 加权 DEM/m -3.5 1.8 -0.6 0.7 1 014 形变速率/(mm\5a-1) -1.07 0.58 -0.12 0.22 -
本文基于时序InSAR的一般函数模型,分别对单主影像和多主影像时序InSAR方法的随机模型进行推导,同时利用模拟数据比较分析了在等权与加权条件下时序InSAR参数估计精度的差异,有以下结论:
(1) 通过对时序InSAR误差模型的推导,将PS InSAR的随机模型扩展到时序InSAR的一般随机模型,使得时序InSAR的数学模型更为完善,丰富了时序InSAR技术的理论体系。
(2) 本文利用模拟数据,比较分析等权与加权条件下参数估计的精度,结果表明,加权条件下参数估计精度有一定的改善,但是计算效率较低。
干涉影像的时序分析是基于像素点的误差分析,若考虑像素与像素之间的相关性,同时考虑同一影像中不同像素点大气误差大小的差异,计算时的内存空间要求成几何倍数增加。目前,在时序InSAR数据处理时,都是基于像素点进行分析,同时每幅影像估计一个误差作为影像中各像素点的误差。在这种条件下,进行时序InSAR的随机模型估计,虽然在理论上更为合理,但是计算复杂度增加,对观测精度改善不显著,而利用等权进行时序InSAR的参数估计能节省计算效率,达到较高的精度,所以利用等权模型进行参数估计是可行的。本研究主要侧重于模型的构建和算法的推导,并以模拟数据进行实验验证。并且,目前实际SAR数据的选择和处理在模型精度评价方面还有一定的困难,这将为下一步的研究方向。
-
摘要: 基于时序InSAR函数模型,分别建立了单主影像和多主影像时序InSAR误差模型,在理论上完善了时序InSAR数学模型。利用构建的随机误差模型,模拟试验研究了随机误差模型对时序InSAR待估参数精度的影响。结果表明,与等权模型相比较,加权模型(随机误差模型)估计的参数精度有一定提高;但由于加权条件下的参数估计模型复杂、计算效率低,目前利用等权方法进行时序InSAR的参数估计更简便易行。Abstract: Based on the function model of TS-InSAR (time series interferometric synthetic aperture radar) technique, this paper establishes stochastic function models (weighted/equal-weight) for single-master and multi-master image TS-InSAR method respectively to enrich the mathematical model of TS-InSAR. A simulation test is used to estimate the precision of model parameters for our equal-weighted function model (stochastic model). Compared to the equal-weight model, the precision of deformation parameters shows little improvement over the weighted model. However, complications in the weighted model means that much disk-space is consumed with low computational efficiency. Most computers therefore cannot undertake TS-InSAR analysis tasks with reasonable hardware configuration. At present, the equal-weight model is feasible for TS-InSAR.
-
Key words:
- TS-InSAR /
- stochastic model /
- parameter estimation
-
表 1 等权与加权模型的参数估计误差
Table 1. Parameters Estimation Error by Equal-Weight Model and Weighted Model
最大值 最小值 平均值 均方差 计算耗时/s 单主影像 等权 DEM/m -1.1 0.5 -0.2 0.3 213 形变速率/(mm\5a-1) 0.55 -0.36 0.1 0.16 加权 DEM/m -0.83 0.80 0.02 0.29 731 形变速率/(mm\5a-1) 0.54 -0.33 0.14 0.15 多主影像 等权 DEM/m -3.5 1.5 -0.7 0.7 329 形变速率/(mm\5a-1) -1.1 0.54 -0.17 0.23 加权 DEM/m -3.5 1.8 -0.6 0.7 1 014 形变速率/(mm\5a-1) -1.07 0.58 -0.12 0.22 -
[1] Sandwell D T, Sichoix L. Topographic Phase Recovery From Stacked ERS Interferometry and a Low-Resolution Digital Elevation Model[J]. Journal of Geophysical Research:Solid Earth (1978-2012), 2000, 105(B12):28211-28222 [2] Ferretti A, Prati C, Rocca F. Nonlinear Subsidence Rate Estimation Using Permanent Scatterers in Differential SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(5):2202-2212 [3] Usai S. A New Approach for Long Term Monitoring of Deformations by Differential SAR Interferometry[D]. Netherlands:Delft University of Technology, 2001 [4] Berardino P, Fornaro G, Fusco A, et al. A New Approach for Analyzing the Temporal Evolution of Earth Surface Deformations Based on the Combination of DIFSAR Interferograms[C]. IGARSS, Sydney, 2001 [5] Blanco-Sanchez P, Mallorqui J J, Duque S, et al. The Coherent Pixels Technique (CPT):An Advanced Dinsar Technique for Nonlinear Deformation Monitoring[J]. Pure and Applied Geophysics, 2008, 165(6):1167-1193 [6] Zhang L, Ding X, Lu Z. Ground Settlement Monitoring Based on Temporarily Coherent Points Between Two SAR Acquisitions[J]. Journal of Photogrammetry and Remote Sensing, 2011, 66(1):146-152 [7] GonzaLez P J, FernaNdez J. Error Estimation in Multitemporal InSAR Deformation Time Series, with Application to Lanzarote, Canary Islands[J]. Journal of Geophysical Research:Solid Earth (1978-2012), 2011, 116, DOI: 10.1029/2011JB008412 [8] 何平,温扬茂,许才军,等.多时相InSAR技术的廊坊地区地下水体积变化研究[J].武汉大学学报·信息科学版,2012,37(10):1181-1185 He Ping, Wen Yangmao, Xu Caijun, et al. Volume Change of Groundwater in Langfang Region Derived from Multi-temporal InSAR[J]. Geomatics and Information Science of Wuhan University, 2012, 37(10):1181-1185 [9] Hooper A, Zebker H, Segall P, et al. A New Method for Measuring Deformation on Volcanoes and Other Natural Terrains Using InSAR Persistent Scatterers[J]. Geophysical Research Letters, 2004, 31(23):275-295 [10] Kwoun O I, Lu Z, Neal C, et al. Quiescent Deformation of the Aniakchak Caldera, Alaska, Mapped by InSAR[J]. Geology, 2006, 34(1):5-8 [11] Anderssohn J, Motagh M, Walter T R, et al. Surface Deformation Time Series and Source Modeling for a Volcanic Complex System Based on Satellite Wide Swath and Image Mode Interferometry:the Lazufre System, Central Andes[J]. Remote Sensing of Environment, 2009, 113(10):2062-2075 [12] Ketelaar V B H. Satellite Radar Interferometry:Subsidence Monitoring Techniques[M]. Netherlands:Springer, 2009 [13] Kampes B M, Hanssen R F. Ambiguity Resolution for Permanent Scatterers Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(11):2446-2453 [14] Hanssen R F. Radar Interferometry:Data Interpretation and Error Analysis[M]. New York:Springer, 2001 [15] Kampes B M. Radar Interferometry:Persistent Scatterer Technique[M]. New York:Springer, 2006 [16] Sudhaus H, Sigurjon J. Improved Source Modelling Through Combined Use of InSAR and GPS Under Consideration of Correlated Data Errors:Application to the June 2000 Kleifarvatn Earthquake, Iceland[J]. Geophysical Journal International, 2009, 176(2):389-404 [17] 许才军,何平,温扬茂.利用PSInSAR研究意大利Etna火山的地表形变[J].武汉大学学报·信息科学版,2011,36(9):1012-1016 Xu Caijun, He Ping, Wen Yangmao. Surface Deformation of Mt. Etna, Italy from PSInSAR[J]. Geomatics and Information Science of Wuhan University, 2011, 36(9):1012-1016 -