一種基于Sentinel-1連續干涉像對的PSInSAR算法

劉 飛，潘 斌

(武漢大學遙感信息工程學院，湖北 武漢 430079)

地面沉降是一種緩慢變化的地質災害，可引起城市建筑物傾斜，破壞地基的穩定性，給生產和生活帶來很大影響。天津市在過去幾十年中由于地下水開采過度遭受了嚴重的地表沉降災害，研究如何有效地監測該地表沉降具有重要的理論價值和實際意義。

永久散射體差分干涉技術(PSInSAR)是近十幾年來新興的遙感技術手段。該技術具有對地表進行長時間、大面積、高精度形變監測的能力，是近些年來沉降監測的有效手段[1-4]。PSInSAR技術應用于多顆星載SAR衛星(如ERS1/2、Envisat、TerraSAR-X等)的數據處理上，在對地表形變監測的應用中可取得與地表數據相符合的測量精度。自2014年和2016年Sentinel-1A和Sentinel-1B發射成功后，其數據由于共享性，且具有更先進的成像方式、更短的重訪周期及更精密的軌道控制等優點[5]，而成為新的研究熱點[6-10]。

針對以往PSInSAR算法受限于數據條件，采用單主影像干涉組合方法而導致的低相干性、解算流程復雜等問題，本文研究針對Sentinel-1的數據特點，提出連續干涉像對的組合方法，利用相鄰成像時間的2幅SAR影像生成干涉圖，使時間基線導致的去相干最小化，提高相干性。在形變解算流程中，利用重訪周期固定的特點，通過干涉圖差分消去線性形變，在不影響形變解算的情況下，單獨解算地形誤差相位。本文以天津市為研究區，基于2016年11月至2018年1月32幅Sentinel-1B數據，采用新算法生成該地區的沉降速率圖，并通過與水準數據比較驗證新算法的有效性。

1 連續干涉像對與PS點選取

1.1 連續干涉像對

區別于傳統PSInSAR單主影像干涉組合方法，連續干涉像對按照成像時間順序，選擇成像時間相鄰的2幅影像兩兩生成干涉像對，即對于按成像時間獲取的每幅影像而言，其既是下一幅影像的主影像，又是上一幅影像的輔影像。對于由N+1幅SAR影像生成的N幅連續干涉像對而言，其第i張干涉相位ψi有

ψi=φi-φi+1

(1)

式中，i=1,2,…,N；φi、φi+1分別表示第i張和第i+1張影像的相位值。

1.2 配準主影像選擇

盡管對于任意一幅SAR影像而言，其選擇時間上對應的上一幅影像作為主影像，但仍然需要從所有影像中選擇一幅配準主影像，將其他所有影像與之配準，以滿足后續的干涉圖生成。配準主影像的選擇需要使得總的配準去相干噪聲最小化，其選擇方式與傳統PSInSAR選擇主影像相同，配準去相干噪聲ρtotal表達式為[11]

ρtotal=ρtemporal·ρspatial·ρdoppler·ρthermal

(2)

式中，ρtemporal、ρspatial、ρdoppler和ρthermal分別表示由時間基線、空間基線、多普勒中心頻移及熱噪聲導致的去相干。

盡管理論上在選擇配準主影像時要考慮上述所有因素，但對于Sentinel-1數據而言，由于具有更精密的軌道控制，無論選擇哪幅影像作為配準主影像，空間基線都保持在一個較小的范圍內，同時其數據的多普勒中心頻移均很小，而熱噪聲在這個過程中被認為是常量，因此，時間基線便成了關鍵因素。為了盡可能減少由時間基線導致的去相干，通常選擇靠近成像時間中心的影像作為配準主影像。

1.3 PS點選取

在配準生成連續干涉像對后，需要從干涉像對中提取相位穩定且相干性較高的點，即PS點，用于后續的相位解纏和形變解算。本文采用基于幅度離差[12]的方法先進行PS點的初步篩選，再根據PS點的相位特征[13]進行精選。

連續干涉像對和單主影像干涉像對在PS點的選取算法上大致相同。但由于連續干涉像對的相干性提高，生成的干涉圖質量更優，使得穩定PS點的噪聲降低，相位殘差更小，PS點的選取更精準。同時單主影像干涉像對方法會引入由于主影像存在于所有干涉圖上而產生的主影像誤差，解算中需要額外求解此誤差項，而連續干涉像對方法則避免了這一誤差項，提高其解算精度。

2 形變相位解算

2.1 相位組成分析

在對選取的PS點進行三維相位解纏[14]后，對于給定像素點x，其在第i張干涉圖上的解纏相位ψx,i為

ψx,i=ψD,x,i+ψA,x,i+ψS,x,i+ψθ,x,i+ψN,x,i+2kx,iπ

(3)

式中，ψD,x,i、ψA,x,i、ψS,x,i、ψθ,x,i和ψN,x,i依次對應形變相位、大氣相位、軌道誤差相位、地形誤差項相位及噪聲；2kx,iπ為解纏相位的模糊度，對于同一干涉圖上不同PS點而言，在相位解纏精度足夠高的情況下，其數值相同[15]。

形變相位ψD,x,i既空間相關，又時間相關。假設像素點的相位值可以用時間上的三次函數進行擬合

ψx(t)=a0,x+a1,xt+a2,xt2+a3,xt3

(4)

式中，a0,x和a1,x為常量和平均速率值，代表線性形變；a2,x和a3,x為加速度和加速度變化值，代表非線性形變。

由于干涉圖連續生成，且重訪周期固定，給定第一張SAR影像的成像時間t1=0，則第i張影像的成像時間ti為

ti=Δt·(i-1)

(5)

式中，i=1,2,…,N+1；Δt為重訪周期，即時間基線。則ψD,x,i為

ψD,x,i=ψx(ti)-ψx(ti+1)=

(6)

地形誤差相位ψθ,x,i是部分空間相關的，其大小與垂直基線B⊥,i成正比，可寫為

ψθ,x,i=B⊥,i·θx

(7)

式中，θx為需要求解的未知參數，對于一個PS點而言，θx在所有干涉圖上的值均相同。

大氣相位ψA,x,i是相鄰兩次成像大氣延遲相位之差，為空間相關但時間不相關。軌道誤差相位ψS,x,i可通過Deramp處理單獨估算并去除。而2kiπ則可通過選取參考點，在所有干涉圖上減去參考點對應的相位值的方法而消除。

2.2 地形誤差相位解算

為了消除地形誤差相位對形變解算的影響，本文利用重訪周期固定的特點，在干涉圖之間作差分處理，在單獨去除軌道誤差相位后，由式(3)得

Δψx,i=ΔψD,x,i+Δψθ,x,i+ΔψA,x,i+ΔψN,x,i+2Δkx,iπ

(8)

式中，Δ表示差分運算符。則ΔψD,x,i為

ΔψD,x,i=ψD,x,i-ψD,x,i+1=

[ψx(ti)-ψx(ti+1)]-[ψx(ti+1)-ψx(ti+2)]=

2Δt2·a2,x+6Δt2·ti+1·a3,x

(9)

由式(9)可知，經過干涉圖差分后，線性形變被消除，意味著對于在時間上只表現線性形變的點而言，其ΔψD,x,i為零，而對于含有非線性形變的點而言，在去除線性形變后，并隨著方程階數的降低，其數值也通常較小。

為了進一步消除其他項對地形誤差解算的影響，在式(8)的基礎上，對相鄰PS點之間進行差分，消去2Δkx,iπ，并代入式(7)和式(9)

(10)

式(10)可改寫為矩陣形式

Aα=δψ

(11)

式中，A為(N-1)×3的矩陣

(12)

α為3×M的矩陣，其中M為式(10)的個數

(13)

δψ為(N-1)×M的矩陣，表示差分干涉圖的相位值。

這種新的地形誤差解算方法利用兩次差分，能避免地形誤差相位解算過程中對形變相位的影響，并且通過三角網平差的方法提高解算精度，在整個解算過程中避免使用空間濾波。

2.3 形變相位解算

在去除地形誤差相位后，由式(3)得

(14)

大氣相位的存在仍會對形變解算精度產生較大的影響，為了消去大氣相位，通過與高斯函數進行卷積運算，對相位進行時間上的低通濾波，則

(15)

將式(15)改寫為矩陣形式，為

Bv=ψ

(16)

式中，B為N×3的矩陣

(17)

v為3×K的矩陣，K表示PS點個數

(18)

式(16)同樣可采用最小二乘方法求解，最終通過乘以系數λ/4π得到PS點在雷達視線(light of sight,LOS)方向的形變量。

3 試驗結果與分析

3.1 連續干涉像對

本文使用的DEM為3″的SRTM數據，SAR影像數據為歐空局(ESA)Sentinel-1B降軌SLC影像，成像時間自2016年11月至2018年1月，共32幅影像。通過連續干涉像對的方法，選擇20170624影像作為配準主影像，生成31幅干涉圖，見表1。

表1 連續干涉像對

表1統計了連續干涉像對的基線及多普勒中心頻移，可以看出連續干涉像對的時間基線大多為12 d，垂直基線均小于200 m，相比于傳統PSInSAR使用單主影像的方法，連續干涉像對方法能在保證垂直基線較短的同時，使時間基線導致的去相干最小化，具有更高的相干性(如圖1所示)。

圖1分別為使用連續干涉像對和單主影像干涉像對方法，對像素點在所有干涉圖上平均相干性的統計。由圖中可知，連續干涉像對的方法有效地提升了相干性，并且消除了單主影像干涉像對中存在像素點完全失相干的情況。

圖1 像素平均相干性對比

相干性的提升也提升了PS選點的質量，圖2為對PS選點過程中求解的平均相位殘差的統計。由圖2可知，在無需估算單主影像干涉像對存在主影像誤差的情況下，絕大部分提取的PS點的相位殘差在-0.1～0.1間波動，其具有較高的質量。

圖2 PS點相位殘差

3.2 地形誤差相位解算

理論上采用連續干涉像對的方法，所有干涉圖均應具有相同的時間基線，但由于多種因素，如雷達天線問題、緊急任務或衛星存儲限制等，會存在預期的影像沒有成像的情況，導致干涉圖差分存在斷點。本文研究從表1中選出3個連續時間段的干涉像對(20170107-20170401，20170425-20170823和20170916-20180102)，生成23張差分干涉圖，對地形誤差相位進行求解。

圖3顯示了單個PS點地形誤差相位由式(11)解算的例子，由圖中可以看出擬合值能較好地還原觀測值的變化趨勢，較精確地解算地形誤差相位，證明了兩次差分方法求解地形誤差相位的可靠性。

圖3 地形誤差相位解算

在解算出地形誤差相位后，可在去除地形誤差相位的基礎上重新進行相位解纏，以提高解纏精度，再用改善的解纏結果重新解算地形誤差相位。此步驟可重復數次直至地形誤差相位解算結果不再改變。

3.3 試驗結果驗證

在由式(16)解算出LOS方向的形變后，通過除以入射角的余弦值將其轉化為豎直方向上的形變，生成天津市的沉降速率圖，如圖4所示。

圖4 平均沉降速率

圖4中采用同時期的Landsat 8灰度影像作為底圖，測量出沉降值的區域即為選取PS點的位置，標出的十字點為用于驗證結果的水準點。由圖4可知，PS點準確地還原了天津市的建筑、橋梁、道路及城鎮的位置，測量出沉降漏斗所在的位置，清晰地顯示了地面沉降在空間上的相關性。

為了驗證測量結果的可靠性，將測量的結果與水準數據進行對比，水準點的位置如圖4所示。本文采用克里金插值得到水準點測量值，其對比結果見表2。

表2 試驗結果對比 mm

表2中各點的真實沉降值分別為2016年11月及2017年10月測量的水準差，為絕對值，而測量值則是選擇SBM26作為參考點(因此該點的相對誤差為零)，由此得出的相對值。由表2可以看出，新算法整體測量精度較高，證明了新算法的可行性。

4 結 語

本文研究結合Sentinel-1數據的特點，提出了基于連續干涉像對的PSInSAR算法。相比傳統單主影像干涉組合方法，新算法能夠使時間導致的去相干最小化，獲取更高的相干性，從理論上整體提高PSInSAR算法的精度。同時利用Sentinel-1數據具有重訪周期固定的特點，通過干涉圖差分去除線性形變，使用三角網平差單獨求解地形誤差相位，達到避免使用空間濾波、提高解算精度的目的。最終通過監測天津市的地表沉降，將測量結果與水準數據進行對比，證明了新算法在該地區的可行性，具有后續的研究價值。

致謝：感謝ESA提供的Sentinel-1數據及開源數據處理軟件SNAP。