更新日志
- ver.0613 前四部分更新完成,InSAR部分 参考面模型结束
- ver.0614 高程模型结束
- ver.0617 差分干涉以及干涉误差分析完成
- ver.0619 配准,基线估计,滤波解缠,时序InSAR完成
- ver.0620 地基雷达基本完成
第一部分 InSAR测绘技术(微波、雷达基础了解)
本章要求:1.熟悉典型SAR系统(星载系统)2.掌握SAR的发展趋势
微波遥感
1.微波成像:通过微波传感器获取目标地物发射或反射的微波辐射成像
2.1952年,Carl Wiley提出了多普勒波束锐化可以提高成像分辨率,由此诞生合成孔径雷达即SAR。优点是技术上使小天线实现更大的天线的有效分辨率。缺点是必须处理高计算负载的复杂算法。
雷达波段
波段 | 波长 |
---|---|
Ka | <1.5cm |
K | 1.5cm |
Ku | >1.5cm |
X | 3cm |
C | 5cm |
S | 10cm |
L | 23cm |
P | <1m |
SAR发展历程
Satellite | Launch Time | Resolution | Wave Length | Incident Angle |
---|---|---|---|---|
Seasat | June 26,1978 | 25m×25m | 0.235m | ~23° |
RadarSat2 | Dec 24,2007 | 18m×24.7m | 0.05m | 20°~49° |
Gaofen-3satellite | Aug 10,2016 | 1m~500m | 0.05m | 20°~50° |
特点:分辨率逐渐提高;成像更灵活;频段设置更合理
SAR发展趋势
单基到多基
优势
1.收发分置,视角可变。可以获取目标多角度的散射信息,实现多角度融合。
2.基线配置机动灵活。根据不同观测区域的测量精度要求,灵活调整基线长度,获取地面高程和地面运动目标信息。
3.静默接受,隐蔽性强。接收机无需发射电磁波,在现代战争中不易被对方的侦查装置侦查到。
4.系统构型多样,收发系统可搭载于卫星,飞机,地面装置等。未来可实现“一星发射多星接收”的分布式多基SAR系统,可同时实现高分辨率和超宽幅成像。
高分辨率宽幅
多维度
SAR成像利用了电磁波信息具有的频率、幅度、极化、相位等特性。可以从SAR图像中获取观测目标的多维度信息。极化SAR技术可以获取地物取向,形状,粗糙度,介电常数等物理特性;干涉SAR技术可以获取场景的高精度数字高程模型、洋流测速、冰川位移、地表形变等。另外还有极化干涉SAR技术和SAR层析技术等。
小型化
略
第二部分 雷达概述
雷达定义
雷达是Radar的音译,Radar是Radio detection and ranging的简称,即无线电探测和测距。
由IEEE(电气和电子工程师协会)对雷达的标准定义为
雷达是通过发射电磁波信号,接收来自其威力覆盖范围内目标的回波,并从回波信号中提取位置和其他信息,以用于探测、定位,以及有时进行目标识别的电磁系统
雷达原理
基本工作原理
由雷达发射机产生的电磁能,经收发开关后传输给天线,再定向辐射于大气中,如果目标位于定向天线波束内,截取一部分电磁能,再将这些截取能量向各方向散射,部分能量进入到雷达接收机。接收机将散射回波信号经信号处理送终端显示。
测速、测方位
略
测距
1.脉冲测距法:测时间
R
=
c
∗
Δ
t
2
R=c*\dfrac{\Delta t}{2}
R=c∗2Δt
2.相位测距法:测相位
略
雷达方程
P
r
=
P
t
G
t
2
λ
2
σ
(
4
π
)
3
R
4
P_r=\dfrac{P_tG_t^2\lambda ^2\sigma}{(4\pi)^3R^4}
Pr=(4π)3R4PtGt2λ2σ
P
r
P_r
Pr:反射功率;
P
t
P_t
Pt:发射功率;
G
t
G_t
Gt:天线增益;
λ
\lambda
λ:反射功率;
σ
\sigma
σ:目标有效面积;
R
R
R:距天线距离;
由公式可求雷达最大作用距离。
第三部分 微波散射特征
散射的定义和理论
定义:电磁波辐射在非均匀介质或各向异性介质中传播时,多方位、多角度地改变原来的传播方向的现象。也即目标对入射电磁波能量的重定向。反射、折射和衍射都是散射的基本形式。
理论:
1.瑞丽散射( a < 0.1 λ a<0.1\lambda a<0.1λ):散射粒子直径远远小于入射光波长,产生的一种没有频率位移(即无能量变化,波长相同)的弹性光散射,其散射强度与方向有关。
2.米氏散射( 0.1 λ < a < 10 λ 0.1\lambda<a<10\lambda 0.1λ<a<10λ):散射粒子直径与入射光波长相当时发生的一种有更强方向性的散射。表现为前向散射强度远高于后向散射。也称为谐振散射
地物的散射特性
地物的散射特性首先取决于地物尺寸(粗糙度)和雷达波长间的关系,此外还有,入射角,介电特性和极化特性。
- 粗糙度:
- 瑞利区:目标尺寸远小于波长
- 米氏区:目标尺寸与波长相当
- 光学区:目标尺寸远大于波长
因此同样的粗糙程度的地表,长波相比短波,短波会认为地表更粗糙,从而背向散射更强,图像上也显示得更亮。
粗糙
光滑
- 入射角:雷达入射波束与地表发现的夹角。大入射角会使得回波能量低,图像显示暗,小角度反之。
- 介电常数:描述材料的电性质(电容,传导率,反射率)。通常定义为物体电容与真空电容之比。介电常数越大反射越强。
- 干燥物体介电常数在3~5之间
- 水体接近80
- 岩石之间介电常数很小
- 波长:不同工作频段会产生不同的透射效应。频段越低,投射效应越强(可以类比衍射来理解)
第四部分 侧视雷达(SLAR/SAR了解)
侧视雷达成像
侧视雷达SLR,side looking radar。成像中心与竖直方向呈一定角度。
其成像目标在像平面的位置,与目标距天线的距离相关
侧视雷达分类
真实孔径雷达
- 两个重要概念,斜距(Slant Range),地距(Ground Range)
- 斜距是雷达观测的自然观测系下,雷达视线方向到目标的距离
- 地距是斜距在地球大地水准面上的投影
- 经过重采样等处理,可以将斜距坐标改变成地距坐标
- 地距影像与斜距影像差异。可以发现相等的地面距离,在入射角小的情况下,斜距差距小,在图像反映为压缩。(可以用余弦定理推导)
-
两个重要参数,方位向分辨率,距离向分辨率
雷达侧视扫描方向称为距离向,行进方向称为方位向。
- 距离向分辨率:即在垂直于飞行方向上,能够分辨目标产生回波的最小距离。其物理意义即要使两目标回波不完全重叠的的脉冲宽度
τ
\tau
τ,要满足
τ
<
2
Δ
R
c
\tau<\dfrac{2\Delta R}{c}
τ<c2ΔR,因此在斜距上就应当满足:
R
d
=
c
τ
2
R_d = \dfrac {c\tau}{2}
Rd=2cτ距离向分辨率考虑的是地距,因此将斜距分辨率转为地距分辨率即为:
R
g
=
c
τ
2
cos
θ
=
c
τ
2
sec
θ
R_g = \dfrac {c\tau}{2\cos\theta}=\dfrac{c\tau}{2}\sec\theta
Rg=2cosθcτ=2cτsecθ式中
θ
\theta
θ为传感器俯角(或是雷达入射角的余角)
+ 即距离向分辨率与距离无关,与脉冲宽度 τ \tau τ有关,要提高分辨率就要减小脉冲宽度。
+ 俯角越大的地方,分辨率越低,俯角为90°,也即是入射角为0°时,雷达中心投影时,无法分辨相邻地面点,这可以解释成像雷达为何采用侧视成像。
- 方位向分辨率:即雷达在行进方向上的最小距离,即
Δ
L
=
λ
R
D
\Delta L = \frac{\lambda R}{D}
ΔL=DλR
式中, D D D为真实孔径大小, R 为斜距 R为斜距 R为斜距, λ \lambda λ为波长
- 距离向分辨率:即在垂直于飞行方向上,能够分辨目标产生回波的最小距离。其物理意义即要使两目标回波不完全重叠的的脉冲宽度
τ
\tau
τ,要满足
τ
<
2
Δ
R
c
\tau<\dfrac{2\Delta R}{c}
τ<c2ΔR,因此在斜距上就应当满足:
R
d
=
c
τ
2
R_d = \dfrac {c\tau}{2}
Rd=2cτ距离向分辨率考虑的是地距,因此将斜距分辨率转为地距分辨率即为:
R
g
=
c
τ
2
cos
θ
=
c
τ
2
sec
θ
R_g = \dfrac {c\tau}{2\cos\theta}=\dfrac{c\tau}{2}\sec\theta
Rg=2cosθcτ=2cτsecθ式中
θ
\theta
θ为传感器俯角(或是雷达入射角的余角)
可以发现,方位向分辨率的提升需要选择更短波长、更大孔径、更小观测距离。无论是在卫星还是飞机上都难以实现。因此合成孔径雷达技术闪亮登场
合成孔径雷达
定义:雷达安装在运动平台(卫星、飞机和滑轨等)上,以一定的频率不断地发射和接受电磁波脉冲,每发射一次脉冲时的天线位置视为阵列天线的单元阵子位置,最后将这些位置上不同时刻存在的单元阵子组合起来,形成一个等效的大孔径天线,从而得到提高的方位向分辨率。
用人话说就是一边移动对一个目标多次观测,信息相干相加,即模拟出了一个超长天线,以此来达到提高方位向分辨率的目的,举例说明比如,昆虫复眼、夜晚多张照片合成一张白天的照片、多个小雷达合成大雷达
那么合成孔径技术模拟出来的超长天线长度为
L
s
=
β
R
=
λ
D
R
L_s=\beta R = \frac{\lambda}{D} R
Ls=βR=DλR新的波束宽度为
β
s
=
λ
2
⋅
L
s
\beta_s=\frac{\lambda}{2\cdot L_s}
βs=2⋅Lsλ则新的**方位向分辨率**为
Δ
L
s
=
β
s
⋅
R
=
λ
⋅
R
2
⋅
L
s
=
λ
2
⋅
β
=
D
2
\Delta L_s=\beta_s\cdot R=\frac{\lambda\cdot R}{2\cdot L_s}=\frac{\lambda}{2\cdot\beta}=\frac{D}{2}
ΔLs=βs⋅R=2⋅Lsλ⋅R=2⋅βλ=2D即方位向分辨率仅与孔径大小有关,解决了之前需要更短波长、更大孔径、更小观测距离 的问题。
- 星载SAR的成像模式
- 条带成像
距离向上通过线性调频信号提高分辨率,方位向上用合成孔径原理提高分辨率。 - ScanSAR
有多个子测绘带 - 聚焦模式(SPOT SAR)
瞄准观测,以达到高分辨率,比条带高 - TOPS SAR
照射方向从后向前,与聚焦模式相反,故得名 - 滑动聚焦
解决了传统聚焦模式的照射范围有限以及方位向照射不均匀,中心与边缘目标能量差异大的缺点
- 条带成像
SAR影像的几何特征
结合侧视雷达成像原理学习
-
距离向压缩
相等的地距,大于斜距,且随着入射角增大而增大
斜距与地距有如下关系: R s = R g sin θ , θ R_s=R_g\sin\theta,\theta Rs=Rgsinθ,θ 为雷达入射角(incident angle) -
透视收缩(foreshortening)
针对斜坡的几何特征
可以从图中直观的看出,面向雷达的斜坡,其底部 A A A点和顶部 B B B的斜距差 Δ R \Delta R ΔR远小于地距差 Δ X \Delta X ΔX
其比值称为压缩比 l = Δ R Δ X = sin ϕ cos α l=\frac{\Delta R}{\Delta X}=\frac {\sin\phi}{\cos\alpha} l=ΔXΔR=cosαsinϕ。- 透视压缩的形成条件:1.侧视角大于坡度角,二者相等时,压缩比达到极小值。2.雷达俯角加坡度角小于等于90°
eg:
- 透视压缩的形成条件:1.侧视角大于坡度角,二者相等时,压缩比达到极小值。2.雷达俯角加坡度角小于等于90°
-
顶底位移(layover)
当透视压缩条件2不满足,即雷达俯角与地面坡度角之和大于90°时,变回出现,坡度顶部斜距小于底部斜距,从而在图像呈现斜坡顶部与底部颠倒的情况,这种情况称为顶底位移,也称为雷达叠掩
易知,此现象与雷达俯角(或入射角)和 局部入射角(雷达视线与实际地面的夹角)密切相关。入射角越小,越有可能发生顶底位移,因此顶底位移多是近距离现象。而与 局部入射角 的关系如下表:
局部入射角 θ \theta θ | 现象 |
---|---|
θ < 0 \theta<0 θ<0 | 顶底位移 |
0 < θ < 90 ° 0<\theta<90° 0<θ<90° | 透视压缩 |
θ > 90 ° \theta>90° θ>90° | 阴影 |
- 雷达阴影(shadow)
如上述,当局部入射角 θ > 90 ° \theta>90° θ>90°时,就会出现雷达阴影情况,意味着传感器没有接收到回波,图像显示黑色区域。
注:因为,雷达阴影的产生除了坡度角问题,所以判断阴影产生还要考虑山脊走向或是卫星航向的关系
第五部分 雷达干涉测量技术(InSAR!)
物理背景
干涉理论基础是杨氏双缝干涉实验,即满足一定条件的两列相干光相遇叠加,在叠加区域出现某些点振动始终加强和始终减弱的情况,即在相干区域内振动强度有稳定的空间分布 振动加强和减弱的结果本质在于波的相位叠加结果
其中公式先要牢记,后面会用的很频繁
δ
=
2
π
λ
(
r
2
−
r
1
)
=
2
π
λ
d
x
D
\delta=\frac{2\pi}{\lambda}(r_2-r_1)=\frac{2\pi}{\lambda}d\frac{x}{D}
δ=λ2π(r2−r1)=λ2πdDx式中
δ
\delta
δ是相位差,相当重要
测量概述
SAR成像过程中,通过相位式测距法,精度在cm-mm级,并记录回波中的强度信号,将之记录在像平面上,在一个平面上记录一个二维信息,采用复数形式记录,即SLC(single looking complex)单数复数数据。
同一区域如果拥有两张SLC,就可以做干涉,从而获取地面高程、形变等信息。
两幅图像上同一点的数据可以表示为
u
1
=
∣
A
1
∣
e
j
ϕ
1
u_1 = \vert A_1 \vert e^{j\phi_1}
u1=∣A1∣ejϕ1
u
2
=
∣
A
2
∣
e
j
ϕ
2
u_2=\vert A_2 \vert e^{j\phi_2}
u2=∣A2∣ejϕ2
要求复数间的幅角差可以采用共轭相乘的办法,即:
u
1
⋅
u
2
‾
=
∣
A
1
∣
∣
A
2
∣
e
j
(
ϕ
1
−
ϕ
2
)
u_1\cdot \overline{u_2}=\vert A_1 \vert\vert A_2 \vert e^{j(\phi_1-\phi_2)}
u1⋅u2=∣A1∣∣A2∣ej(ϕ1−ϕ2)进而可以求得
u
1
⋅
u
2
‾
u_1\cdot \overline{u_2}
u1⋅u2的主幅角
Φ
\Phi
Φ。注:
ϕ
1
、
ϕ
2
、
Φ
\phi_1、\phi_2、\Phi
ϕ1、ϕ2、Φ作为主幅角,取值范围为
(
−
π
,
π
]
(-\pi,\pi]
(−π,π]
解决了数学问题,现在应当解决如何获得两张SLC的问题了
- 可以对着同一片区域拍两次。有点类似SAR成像模式中讲的聚焦模式(SPOT SAR)。这种称为重复轨道单天线模型
- 基本假设:两次拍摄间距,即空间基线 B B B应远远小于斜距 R R R这样可以保证两次成像时,目标视线几乎平行
- 不足在于:在星载SAR中,两次拍摄的时间基线长,容易影像干涉质量,即失相干。
- 还可以有两个天线,同一时间对着同一目标拍两次,这样时间基线为零,这种模式称为单轨双天线横向模型
- 优势在于影像配准比较容易,缺点是空间基线选择余地有限,受飞行平台几何条件限制。
函数模型
要求掌握:1.理解并推导四种函数模型。2.重要名词:干涉基线、平地效应、高程模糊度。
- 参考面模型:假设目标点
P
P
P位于大地水准面,忽略形变,大气和噪声
如图, S 1 S_1 S1, S 2 S_2 S2,分别为两天线相位中心, R 1 R_1 R1, R 2 R_2 R2为对应斜距, B B B为基线 , B ∥ B_\parallel B∥, B ⊥ B_\bot B⊥分别是基线的延视线方向的水平分量和垂直于视线方向的垂直分量, θ \theta θ是雷达入射角, α \alpha α是基线水平角。
在 △ P S 1 S 2 \triangle PS_1S_2 △PS1S2中,基于余弦定理有如下等式: cos ∠ P S 1 S 2 = R 1 2 + B 2 − R 2 2 2 B R 1 = ( R 1 + R 2 ) ( R 1 − R 2 ) + B 2 2 B R 1 = δ R ⋅ ( R 1 + R 2 ) B ⋅ 2 R 1 + B 2 R 1 \cos\angle PS_1S_2 =\frac{R_1^2+B^2-R_2^2}{2BR_1}=\frac{(R_1+R_2)(R_1-R_2)+B^2}{2BR_1}=\frac{\delta R\cdot(R_1+R_2)}{B\cdot2R_1}+\frac{B}{2R_1} cos∠PS1S2=2BR1R12+B2−R22=2BR1(R1+R2)(R1−R2)+B2=B⋅2R1δR⋅(R1+R2)+2R1B又因为 B ≪ R B\ll R B≪R,即 R 1 + R 2 ≈ 2 R 1 R_1+R_2\approx2R_1 R1+R2≈2R1, B R 1 ≈ 0 , \frac{B}{R_1}\approx0, R1B≈0,进一步可得 cos ∠ P S 1 S 2 ≈ δ R B \cos\angle PS_1S_2 \approx\frac{\delta R}{B} cos∠PS1S2≈BδR而 ∠ P S 1 S 2 = π 2 + α − θ \angle PS_1S_2=\frac{\pi}{2}+\alpha-\theta ∠PS1S2=2π+α−θ,有: cos ∠ P S 1 S 2 = cos ( π 2 + α − θ ) = sin ( θ − α ) \cos\angle PS_1S_2=\cos(\frac{\pi}{2}+\alpha-\theta)=\sin(\theta-\alpha) cos∠PS1S2=cos(2π+α−θ)=sin(θ−α),联立上式: R 1 − R 2 = δ R ≈ B ⋅ sin ( θ − α ) = B ∥ R_1-R_2 = \delta R\approx B\cdot\sin(\theta-\alpha)=B_\parallel R1−R2=δR≈B⋅sin(θ−α)=B∥
通过上述推导,我们知道了基线为
B
B
B的InSAR对于目标点
P
P
P的斜距差
δ
R
\delta R
δR,再通过前文InSAR物理背景中的相位差公式(考虑回波,相位差所对应的距离,是斜距差的两倍),我们可以得出
P
P
P点相位差
ϕ
p
\phi_p
ϕp,为:
ϕ
p
=
2
⋅
2
π
λ
(
R
2
−
R
1
)
=
−
4
π
λ
B
sin
(
θ
−
α
)
\phi_p=2\cdot\frac{2\pi}{\lambda}(R_2-R_1)=-\frac{4\pi}{\lambda}B\sin(\theta-\alpha)
ϕp=2⋅λ2π(R2−R1)=−λ4πBsin(θ−α)
现在,解决了大地水准面上的单一点相位差,那水准面(同一高度)上的两点间对应的干涉相位是
ϕ 1 = − 4 π λ B sin ( θ − α ) \phi_1=-\frac{4\pi}{\lambda}B\sin(\theta-\alpha) ϕ1=−λ4πBsin(θ−α) ϕ 2 = − 4 π λ B sin ( θ + Δ θ − α ) \phi_2=-\frac{4\pi}{\lambda}B\sin(\theta+\Delta\theta-\alpha) ϕ2=−λ4πBsin(θ+Δθ−α)
如图
二者差分可得:
Δ
ϕ
=
−
4
π
λ
[
B
sin
(
θ
0
+
Δ
θ
−
α
)
−
B
sin
(
θ
0
−
α
)
]
≈
−
4
π
λ
Δ
θ
B
cos
(
θ
0
−
α
)
\Delta\phi=-\frac{4\pi}{\lambda}[B\sin(\theta_0+\Delta\theta-\alpha)-B\sin(\theta_0-\alpha)]\approx-\frac{4\pi}{\lambda}\Delta\theta B\cos(\theta_0-\alpha)
Δϕ=−λ4π[Bsin(θ0+Δθ−α)−Bsin(θ0−α)]≈−λ4πΔθBcos(θ0−α)此处,因
Δ
θ
\Delta\theta
Δθ极小,故而用到了
Δ
θ
≈
sin
Δ
θ
\Delta\theta\approx\sin\Delta\theta
Δθ≈sinΔθ,
cos
Δ
θ
≈
1
\cos\Delta\theta\approx1
cosΔθ≈1。且有上图几何关系可知,
R
Δ
θ
≈
R
sin
Δ
θ
≈
Δ
R
tan
θ
0
R\Delta\theta\approx R\sin\Delta\theta\approx\frac{\Delta R}{\tan\theta_0}
RΔθ≈RsinΔθ≈tanθ0ΔR
Δ
R
\Delta R
ΔR为斜距差。进而将上式进一步推导得:
Δ
ϕ
=
−
4
π
λ
B
cos
(
θ
0
−
α
)
Δ
R
R
tan
θ
0
≈
4
π
B
⊥
Δ
R
λ
R
tan
θ
0
\Delta\phi=-\frac{4\pi}{\lambda}\frac{B\cos(\theta_0-\alpha)\Delta R}{R\tan\theta_0}\approx\frac{4\pi B_\bot\Delta R}{\lambda R\tan\theta_0}
Δϕ=−λ4πRtanθ0Bcos(θ0−α)ΔR≈λRtanθ04πB⊥ΔR
可以知道,当地面两个点高度相同时,对应的相位差仍然存在,并与斜距差
Δ
R
\Delta R
ΔR成正比,因此即使是水平地表,干涉相位图也表现为与方位向平行的周期性变化的条纹,这部分条纹称之为平地效应
-
高程函数模型
实际情况,目标点 P P P不会是大地水准面上点,而是有高程 h h h,其对应几何关系如图:
依据相位差公式,可知目标点 P P P的相位差为: ϕ p = − 4 π λ B sin ( θ 1 − α ) \phi_p=-\frac{4\pi}{\lambda}B\sin(\theta_1-\alpha) ϕp=−λ4πBsin(θ1−α)
做一条斜距与 R 1 R_1 R1相等的直线与地面相交于点 C C C。由于 B ≪ R B\ll R B≪R(或是 h ≪ H h\ll H h≪H),所以 ∠ C S 1 P = Δ θ \angle CS_1P=\Delta\theta ∠CS1P=Δθ为极小值,所以 ∠ S 1 C P = ∠ S 1 P C ≈ 90 ° \angle S_1CP=\angle S_1PC\approx90° ∠S1CP=∠S1PC≈90°进而得, R sin Δ θ ≈ h sin θ R\sin\Delta\theta\approx\frac{h}{\sin\theta} RsinΔθ≈sinθh。那么P点的相位差就可以进一步写为: ϕ p = − 4 π λ B sin ( θ + Δ θ − α ) = − 4 π λ B sin ( θ − α ) − 4 π λ B cos ( θ − α ) sin Δ θ \phi_p=-\frac{4\pi}{\lambda}B\sin(\theta+\Delta\theta-\alpha)=-\frac{4\pi}{\lambda}B\sin(\theta-\alpha)-\frac{4\pi}{\lambda}B\cos(\theta-\alpha)\sin\Delta\theta ϕp=−λ4πBsin(θ+Δθ−α)=−λ4πBsin(θ−α)−λ4πBcos(θ−α)sinΔθ然后: ϕ p = − 4 π λ B ∥ − 4 π λ B ⊥ h R sin θ \phi_p=-\frac{4\pi}{\lambda}B_\parallel-\frac{4\pi}{\lambda}\frac{B_\bot h}{R\sin\theta} ϕp=−λ4πB∥−λ4πRsinθB⊥h
最终,我们得到了,不考虑大气误差,噪声,地面地物形变引起的干涉相位的**高程模型,在大地水准面上高程为 h h h的以目标点 P的干涉相位差等于入射角为 θ \theta θ( C C C点)对应的参考面相位差加上高程 h h h带来高程相位差**
推导过程中用了很多涉及到
Δ
θ
\Delta\theta
Δθ的近似关系,使用matlab进行了简易计算发现,
Δ
θ
\Delta\theta
Δθ对相位差的影响在10-8数量级,足够小,可以接受。
*注:多说几句唠叨话
1.需要明确的是
C
C
C与
P
P
P点在SAR图像上因其斜距相同,应当是同一个点,很容易寻找。
2.因为1中所言,所以做
S
1
C
S_1C
S1C这条线是为了将模型推导为两个部分,一是参考面相位,另一是高程引起的相位。实际上不存在如此一条线。
*
高程模型:目标点
P
P
P相位是有参考面相位和高程相位组成:
ϕ
p
=
ϕ
R
e
f
e
r
e
n
c
e
+
ϕ
e
l
e
v
a
t
i
o
n
\phi_p=\phi_{Reference}+\phi_{elevation}
ϕp=ϕReference+ϕelevation可以据此消除SAR影响的平地效应,从而放大研究区域高程带来的相位信息,如图:
消除平地效应后,高程与相位的关系:
ϕ
e
l
e
v
a
t
i
o
n
=
−
4
π
λ
B
⊥
R
sin
θ
H
p
\phi_{elevation}=-\frac{4\pi}{\lambda}\frac{B_\bot }{R\sin\theta}H_p
ϕelevation=−λ4πRsinθB⊥Hp当
ϕ
e
l
e
v
a
t
i
o
n
=
2
π
\phi_{elevation}=2\pi
ϕelevation=2π,即一个条纹间距所对应的高程
H
2
π
=
−
λ
R
sin
θ
2
B
⊥
H_{2\pi}=-\frac{\lambda R\sin\theta}{2B_\bot}
H2π=−2B⊥λRsinθ。我们将
H
2
π
H_{2\pi}
H2π称为高程模糊度,从式中可以看出,高程模糊度与垂直基线
B
⊥
B_\bot
B⊥成反比,因此,要想获得一个高精度的DEM,即高程模糊度小的DEM,就需要长基线。注:但基线不能无限长,后续再提
- 形变模型
假设一目标点 P P P,两次SAR数据收集间隔,无高程变化下,产生了微小形变 r r r,如图:
第一次采集的 P P P点斜距为 R 1 ′ R_1^{'} R1′,第二次为 R 2 ′ ′ R_2^{''} R2′′,两次 P P P点的相位差为: ϕ p = − 4 π λ ( R 1 ′ − R 2 ′ ′ ) \phi_p=-\frac{4\pi}{\lambda}(R_1^{'}-R_2^{''}) ϕp=−λ4π(R1′−R2′′)因为InSAR的基本假设(见测量概述部分),以及图中圆弧线,可以认为: R 2 ′ ′ = R 1 ′ + Δ r = R 2 ′ + Δ r R_2^{''}=R_1^{'}+\Delta r=R_2^{'}+\Delta r R2′′=R1′+Δr=R2′+Δr因此 ϕ p = − 4 π λ ( R 1 ′ − R 2 ′ ) + 4 π λ Δ r = ϕ p = − 4 π λ B ∥ − 4 π λ B ⊥ h R sin θ + 4 π λ Δ r \phi_p=-\frac{4\pi}{\lambda}(R_1^{'}-R_2^{'})+\frac{4\pi}{\lambda}\Delta r=\phi_p=-\frac{4\pi}{\lambda}B_\parallel-\frac{4\pi}{\lambda}\frac{B_\bot h}{R\sin\theta}+\frac{4\pi}{\lambda}\Delta r ϕp=−λ4π(R1′−R2′)+λ4πΔr=ϕp=−λ4πB∥−λ4πRsinθB⊥h+λ4πΔr即 ϕ p = ϕ R e f e r e n c e + ϕ e l e v a t i o n + ϕ d e f o r m a t i o n \phi_p=\phi_{Reference}+\phi_{elevation}+\phi_{deformation} ϕp=ϕReference+ϕelevation+ϕdeformation这就是InSAR的形变函数模型
注:可以比较 ϕ e l e v a t i o n \phi_{elevation} ϕelevation与 ϕ d e f o r m a t i o n \phi_{deformation} ϕdeformation的系数,因为基线远远小于斜距的原因,前者系数远远小于后者,这也是为什么InSAR对形变的敏感度远高于高程敏感度。因此,InSAR技术在形变监测上应用广泛。
-
一般模型
除了参考面,高程,形变,往往还要考虑大气延迟误差和噪声对相位数据采集的影响-
大气延迟误差
- 电离层
内容略,给个大概公式: ϕ ( f ) = − 2 π f 2 1 0 6 ∫ 0 H N i o n o ( f , h ) c d h ≈ 4 π λ K ⋅ T E C \phi(f)=-2\pi f\frac{2}{10^6}\int_0^H\frac{N_{iono}(f,h)}{c}dh\approx\frac{4\pi}{\lambda}K\cdot TEC ϕ(f)=−2πf1062∫0HcNiono(f,h)dh≈λ4πK⋅TEC- 对流层
可以基于InSAR数据进行改正,或是地面模型,或是大气数据
- 对流层
- 电离层
-
噪声误差
略 -
一般模型
综上,InSAR的完整函数模型包含参考面相位,高程相位,形变相位,大气误差以及噪声影响 ϕ p = ϕ R e f e r e n c e + ϕ e l e v a t i o n + ϕ d e f o r m a t i o n + ϕ a t m + ϕ n o i s e \phi_p=\phi_{Reference}+\phi_{elevation}+\phi_{deformation}+\phi_{atm}+\phi_{noise} ϕp=ϕReference+ϕelevation+ϕdeformation+ϕatm+ϕnoise或 ϕ = − 4 π λ B ∥ − 4 π λ B ⊥ h R sin θ + 4 π λ Δ r + 4 π λ S t 1 − 4 π λ S t 2 + 2 π w k \phi=-\frac{4\pi}{\lambda}B_\parallel-\frac{4\pi}{\lambda}\frac{B_\bot h}{R\sin\theta}+\frac{4\pi}{\lambda}\Delta r+\frac{4\pi}{\lambda}S^ {t_1}-\frac{4\pi}{\lambda}S^ {t_2}+2\pi w_k ϕ=−λ4πB∥−λ4πRsinθB⊥h+λ4πΔr+λ4πSt1−λ4πSt2+2πwk
若是将各点收集到的SAR数据写成矩阵形式,并求解个未知数 [ ϕ 1 w ϕ 2 w ⋮ ϕ n w ] = [ A 1 A 2 ⋱ A n ] [ X 1 X 2 ⋮ X n ] \begin{bmatrix} \phi_1^w\\ \phi_2^w\\ \vdots \\ \phi_n^w \end{bmatrix}= \begin{bmatrix} A_1&&& \\ &A_2&&\\ &&\ddots&\\ &&&A_n \end{bmatrix} \begin{bmatrix} X_1\\X_2\\\vdots\\X_n \end{bmatrix} ϕ1wϕ2w⋮ϕnw = A1A2⋱An X1X2⋮Xn
系数矩阵中的 A k = [ − 4 π λ B ⊥ h R sin θ 4 π λ 4 π λ − 4 π λ 2 π ] A_k=\begin{bmatrix}-\frac{4\pi}{\lambda}\frac{B_\bot h}{R\sin\theta}&\frac{4\pi}{\lambda}&\frac{4\pi}{\lambda}&-\frac{4\pi}{\lambda}&2\pi\end{bmatrix} Ak=[−λ4πRsinθB⊥hλ4πλ4π−λ4π2π], X k = [ H Δ r S 1 k S 2 k 2 π ] X_k=\begin{bmatrix}H&\Delta r&S_1^k&S_2^k&2\pi\end{bmatrix} Xk=[HΔrS1kS2k2π](参考面相位为常数)
可以发现, n n n个观测值需要求解 5 n 5n 5n个未知数,没有唯一解。因此需要采取其他办法,消除未知数,这个方法就差 分 干 涉
-
差分干涉(DInSAR)
根据InSAR的函数模型,一张干涉影像上,没有办法对所有未知参数做到求解。解决这一问题的思路有两个:
1.加入更多观测值。
2.借助于先验信息。
二轨差分干涉测量技术(借助先验知识)
- 二轨可以理解为卫星两次经过同一区域所采集的SAR影像做差分干涉,即两幅不同一时间的SAR影像进行地形形变求解。
x k = [ H k D S k 1 S k 2 w k ] x_k=\begin{bmatrix}\boxed {H_k}&D &\boxed {S^1_k} &\boxed {S^2_k} &\boxed {w_k} \end{bmatrix} xk=[HkDSk1Sk2wk]借助先验知识,未知量中 H k H_k Hk用已知DEM计算,对大气误差采取忽略和减弱的策略,对噪声影响采用建模消除的策略。从而实现未知量从五个变为一个的转换,最终实现地形形变的测量计算。大致流程图如下:
三轨差分干涉测量技术(加入更多观测值)
- 在二轨差分基础上,增加一幅新的SAR影像,新影像要求与另外两张中的一张时间基线为零,即同时拍摄。
人话就是:三张SAR影像,其中两张同时拍摄,另一张不同时。
两张时间基线为0的SAR影像可以组成干涉对,其相位表示为(
m
=
4
π
λ
m=\frac{4\pi}{\lambda}
m=λ4π):
Φ
t
o
p
o
=
−
m
B
⊥
,
k
t
o
p
o
R
k
1
sin
θ
k
1
H
k
+
m
S
k
1
−
m
S
k
2
\Phi_{topo}=-m\frac{B^{topo}_{\bot,k}}{R^1_k\sin\theta^1_k}H_k+mS^1_k-mS^2_k
Φtopo=−mRk1sinθk1B⊥,ktopoHk+mSk1−mSk2
此干涉对称为地形干涉对,特点在于时间基线为零,其相位差不包含形变相位,噪声影响也通过差分消除或减弱了。
再加入剩下的那张SAR影像组成干涉对,其相位表示为:
Φ
d
e
f
o
=
−
m
B
⊥
,
k
d
e
f
o
R
k
2
sin
θ
k
2
H
k
+
D
k
+
m
S
k
1
−
m
S
k
3
−
2
π
w
k
\Phi_{defo}=-m\frac{B^{defo}_{\bot,k}}{R^2_k\sin\theta^2_k}H_k+D_k+mS^1_k-mS^3_k-2\pi w_k
Φdefo=−mRk2sinθk2B⊥,kdefoHk+Dk+mSk1−mSk3−2πwk
此干涉对称为形变干涉对,其相位组成包含了形变相位。将上述两式简要结合,可得
Φ
d
e
f
o
=
B
⊥
,
k
d
e
f
o
B
⊥
,
k
t
o
p
o
Φ
t
o
p
o
+
(
1
−
B
⊥
,
k
d
e
f
o
B
⊥
,
k
t
o
p
o
)
m
S
k
1
+
B
⊥
,
k
d
e
f
o
B
⊥
,
k
t
o
p
o
m
S
k
2
−
m
S
k
3
+
m
D
k
−
2
π
w
k
\Phi_{defo}=\frac{B^{defo}_{\bot,k}}{B^{topo}_{\bot,k}}\Phi_{topo}+(1-\frac{B^{defo}_{\bot,k}}{B^{topo}_{\bot,k}})mS^1_k+\frac{B^{defo}_{\bot,k}}{B^{topo}_{\bot,k}}mS^2_k-mS^3_k+mD_k-2\pi w_k
Φdefo=B⊥,ktopoB⊥,kdefoΦtopo+(1−B⊥,ktopoB⊥,kdefo)mSk1+B⊥,ktopoB⊥,kdefomSk2−mSk3+mDk−2πwk
地形干涉对选择要求:1. B t o p o ≥ B d e f o B^{topo}\ge B^{defo} Btopo≥Bdefo,长基线可以保证较高的高程精度。2. B t o p o ≤ 0.7 B c r i t B^{topo}\le 0.7B^{crit} Btopo≤0.7Bcrit,可以保证较高的相干值数据。
干涉测量的误差源及误差分析
1.相干系数和相干值的区别
2.影响误差的因素
SAR采集的某点复数数据统计量分别用 u 1 u_1 u1和 u 2 u_2 u2表示,对应的干涉相位 v = u 1 u 2 ∗ v=u_1u_2^* v=u1u2∗,其相关系数为 γ = E [ u 1 u 2 ∗ ] E [ ∣ u 1 ∣ 2 ] E [ ∣ u 2 ∣ 2 ] \gamma=\frac{E[u_1u_2^*]}{\sqrt{E[|u_1|^2]E[|u_2|^2]}} γ=E[∣u1∣2]E[∣u2∣2]E[u1u2∗]此系数,称为相干系数,相干系数越低说明干涉相位中噪声越多,也说明干涉影像的质量越低,反之同理。
相干值
相干系数求解结果是一个复数,因此提出了一种估计相干系数的办法,即相干值的概念,定义:
γ
^
[
i
,
j
]
=
∣
∑
W
u
1
[
i
,
j
]
u
2
∗
[
i
,
j
]
∣
∑
W
∣
u
1
[
i
,
j
]
∣
2
∑
W
∣
u
2
[
i
,
j
]
∣
2
\hat \gamma[i,j]=\frac {{\bigg|\sum\limits_{W}u_1[i,j]u_2^*[i,j]\bigg|}}{\sqrt{\sum\limits_W|u_1[i,j]|^2\sum\limits_W|u_2[i,j]|^2}}
γ^[i,j]=W∑∣u1[i,j]∣2W∑∣u2[i,j]∣2
W∑u1[i,j]u2∗[i,j]
W
W
W是估计窗口,一般包括5~100像素。且要求像素内平稳。需要注意的是,Touzi(1999)表明相干值对相干系数的估计是有偏的,主要表现为估计窗口越大,偏差越低;相干系数幅度越大,偏差越低;最大似然估计可以做到无偏。
新概念——多视视数:SAR的单视复数数据(SLC)是原始的最高分辨率数据,但是从单个像元散射的雷达回波信号的相干叠加,导致强度信息有很多噪声。多视处理是对SLC数据方位向和/或距离向做平均,得到的结果是多视后的强度数据。通过多视处理的SLC数据,空间分辨率降低,提升了数据的辐射分辨率,也就是强度信息。(可以简单的人话概括,原本分辨率距离向和方位向不等,类似矩形,现在多来几个矩形拼在一起拼成一个“正方形”,比如1.5 * 5的分辨率,来凑3个,凑成4.5 * 5的分辨率,而3就是多视视数。这样处理可以消除原本单个像元里的强度噪声,使得图像的辐射分辨率提高,代价就是空间分辨率降低。)
了解完相干值和多视视数,应知道干涉相位密度函数
p
d
f
(
ϕ
)
pdf(\phi)
pdf(ϕ)是相干值
γ
\gamma
γ和多视视数
L
L
L的函数。(具体公式这里不写)也就是说,相干值和多视视数可以反应干涉相位的噪声水平。在大多视视数和高相干的情况下,相位误差的近似值为:
σ
2
=
1
2
L
1
−
γ
2
γ
2
360
°
2
π
\sigma^2=\frac{1}{2L}\frac{1-\gamma^2}{\gamma^2}\frac{360°}{2\pi}
σ2=2L1γ21−γ22π360°
因此高相干值或增加多视视数可以取得较小的相位噪声。
影响相干系数的因素
热噪声、时间失相干、空间失相干、体失相干、数据处理的失相干、频谱偏移引起的失相干
-
体失相干
某些地物垂直结构存在的散射差异,比如冰川,森林等。导致同一像素获取的信号不包含相同散射特征,从而引起的失相干。 -
数据处理引起的失相干
配准精度低、重采样 -
热噪声
信号接受存在噪声,即信噪比。每次成像噪声失相干的贡献: γ s n r = 1 1 + N S \gamma_{snr}=\frac{1}{\sqrt{1+\frac{N}{S}}} γsnr=1+SN1 N N N是噪声功率, S S S是信号功率。 -
时间失相干
在两幅SAR获取期间,地面目标散射相位的变化和大气引起的相位误差,导致相干性的降低。- 人工翻耕
- 人工建筑
- 大气状态变化
- 季节导致的植被、水域变化
- 长波段可以减弱地距方向的随机运动误差
L波段下,相干值的时序变化可以用指数模型拟合。 γ ( t ) = ( γ 0 − γ k ) e − t τ + γ k \gamma(t)=(\gamma_0- \gamma_k)e^{-\frac{t}{\tau}}+\gamma_k γ(t)=(γ0−γk)e−τt+γk, γ 0 \gamma_0 γ0为起始相干值, γ k \gamma_k γk为结束相干值, τ \tau τ为时间常数。
- 长波段可以减弱地距方向的随机运动误差
-
几何失相干
平地效应在入射角越小的情况下,条纹频率越大,条纹间距越小。条纹频率为 Δ ϕ 2 π Δ r = 2 B ⊥ λ R tan θ \frac{\Delta\phi}{2\pi\Delta r}=\frac{2B_\bot}{\lambda R\tan\theta} 2πΔrΔϕ=λRtanθ2B⊥
当入射角过小,那意味着条纹间距极小,条纹就没有意义了。或者是当平地效应需要考虑坡度 β \beta β的情况下,条纹间距也为变得很小,即 β = θ \beta=\theta β=θ, Δ ϕ 2 π Δ r = 2 B ⊥ λ R tan ( θ − β ) \frac{\Delta\phi}{2\pi\Delta r}=\frac{2B_\bot}{\lambda R\tan(\theta-\beta)} 2πΔrΔϕ=λRtan(θ−β)2B⊥ -
频谱偏移引起的失相干
即同一反射频谱,因入射角不一致,导致接受频谱有偏移所引起的失相干
- 想要引起频谱偏移导致失相干,那么入射角差异就要足够大,也意味着空间基线足够长。也就意味着存在着不导致失相干的极限基线 公式不给出,给结论
- 坡度角为正,极限基线短,坡度角为负,极限基线长。
- 想要引起频谱偏移导致失相干,那么入射角差异就要足够大,也意味着空间基线足够长。也就意味着存在着不导致失相干的极限基线 公式不给出,给结论
图像配准
掌握配准(registration)的最大谱方法和相关系数法。
InSAR影像配准方法:1.基于强度互相关的配准。 2.基于干涉条纹清晰度配准方法。
InSAR影像配准的大致流程是先采用轨道信息或图像强度信息进行粗配,然后进一步依据强度信息或者干涉条纹的清晰度来进行精配准,同理采样多个点,利用最小二乘法确定配准多项式,最后根据多项式进行重采样。
配准要求:要使干涉图的相关性不低于5%,配准精度必须优于0.2个像元。要求很苛刻。
基于卫星轨道参数粗配准方法
基于卫星轨道参数的粗配准方法认为主辅影像之间存在偏差offset,即:
P
s
u
p
(
l
,
p
)
=
P
m
a
i
n
(
l
,
p
)
+
o
f
f
s
e
t
(
l
,
p
)
P_{sup}(l,p)=P_{main}(l,p)+offset(l,p)
Psup(l,p)=Pmain(l,p)+offset(l,p)
借助多普勒方程和斜距方程计算偏移值。
- 1.计算主影像中心像元在轨道坐标系中的位置
- 2.基于多普勒方程计算该点在辅影像上的行列值
- 3.计算偏移量
互相关系数以及最大谱精配准方法
- 互相关系数法的基本思想是:通过寻找两幅图像互相关系数的最大值来配准,当同一目标的SAR复图像之间精确配准时,其互相关系数取得最大值。
- 最大谱法的基本思想是:当两幅复图像配准精度越高时,图像之间的干涉条纹质量越好,其频谱值最大,找到频谱最大值,就可以找到两幅复图像配准的位置。
- 各点干涉相位可以用两复数共轭相乘表示
c
=
c
1
×
c
2
∗
c=c_1\times c_2^*
c=c1×c2∗对干涉图像各点干涉复数c做傅里叶变换得到它的二维频谱
f
(
c
)
f(c)
f(c)。进而找其最大值。
- 各点干涉相位可以用两复数共轭相乘表示
c
=
c
1
×
c
2
∗
c=c_1\times c_2^*
c=c1×c2∗对干涉图像各点干涉复数c做傅里叶变换得到它的二维频谱
f
(
c
)
f(c)
f(c)。进而找其最大值。
几何变换模型确定多项式及重采样
- 确定多项式
多项式按照阶数N不同可以写为: x ′ = ∑ i = 0 N ∑ j = i N a i j x i y j i x^{'}=\sum\limits_{i=0}^N\sum\limits_{j=i}^Na_{ij}x^iy^{ji} x′=i=0∑Nj=i∑Naijxiyji y ′ = ∑ i = 0 N ∑ j = i N b i j x i y j − i y^{'}=\sum\limits_{i=0}^N\sum\limits_{j=i}^Nb_{ij}x^iy^{j-i} y′=i=0∑Nj=i∑Nbijxiyj−i
- 重采样
重采样的方法可以使用:双线性插值 或 三次卷积插值。
基线计算
基线对于相位的计算至关重要,实际InSAR技术使用过程中,基线难免出现精度不高或者不知道的情况。那么就需要学习基线估计的方法
基于轨道信息的基线估算方法
当卫星星历数据比较精确,地面起伏大,干涉图相干性差的情况下,可以采用轨道信息进行基线估计。主要思想是采用星历数据中的卫星天线中心的瞬时卫星坐标
(
X
s
,
Y
s
,
Z
s
)
(X_s,Y_s,Z_s)
(Xs,Ys,Zs),和瞬时速度向量
(
V
x
s
,
V
y
s
,
V
z
s
)
(V_{xs},V_{ys},V_{zs})
(Vxs,Vys,Vzs),以及对应的时刻信息
t
t
t进行多项式估计。
[
X
s
Y
s
Z
s
]
=
[
a
0
+
a
1
t
+
a
2
t
2
+
a
3
t
3
b
0
+
b
1
t
+
b
2
t
2
+
b
3
t
3
c
0
+
c
1
t
+
c
2
t
2
+
c
3
t
3
]
\begin{bmatrix}X_s\\Y_s\\Z_s \end{bmatrix} =\begin{bmatrix}a_0+a_1t+a_2t^2+a_3t^3\\b_0+b_1t+b_2t^2+b_3t^3\\c_0+c_1t+c_2t^2+c_3t^3 \end{bmatrix}
XsYsZs
=
a0+a1t+a2t2+a3t3b0+b1t+b2t2+b3t3c0+c1t+c2t2+c3t3
[
V
x
s
V
y
s
V
z
s
]
=
[
a
1
+
2
a
2
t
+
3
a
3
t
2
b
1
+
2
b
2
t
+
3
b
3
t
2
c
1
+
2
c
2
t
+
3
c
3
t
2
]
\begin{bmatrix}V_{xs}\\V_{ys}\\V_{zs} \end{bmatrix} =\begin{bmatrix}a_1+2a_2t+3a_3t^2\\b_1+2b_2t+3b_3t^2\\c_1+2c_2t+3c_3t^2 \end{bmatrix}
VxsVysVzs
=
a1+2a2t+3a3t2b1+2b2t+3b3t2c1+2c2t+3c3t2
基于干涉图信息的基线估算方法
当卫星轨道信息不太精确,但干涉区域有较大地形平坦区域,且干涉图相干性好时,可采用这种方法。其好处在于这是一种盲估计方法,只需要干涉相位,干涉条纹等信息,不需要地面控制点。
B
⊥
=
r
1
tan
(
θ
−
α
)
f
0
f
f
r
i
n
g
e
B_\bot=\frac{r_1\tan(\theta-\alpha)}{f_0}f_{fringe}
B⊥=f0r1tan(θ−α)ffringe
f
0
f_0
f0是载波频率,
f
f
r
i
n
g
e
f_{fringe}
ffringe是条纹频率。
基于地面控制点的基线估计方法
若地面在同一距离向上有两个以上控制点时,可以通过干涉相位和几何地面几何关系来求解基线长度
B
B
B和基线水平角
α
\alpha
α。
InSAR的相位滤波和解缠
滤波
该部分主要了解InSAR干涉条纹图的评价方法和主要滤波方法
-
评价方法
- 相干值图像
略 - 残差点数
解缠部分阐述 - 相位标准偏差指标(PSD)
即整幅图像或局部窗口图像的各位相位值的标准差。这个值越小,说明相位噪声越小,干涉图质量越好。
σ ϕ = ∑ N ( ϕ ( i , j ) − ϕ ˉ ( i , j ) ) 2 N − 1 \sigma_\phi=\sqrt{\frac{\sum\limits^N(\phi(i,j)-\bar\phi(i,j))^2}{N-1}} σϕ=N−1∑N(ϕ(i,j)−ϕˉ(i,j))2 - 相位梯度的和值指标(SPD)
略,也是值越小越好。
- 相干值图像
-
滤波方法
- 空间域滤波:均值滤波和中值滤波,条纹自适应滤波,多视滤波(容易造成边缘图像缺失)。
- 频率域滤波:Goldstein滤波(基于功率谱的卷积滤波方法)。
- 时频分析滤波:矢量分离式小滤波,零中频矢量滤波
解缠
1.了解什么是解缠,一维解缠的难点。
2.了解一维解缠到二维解缠的推广,二维解缠需要兼顾的两个方面。
3.了解InSAR领域中相位解缠
SAR影像接受到的相位是一个处在
(
−
π
,
π
]
(-\pi,\pi]
(−π,π],表现在相位图上就有突变。
我们将这种现象描述为, ϕ \phi ϕ是一个缠绕在 ( − π , π ] (-\pi,\pi] (−π,π]之间的非线性信号。所以将相位由主值或相位差值恢复为真实值的过程统称为相位解缠。
-
解缠条件:采样率满足 Nyquist 采样定理,采样频率必须大于信号最高频率的两倍,解缠绕的干涉相位中相邻像素点之间的相位差值不可能超过半个周期(一个π)。
- 采样数??????????????????????? -
解缠步骤:
- 从数据序列左边第二个开始,求得相位差。
- 如果相位差大于 π \pi π,而从第二个数据开始,减去 − 2 π -2\pi −2π
- 如果相位差小于 π \pi π,而从第二个数据开始,加上 − 2 π -2\pi −2π
- 接着判断第三和第二个数据的相位差,如果相位差绝对值大于 π \pi π,执行前两步,反之保留该值。直到所有数据完成。
-
~~一维解缠的难点:噪声水平和采样数~
-
一维解缠推广到二维解缠:首先提取垂直向和水平向的相位偏导数,然后沿着垂直方向和水平方向积分,就可以达到相位解缠的目的。
-
二维解缠需要兼顾的两个方面
- 一致性:即解缠后,任意两点间的相位差与路径无关。
- 精确性:即解缠后的相位要忠实地恢复原始相位。
-
InSAR领域中的相位解缠
- 路径跟踪法
- 路径跟踪方法就是通过积分(积分可以类比加法)相邻缠绕相位的差分值来恢复相位的真实值的。
即如果已知像元 r 0 r_0 r0上的相位,那么在其它像元 r r r上的相位可以表示为 ϕ ( r ) = ∫ c ∇ ϕ d r + ϕ ( r 0 ) \phi(r)=\int_c\nabla\phi dr+\phi(r_0) ϕ(r)=∫c∇ϕdr+ϕ(r0)( ∇ \nabla ∇是微分符号。) - 根据一致性,我们知道相位解缠后要保证两点相位差与路径无关,且在高数中我们知道,积分要与路径无关,需要保证闭合积分成立。
人话说就是路径要在单连通区域,而一旦区域内有洞,就变成复连通区域了 - 在实际的缠绕相位数据中,常常因为受到噪声的影响导致像元相位闭合路径不满足条件(即有洞)。这样的像元点就成为残差点。在路径跟踪的相位解缠算法中。而过能找到这样的残差点,并通过正负相连(将复连通区域中的洞连起来,即分割成了多个单连通区域)从而做到防止积分路径穿过这些连线(“分枝”)。将残差点连起来的方法也称为枝切法。分枝将图像分成各个区域后,再进行二维解缠。
- 路径跟踪方法就是通过积分(积分可以类比加法)相邻缠绕相位的差分值来恢复相位的真实值的。
- 还有Goldstein枝切法,最小范数法,最优估计算法。
- Goldstein枝切法:优点是运算速度快,噪声低,残差点少的情况,精确度高。缺点是一旦残差点多就难以准确链接枝切线。
- 最小范数法:主要是通过数学方法约束,使得解缠后的相位最接近真实值,是一种全局解缠方法。
- 路径跟踪法
时序InSAR
1.掌握两种方法的原理和优劣
InSAR技术常常用来关注形变。那么形变量,形变速率就是关注的具体指标了。因此时序InSAR技术逐渐发展。
小基线集(SBAS)
-
Stacking技术
SAR影像中,相位干涉模型可以表示为:
ϕ = − 4 π λ B ⊥ h R sin θ + 4 π λ Δ r + 4 π λ S t 1 − 4 π λ S t 2 + 2 π w k \phi=-\frac{4\pi}{\lambda}\frac{B_\bot h}{R\sin\theta}+\frac{4\pi}{\lambda}\Delta r+\frac{4\pi}{\lambda}S^ {t_1}-\frac{4\pi}{\lambda}S^ {t_2}+2\pi w_k ϕ=−λ4πRsinθB⊥h+λ4πΔr+λ4πSt1−λ4πSt2+2πwk
在一个长时间跨度下,同一区域获取 N N N张SAR影像组成共 M M M张干涉影像集(stack)。其中的高程相位,大气延迟相位,噪声相位,(参考面相位直接去除)在长时间跨度下可以视作随机噪声。即噪声符合正态分布,期望为零。且将区域形变视为线性形变: d ( x ) = v ( x ) × t d(x)=v(x)\times t d(x)=v(x)×t那么形变就容易求得为: v ˉ = − λ 4 π ∑ i = 1 M ϕ i ∑ i = 1 M t i \bar v=-\frac{\lambda}{4\pi}\frac{\sum\limits_{i=1}^{M}\phi_i}{\sum\limits_{i=1}^{M}t_i} vˉ=−4πλi=1∑Mtii=1∑Mϕi
对多幅干涉图像各个点位的相位写成矩阵形式做最小二乘求解,从而求得图像每个点的形变速率。
这种方法可以求得的形变速率,优点是模型简单。缺点也很明显:- 1.干涉影像要求时间基线选择足够合适,使得高程残差相位,大气以及噪声相位的统计分布符合高斯分布。因此基于SAR影像选取干涉对要足够细心。
- 2.形变模型为线性形变,不能计算类似雪山,火山地形这种非线性形变。
- 3.求得的形变速率是一个全局值,无法求得监测区域的短时形变速率或瞬时形变速率。
- 4.干涉对的空间基线要短,从而降低高程带来的相位影响。
为了解决上述缺点,引出了小基线集技术
- 小基线集(SBAS)
在Stacking技术的基础上,不再将大气延迟,噪声和轨道等因素考虑为随机噪声,然后选择时间基线相对短的SAR影像组成干涉对。
每两张SAR影像之间的形变速率设为 v i v_i vi
从而可以将干涉相位如下表示:
ϕ d e f o k = − 4 π λ Δ r x , y k = − 4 π λ ∑ ( t j − t j − 1 ) v j = β k V \phi_{defo}^k=-\frac{4\pi}{\lambda}\Delta r_{x,y}^k=-\frac{4\pi}{\lambda}\sum(t_j-t_{j-1})v_j=\beta_kV ϕdefok=−λ4πΔrx,yk=−λ4π∑(tj−tj−1)vj=βkV ϕ t o p o k = − 4 π λ B ⊥ , x , y r x , y k sin θ x , y k Δ h x , y = α x , y k Δ h x , y \phi_{topo}^k=-\frac{4\pi}{\lambda}\frac{B_{\bot,x,y}}{r^k_{x,y}\sin\theta^k_{x,y}}\Delta h_{x,y}=\alpha^k_{x,y}\Delta h_{x,y} ϕtopok=−λ4πrx,yksinθx,ykB⊥,x,yΔhx,y=αx,ykΔhx,y那么,总的相位就可以表示为 Δ Φ = α Δ h + β V + w k \Delta \Phi=\alpha\Delta h+\beta V+w_k ΔΦ=αΔh+βV+wk其中, w k = ϕ o r b i t + ϕ a t m o + ϕ n o i s e w_k=\phi_{orbit}+\phi_{atmo}+\phi_{noise} wk=ϕorbit+ϕatmo+ϕnoise
因此,利用上述式子,列成矩阵求解,可以解出该点的短时形变速率。- 上述考虑的情况都是基线连通和各个干涉相位数据相干性都较好的的情况,当出现基线不连通的情况时。
首先可以考虑增加数据来源。利用先验数据。而后考虑采用添加虚拟观测,即认为速率呈线性变化: Δ v k = 0.5 Δ v k − 1 + 0.5 Δ v k + 1 \Delta v_k=0.5\Delta v_{k-1}+0.5\Delta v_{k+1} Δvk=0.5Δvk−1+0.5Δvk+1其矩阵形式如下:
- 当遇到某些基线上的点相干性差的情况,可以考虑去掉该观测值。那么就会出现相干点间断的情况,如图中的
B
B
B:
称这样的技术为间断相干小基线集技术 - 如果解算过程中出现了秩亏的情况,可以采用奇异值分解的方法求得所需的线性形变速率解。
- 上述考虑的情况都是基线连通和各个干涉相位数据相干性都较好的的情况,当出现基线不连通的情况时。
PSInSAR
PSInSAR即在SAR影像集中,选择一幅主影像,依据主影像与其他影像组成干涉对形成干涉影像,进而解求形变速率。
PS(Persistent Scatter)稳定散射回波点,也意味着这种方法,需要寻找良好相干性的点。比如,建筑物,桥梁等人工目标或者裸露土地,岩体等均一天然目标。
- 选择PS点的方法
- 时间序列相关系数阈值法
其实就是比较相干值,只不过干涉影像时间基线不为零,所以就是时序相关系数阈值法。
- 振幅离差指数阈值法
在前文讲相关系数中有提到,在相位标准差较小的情况下,振幅离差和相位标准差相差较小,可以作为相干性的指标。需要注意的是,不同时间SAR成像获取的振幅值不同,所以这种方法的使用需要提前定标。
- 时间序列相关系数阈值法
- 主影像的选择方法
主影像选择总体相干性高的影像: ρ t o t a l ≈ ( 1 − f ( T i T c ) ) α ( 1 − f ( B ⊥ i B ⊥ c ) ) β ( 1 − f ( F D C i F D C c ) ) γ \rho_{total}\approx(1-f(\frac{T^i}{T^c}))^\alpha(1-f(\frac{B_\bot^i}{B_\bot^c}))^\beta(1-f(\frac{F_{DC}^i}{F_{DC}^c}))^\gamma ρtotal≈(1−f(TcTi))α(1−f(B⊥cB⊥i))β(1−f(FDCcFDCi))γ α + β + γ = 1 \alpha+\beta+\gamma=1 α+β+γ=1
依据上面的式子,选择不同的 α \alpha α, β \beta β, γ \gamma γ,总体相关性如图:
以上图为例,应当选择第9幅为主影像。
注:SBAS相比于PSInSAR对于失相干的处理更好,PSInSAR因为选出来的都是相干性好的点,可以考虑使用网平差的思路求解高程改正和形变增量。
地基干涉雷达
1.了解星载与地基的区别及其各自优劣
2.了解地基原理
3.了解地基雷达误差源及其分析
地基干涉雷达原理
回顾之前的星载雷达或者机载雷达,它们采用的监测波一般是都是脉冲波,其距离分辨率的大小差异,主要取决于雷达所能分辨的回波间隔 τ \tau τ,即脉冲宽度。
脉冲宽度和频带宽度:前者是时域概念,是指脉冲波一个周期内的两个零点间的时间间隔。而后者是频域概念,在周期信号频谱中,从零频率到需要考虑的最高次谐波频率之间的频段即为该信号的有效占有带宽,亦称频带宽度 B B B。也可以简单理解为最低频率到最高频率的宽度。两者大致关系为: Δ R = c τ 2 = c 2 B \Delta R=\frac{c\tau}{2}=\frac{c}{2B} ΔR=2cτ=2Bc
脉冲雷达
- 优点:
- 具有较高的功率峰值,可作用距离远
- 宽带信号有利于提高系统的抗干扰能力
- 缺点:
- 制作繁琐
- 雷达结构复杂,重量体积大
地基雷达
检测波采用的步进频率连续波:雷达发送的是n组连续频率的电磁波。
在频域上级显示为几个峰值
可以理解为相比脉冲波,增加了带宽,从而达到了提高距离向分辨率的目的。
- 优点:
- 能检测很近的距离,测量精度较高
- 简单、小巧、轻便
- 带宽大,平均发射功率低,被截获的概率小
- 缺点:
- 发射功率小,距离近,收发间难完成隔离
误差源及其改正方法
- 轨道位移误差
按照星载SAR的干涉相位模型,相位分为参考面相位,高程相位,形变相位,大气延迟以及噪声相位 ϕ = − 4 π λ B ∥ − 4 π λ B ⊥ h R sin θ + 4 π λ Δ r + 4 π λ S t 1 − 4 π λ S t 2 + 2 π w k \phi=-\frac{4\pi}{\lambda}B_\parallel-\frac{4\pi}{\lambda}\frac{B_\bot h}{R\sin\theta}+\frac{4\pi}{\lambda}\Delta r+\frac{4\pi}{\lambda}S^ {t_1}-\frac{4\pi}{\lambda}S^ {t_2}+2\pi w_k ϕ=−λ4πB∥−λ4πRsinθB⊥h+λ4πΔr+λ4πSt1−λ4πSt2+2πwk然而地基雷达理论上空间基线为0,因此参考面相位和高程相位为0,因此地基干涉雷达相位模型为 ϕ = 4 π λ Δ r + ϕ a t m + ϕ n o i s e \phi=\frac{4\pi}{\lambda}\Delta r+\phi_{atm}+\phi_{noise} ϕ=λ4πΔr+ϕatm+ϕnoise
在实际应用时,监测站长时间监测难免雷达中心有有所移动。因此,最终相位中还包含着轨道误差
ϕ = 4 π λ Δ r + ϕ B + ϕ a t m + ϕ n o i s e \phi=\frac{4\pi}{\lambda}\Delta r+\boxed{\phi_B}+\phi_{atm}+\phi_{noise} ϕ=λ4πΔr+ϕB+ϕatm+ϕnoise
将轨道误差按照距离向,方位向,竖直方向上做分解- 距离向轨道误差:距离向误差与方位向移动量和方位角相关,且影响较大。当方位向角度为20°时,SAR中心在距离向移动1.06mm,可引起1mm的形变误差。
- 方位向轨道误差:方位向误差与方位向移动量和方位角相关。当方位向角度为20°时,SAR中心在方位向移动2.92mm,可引起1mm的形变误差。
- 竖直方向轨道误差:竖直方向误差与竖直方向移动量,目标点距离和高程相关。当目标点距离为200m,高程为50m,竖直向移动距离为4mm时,可引起1mm形变误差。
- 大气误差
- 可以采用函数模型 Δ r = a 1 + a 2 R + a 3 R 2 \Delta r=a_1+a_2R+a_3R^2 Δr=a1+a2R+a3R2改正。或者采用经验模型改正,公式略。
- 上述函数模型是一维情况,以斜距作为参考量来对大气延迟建模。还可以加入方位角和高程建立二维大气模型甚至三维大气模型。
- 无论哪种函数模型,其考虑的大气状况都相对单一。然后在实际状况中,不同地物区域,甚至不同时间段,产生的大气误差都不尽相同。大气延迟模型复杂是常见状况。因此可以采用地物聚类——分类建模的方法或者其他方法加以改正。
地质灾害
1.崩塌,滑坡,泥石流,采用塌陷
滑坡识别:老滑坡看长相,新滑坡看迹象