徐进军,王振华*,邢 诚,赖慧斌
(1.武汉大学 测绘学院,湖北 武汉430079)
GB-SAR 是一种基于雷达的地面遥感成像系统,具有高空间分辨率、高采样频率以及零空间基线的特点[1],已被广泛应用于桥梁[2]、滑坡[3-4]、大坝[5]、冰川[6]、火山[7]、建筑物[8-9]等变形监测中。研究表明,GB-SAR 变形监测成果易受大气效应影响,因此需要选择适当的方法将其剔除,以提高变形监测精度。
气象数据改正法是常见的GB-SAR 大气改正方法之一,通过测量现场的气象数据(大气压强、温度、相对湿度)建立大气扰动模型来进行大气改正,在小范围区域内精度可达mm 级[10];但该方法受制于气象数据的观测精度,当采集的气象数据精度较低或气象条件变化剧烈时,其大气改正结果不理想,精度有限。针对该情况,本文利用稳定点和大气监测数据进行了相关的研究。
通过对GB-SAR 影像的干涉处理可得到目标点的干涉相位φw,其中包括4 个部分:
式中,φdis为目标点形变引起的相位;φatm为大气变化产生的大气误差相位;φnoise为相位测量噪声;n为相位整周模糊度。
φnoise可通过滤波进行削弱,削弱后的φnoise相对于φdis和φatm可忽略不计,因此φdis可被表示为:
雷达视线向的位移量可由形变相位得到,表达式为:
式中,λ为雷达波的波长。
假设向雷达视线向距离为r的目标发射频率为f的电磁波,并于t时刻返回,在小区域范围内,其受大气影响的回波相位可表示为:
式中,c 为真空中的光速;n(t)为与时间相关的折射指数。
由于气象条件随时间的变化而变化,因此大气折射指数在不同时刻也不相同,同一目标点两个时刻t1和t2之间的回波相位差为:
为了便于计算,一般用折射率N表示折射指数n,其关系为N=(n-1)×106,则式(5)可写为:
折射率N可由大气压强P(单位为hPa)、温度T(单位为K)、水汽压e(单位为hPa)[11]表示:
水汽压e则可由相对湿度RH(单位为%)和温度T表示[12]:
将式(8)代入式(7),得到由t时刻大气压强、温度、相对湿度计算的折射率,即
联合式(6)、式(9),可得到大气误差相位变化为:
实验以意大利IDS 公司生产的GB-SAR 设备为平台,采用IBIS-L 静态测量模式,对三峡库区链子崖危岩体进行监测,测量现场与雷达位置如图1 所示。设备的工作频率为ku 波段,距离向分辨率为0.5 m,方位向分辨率为4.4 mrad。测量时间为2018 年3 月10 日17:05-3 月11 日1:15,采样间隔约为5 min,共采集了94 幅影像。
图1 监测场景
监测现场的气象数据由PH-II 手持式气象仪采集,其技术参数见表1。一共放置了两台气象仪,一台在雷达附近,一台在链子崖变形体上;两台气象仪的观测序列取平均作为现场的气象数据,如图2 所示。经计算可知,相对湿度的平均值为62.0%,温度的平均值为16.15℃,大气压强的平均值为973.6 hPa。
表1 PH-II 手持气象仪的技术参数
采集的雷达影像数据经过IBIS Data Viewer 软件初步处理,可得到监测区域的热信噪比图(图3a)和相位稳定性图(图3b),可以看出,监测区域的热信噪比值大多超过30 dB,相位稳定性值大多超过2.5,因此以这两个数值作为目标点选取的阈值。在图1 的变形区域中选取3 个像元点(P1~P3)作为变形点,从而分析气象数据大气改正的效果。3 个变形点的位置分布如图3c 所示,由于观测期间短,整个区域没有变形,因此所选变形点的变形值理论上均为0 mm。
图2 监测现场气象数据
图3 雷达影像参数图和变形点位置分布图
变形点P1~P3沿雷达视线向未经大气改正的监测结果如图4 所示,变化约在0 ~15 mm 之间。利用采集的气象数据对其进行大气改正,改正结果如图5 所示,3 个点的变化约在-1 ~6 mm 之间。对比图4、5可知,利用原始气象数据进行大气改正有一定的效果,但与实际情况有较大的偏差。由于气象数据存在观测误差,上述结果说明观测的气象数据与实际气象值有一定的偏差,需对其进行修正。
图4 P1~P3 雷达视线向未经大气改正的监测结果
图5 P1~P3 雷达视线向通过原始气象数据大气改正的结果
气象数据存在观测误差,那么由其计算的折射率也存在误差,结合式(9),根据误差传播定律,折射率的标准差σN可表示为:
为了书写简便,将式(9)中第二项的数字用符号α代替,即令α=3.73×105×6.112 1,则式(11)中各项偏导数可表示为:
结合式(10),大气误差相位的标准差σφatm也可表示为:
大气误差相位误差对雷达视线向位移量的误差σΔLOS可表示为:
在本文实验中,雷达视线向的平均观测距离约为800 m,结合气象仪的观测精度和各气象因子观测值的平均值,根据式(11)~(16),计算大气压强、温度、相对湿度各自变化所引起的雷达视线向变化,结果如表2 所示。
表2 单个气象因子变化引起的雷达视线向变化
由表2 可知,相对湿度是影响雷达视线向形变测量精度的主要因素,温度和大气压强的影响相对其可以忽略,因此认为图5 中的形变偏差主要源于现场相对湿度测量的不准确,只对现场测量的相对湿度进行修正即可。
为了消除现场相对湿度测量不准确的影响,本文提出采用稳定点(没有任何变形的点)修正现场气象测量值的方法。具体过程为:
1)在小区域范围内认为气象变化是均匀的,湿度修正模型选择线性模型,可表示为:
式中,RH为气象观测站所测得的相对湿度;为修正后的相对湿度;a1、a0为修正系数。
2)在测区选择稳定点。稳定点没有变形,其干涉相位的变化主要来自大气的影响,利用式(10)建立稳定点实测相位值与修正系数的计算模型。为能求得有效的修正系数,至少需要一个稳定点两个时刻以上的相位观测值。
3)利用式(17)修正气象数据,再代入式(9)和式(10)进行大气改正。
按照§2.2 中选择变形点的阈值选择稳定点,为了对比稳定点数量对大气改正结果的影响,本文选择了3 个稳定点(GSP1 ~GSP3),其位置分布如图6所示。
图6 稳定点位置分布
1)相对湿度整体修正。假设观测期间式(17)中的修正系数不变,即所有时刻的相对湿度观测值都用相同的修正系数进行修正。本文分别通过一个稳定点GSP1 和3 个稳定点GSP1 ~GSP3 对所观测的相对湿度序列作整体修正,再利用修正后的气象数据对P1~P3监测结果进行大气改正,结果如图7 所示。
由图7 可知,将相对湿度整体修正后再进行大气改正,相较于原始气象数据大气改正,其结果有了很大提高,基本都在±2 mm 以内;但其在40 min 和340 min左右都出现了极值点。观察所采集的相对湿度数据发现,整体的湿度序列变化趋势是不一样的,极值点出现在湿度序列趋势陡变的地方,鉴于此,本文将相对湿度序列按照其变化趋势来分时段修正。
2)相对湿度分时段修正。相对湿度变化趋势不同的时段的修正系数是不一样的,结合实验情况,本文将其分为4 个时段来修正,分段情况如图8 所示。仍按上述选择稳定点的情况对湿度进行分时段修正,修正后气象数据大气改正结果如图9 所示。
通过对比图7、9 可知,分时段修正气象数据的大气改正结果更好,波动主要位于±1 mm 之间,并消除了之前出现的极值点情况,符合实际变形。由此可得以下结论:①通过稳定点能实现对气象值的正确修正,所采用的气象改正模型是正确的;②结合相对湿度变化趋势对其进行分时段修正后再进行大气改正,可得到更好的结果;③选择3 个稳定点的大气改正结果稍优于选择一个稳定点的大气改正,但区别不大,实际工作中稳定点选择比较灵活。
图7 相对湿度整体修正后进行大气改正的结果
图8 相对湿度分时段情况
图9 相对湿度分时段修正后进行大气改正的结果
本文探讨了GB-SAR 气象数据大气扰动模型以及气象数据观测精度对其变形监测精度的影响,并在气象数据大气改正方法的基础上,提出了一种基于稳定点分时段修正观测气象数据的大气改正方法。该方法被应用于链子崖危岩体监测实验中,结果表明,其能有效削弱大气效应的影响。
由于现场气象测量条件有限,本次气象数据仅测量了雷达站和变形监测目标处的气象数据。根据链子崖特殊的地形条件(视线经过U 字型山谷地段,测区紧邻长江),两端气象测量并不能准确反映该区域实际的气象情况。采用更精确的气象仪,同时测量雷达和变形目标之间的气象值,以获取更准确的局部气象模型,从而解除对地面稳定点的依赖,实现SAR 数据准确的大气改正,需要继续研究。