劉 飛，潘 斌

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

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

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

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

1 連續(xù)干涉像對(duì)與PS點(diǎn)選取

1.1 連續(xù)干涉像對(duì)

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

ψi=φi-φi+1

(1)

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

1.2 配準(zhǔn)主影像選擇

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

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

(2)

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

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

1.3 PS點(diǎn)選取

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

連續(xù)干涉像對(duì)和單主影像干涉像對(duì)在PS點(diǎn)的選取算法上大致相同。但由于連續(xù)干涉像對(duì)的相干性提高，生成的干涉圖質(zhì)量更優(yōu)，使得穩(wěn)定PS點(diǎn)的噪聲降低，相位殘差更小，PS點(diǎn)的選取更精準(zhǔn)。同時(shí)單主影像干涉像對(duì)方法會(huì)引入由于主影像存在于所有干涉圖上而產(chǎn)生的主影像誤差，解算中需要額外求解此誤差項(xiàng)，而連續(xù)干涉像對(duì)方法則避免了這一誤差項(xiàng)，提高其解算精度。

2 形變相位解算

2.1 相位組成分析

在對(duì)選取的PS點(diǎn)進(jìn)行三維相位解纏[14]后，對(duì)于給定像素點(diǎn)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依次對(duì)應(yīng)形變相位、大氣相位、軌道誤差相位、地形誤差項(xiàng)相位及噪聲；2kx,iπ為解纏相位的模糊度，對(duì)于同一干涉圖上不同PS點(diǎn)而言，在相位解纏精度足夠高的情況下，其數(shù)值相同[15]。

形變相位ψD,x,i既空間相關(guān)，又時(shí)間相關(guān)。假設(shè)像素點(diǎn)的相位值可以用時(shí)間上的三次函數(shù)進(jìn)行擬合

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

(4)

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

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

ti=Δt·(i-1)

(5)

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

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

(6)

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

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

(7)

式中，θx為需要求解的未知參數(shù)，對(duì)于一個(gè)PS點(diǎn)而言，θx在所有干涉圖上的值均相同。

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

2.2 地形誤差相位解算

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

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

(8)

式中，Δ表示差分運(yùn)算符。則Δψ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)可知，經(jīng)過干涉圖差分后，線性形變被消除，意味著對(duì)于在時(shí)間上只表現(xiàn)線性形變的點(diǎn)而言，其ΔψD,x,i為零，而對(duì)于含有非線性形變的點(diǎn)而言，在去除線性形變后，并隨著方程階數(shù)的降低，其數(shù)值也通常較小。

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

(10)

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

Aα=δψ

(11)

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

(12)

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

(13)

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

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

2.3 形變相位解算

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

(14)

大氣相位的存在仍會(huì)對(duì)形變解算精度產(chǎn)生較大的影響，為了消去大氣相位，通過與高斯函數(shù)進(jìn)行卷積運(yùn)算，對(duì)相位進(jìn)行時(shí)間上的低通濾波，則

(15)

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

Bv=ψ

(16)

式中，B為N×3的矩陣

(17)

v為3×K的矩陣，K表示PS點(diǎn)個(gè)數(shù)

(18)

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

3 試驗(yàn)結(jié)果與分析

3.1 連續(xù)干涉像對(duì)

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

表1 連續(xù)干涉像對(duì)

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

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

圖1 像素平均相干性對(duì)比

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

圖2 PS點(diǎn)相位殘差

3.2 地形誤差相位解算

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

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

圖3 地形誤差相位解算

在解算出地形誤差相位后，可在去除地形誤差相位的基礎(chǔ)上重新進(jìn)行相位解纏，以提高解纏精度，再用改善的解纏結(jié)果重新解算地形誤差相位。此步驟可重復(fù)數(shù)次直至地形誤差相位解算結(jié)果不再改變。

3.3 試驗(yàn)結(jié)果驗(yàn)證

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

圖4 平均沉降速率

圖4中采用同時(shí)期的Landsat 8灰度影像作為底圖，測(cè)量出沉降值的區(qū)域即為選取PS點(diǎn)的位置，標(biāo)出的十字點(diǎn)為用于驗(yàn)證結(jié)果的水準(zhǔn)點(diǎn)。由圖4可知，PS點(diǎn)準(zhǔn)確地還原了天津市的建筑、橋梁、道路及城鎮(zhèn)的位置，測(cè)量出沉降漏斗所在的位置，清晰地顯示了地面沉降在空間上的相關(guān)性。

為了驗(yàn)證測(cè)量結(jié)果的可靠性，將測(cè)量的結(jié)果與水準(zhǔn)數(shù)據(jù)進(jìn)行對(duì)比，水準(zhǔn)點(diǎn)的位置如圖4所示。本文采用克里金插值得到水準(zhǔn)點(diǎn)測(cè)量值，其對(duì)比結(jié)果見表2。

表2 試驗(yàn)結(jié)果對(duì)比 mm

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

4 結(jié) 語(yǔ)

本文研究結(jié)合Sentinel-1數(shù)據(jù)的特點(diǎn)，提出了基于連續(xù)干涉像對(duì)的PSInSAR算法。相比傳統(tǒng)單主影像干涉組合方法，新算法能夠使時(shí)間導(dǎo)致的去相干最小化，獲取更高的相干性，從理論上整體提高PSInSAR算法的精度。同時(shí)利用Sentinel-1數(shù)據(jù)具有重訪周期固定的特點(diǎn)，通過干涉圖差分去除線性形變，使用三角網(wǎng)平差單獨(dú)求解地形誤差相位，達(dá)到避免使用空間濾波、提高解算精度的目的。最終通過監(jiān)測(cè)天津市的地表沉降，將測(cè)量結(jié)果與水準(zhǔn)數(shù)據(jù)進(jìn)行對(duì)比，證明了新算法在該地區(qū)的可行性，具有后續(xù)的研究?jī)r(jià)值。

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