一种分布式星载InSAR高程反演克拉美罗界的确定方法

一种分布式星载insar高程反演克拉美罗界的确定方法
技术领域
1.本发明属于合成孔径雷达技术领域，尤其涉及一种分布式星载insar高程反演克拉美罗界的确定方法。


背景技术：

2.数字高程模型(dem)作为重要的地理信息产品，在地质学、气象学等领域发挥了重要作用。分布式星载合成孔径雷达(sar)系统作为一种较新的星载sar体制，具有灵活多变的系统构型和高时空分辨率数据获取能力，因此在实现高精度高分辨率的高程反演方面具有优势。
3.在利用分布式星载sar进行高程反演过程中，大气会对高程反演精度造成严重影响，因此考虑大气特征在内的精度评估方法非常重要。克拉美罗界(crb)是评估参数性能的重要指标，可以用来衡量无偏估计器的性能。文献1(guarnieri a m,tebaldini s.hybridbounds for crustal displacement field estimators in sar interferometry[j/ol].ieee signal processing letters,2007,14(12):1012-1015.doi:10.1109/lsp.2007.904705.)给出了利用混合克拉美罗界(hcrb)评估视线向地表形变的方法，该方法将大气相位屏作为随机参数考虑在内，对性能进行可行的评估；文献2(prats-iraola p,lopez-dekker p,de zan f,等.performance of 3-d surface deformation estimation for simultaneous squinted sar acquisitions[j/ol].ieee transactions on geoscience and remote sensing,2018,56(4):2147-2158.doi:10.1109/tgrs.2017.2776140.)在文献1的基础上做了进一步的研究，利用hcrb并对公式进行了扩展，最终获得三维形变性能。就我们目前所知，现有的考虑大气影响的crb评估方法都将大气相位作为随机变量考虑，并假设其具有高斯特性。
[0004]
然而，事实上，大气相位成分复杂。文献3(mulder g,van leijen f j,barkmeijer j,等.estimating single-epoch integrated atmospheric refractivity from insar for assimilation in numerical weather models[j/ol].ieee transactions on geoscience and remote sensing,2022:1-1.doi:10.1109/tgrs.2022.3177041.)证明了大气并不具有高斯特性，因此高斯假设在crb计算中并不成立。另外，由于气象环境的多变性，风、雨、雪等不同的气象条件都会导致不同的大气统计特性；在空间和时间上，大气相位也并非一个平稳过程。因此，为了对分布式星载insar高程反演的无偏估计器进行评估，考虑大气非平稳非高斯特征的crb计算方法是十分重要的。


技术实现要素：

[0005]
为解决上述问题，本发明提供一种分布式星载insar高程反演克拉美罗界的确定方法，该方法不再对整个干涉图中的大气相位进行高斯和平稳的假设，在空间上采用后验crb的思路，通过高斯混合模型拟合误差分布来获得最终的crb结果。
[0006]
本发明的技术解决方案是：
[0007]
一种分布式星载insar高程反演克拉美罗界的确定方法，该方法的步骤包括：
[0008]
步骤1，获取分布式星载insar观测地区的参考数字高程模型h，并根据获得的参考数字高程模型h得到待估计数字高程模型
[0009]
所述的观测地区的参考数字高程模型h能够从航天飞机雷达地形测绘(srtm)数据中获得；
[0010]
所述的根据获得的参考数字高程模型h得到待估计数字高程模型的方法为：
[0011]
设参考数字高程模型h中第m行n列的像元为[h]
m,n
＝h
m,n
，待估计数字高程模型中第m行n列的像元为则
[0012][0013]
其中，
△hm,n
＝h
m,n+1-h
m,n
；∈为介于0到1的参数，v
m,n
为和h的误差所导致的过程噪声，过程噪声为具有零均值和方差的高斯噪声模型，零均值和方差都和h的精度相关；
[0014]
步骤2，计算分布式星载insar的波数，具体方法为：
[0015]
获取由分布式星载insar得到的步骤1中所述的观测地区的干涉图其中干涉图总数量记为n；第k幅干涉图中，第m行第n列的像元表示为其中k＝1,2,
…
,n；则第m行第n列所有像元处获得的干涉相位组成向量
[0016]
记每幅干涉图对应的雷达波长为λk，首先计算每幅干涉图对应的雷达间的干涉基线，记第k幅干涉图xk对应的干涉基线为bk；再计算每幅干涉图对应的主轨道孔径中心位置到场景中心的斜距，记第k幅干涉图xk对应斜距为rk；最后计算每幅干涉图对应的入射角，第k幅干涉图xk的入射角记为θk；根据以上参数，计算得到干涉图xk对应的波数ξk：
[0017][0018]
步骤3，根据步骤1得到的参考数字高程模型h和步骤2得到的波数ξk计算第k幅干涉图中的第m行第n列像元处的相位差
[0019][0020]
步骤4，对步骤3得到的相位差进行误差相位处理得到符合均值为方差为的高斯分布；具体为：
[0021]
对于干涉图中的第k幅图中的第m行第n列像元，选取它周围w
×
w大小窗口内的相位差干涉图(w为奇数)，记为则有：
[0022][0023]
利用高斯混合模型对进行拟合，得到符合的分布，为使后续计算期望更加便捷，根据中心极限定理，同时由于w
×
w的窗口内像元数量相对于整个干涉图的像元总数很小，假设窗口内大气信号平稳、误差相位近似符合高斯分布，即得到符合的分布为高斯分布，同时，由于对每个像元独立进行窗口内高斯分布的假设，因此也可以表示整幅干涉图大气相位的非平稳特性；
[0024]
设经过拟合后，对于第m行第n列的像元，其的分布符合均值为方差为的高斯分布；
[0025]
步骤5，从待估计数字高程模型的每行第一列开始，逐像元计算crb，具体计算方法为：
[0026]
记待估计数字高程模型的fisher信息量为j,设第m行第n列的fisher信息量为[j]
m,n
＝j
m,n
，在crb计算过程中，将待估参数分解为则相应的fisher信息量定义为：
[0027][0028]
对于第m行第n列，有：
[0029][0030][0031][0032][0033]
递推得到第m行第n+1列的fisher信息量：
[0034][0035]
对fisher信息量求倒数，得到各个像元处的crb：
[0036]
interferometry data interpretation and error analysis.2001.)；根据以上参数，计算得到干涉图xk对应的波数：
[0052][0053]
步骤3，根据实际获得干涉图和计算的地形相位计算相位差。将第k幅干涉图中的第m行第n列像元处的相位差记为则处的相位差可以表示为
[0054]
步骤4，处理干涉图的误差相位。对于第k幅图中的第m行第n列像元，选取它周围w
×
w大小窗口内的相位差(w为奇数)，记为则有：
[0055][0056]
利用高斯混合模型对进行拟合，得到符合的分布。为使后续计算期望更加便捷，根据中心极限定理，同时由于w
×
w的窗口内像元数量相对于整个干涉图像元总数很小，可假设窗口内大气信号平稳、误差相位近似符合高斯分布。同时，由于对每个像元独立进行窗口内高斯分布的假设，因此也可以表示整幅干涉图大气相位的非平稳特性。
[0057]
进行拟合后，对于第m行第n列的像元，其的分布符合均值为方差为的高斯分布。
[0058]
步骤5，从每行第一列开始，逐像元计算crb。记整幅图的fisher信息量为j,设第m行第n列的fisher信息量为[j]
m,n
＝j
m,n
。在crb计算过程中，将待估参数分解为则相应的fisher信息量可以定义为：
[0059][0060]
对于第m行第n列，有：
[0061][0062]
[0063][0064][0065]
可以递推地得到第m行第n+1列的fisher信息量：
[0066][0067]
对fisher信息量求倒数、得到各个像元处的crb：
[0068][0069]
若比较整幅图的高程反演精度，根据均方误差(mse)的定义，可以对比mse与整幅图的平均crb：
[0070][0071]
实施例
[0072]
接下来结合具体参数给出实施的例子。
[0073]
在本实例中，我们考虑一个3星同轨道面分布式星载sar系统。其卫星轨道根数和平台参数如表1所示。假设三颗星构成的系统在6天时间内分别进行5次重轨干涉，共获得15幅干涉图，干涉图大小为182像素
×
182像素，像素间隔为1m。高程反演的场景中心位于东经69
°
、北纬30
°
，场景选自30m分辨率的先进星载热发射和反射辐射仪全球数字高程模型(aster gdem)，并按比例降低高度。仿真场景的图像见图2。
[0074]
表1
[0075][0076]
首先按照步骤1，获取观测地区的30m分辨率dem作为真实高程，等比例缩小至最高高度降为18.3m。参考dem由真实高程经过窗口平均后得到。
[0077]
执行步骤2，计算每幅干涉图对应的波数。
[0078]
执行步骤3，计算实际获得的干涉图和计算得到的地形相位的相位差。
[0079]
执行步骤4，选取窗口大小w＝21，则对于每个像元位置，利用它周围21
×
21窗口内的相位差，拟合得到每个位置处像元对应的方差和均值。
[0080]
执行步骤5，计算高程反演的最终克拉美罗界。
[0081]
图2给出了15幅干涉图情况下两种平稳性的高程反演结果和真实高程，两者mse均接近crb，但平稳大气条件下mse更小，精度更高。图3展示了平稳和非平稳的大气条件下，利用不同数量的干涉图，本方法计算出的crb和高程反演得到的mse。可以看出本方法计算得到的crb随着干涉图数量上升而降低；非平稳特性下的crb均要高于平稳大气条件下的结果；且实验结果证明两种情况下高程反演得到的mse均接近crb，从而证明了crb的正确性。
[0082]
当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当可根据本发明做出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。