本發(fā)明屬于地表沉降監(jiān)測(cè)技術(shù)領(lǐng)域，具體地說(shuō)，涉及一種適用于大梯度地表沉降監(jiān)測(cè)的時(shí)序差分雷達(dá)干涉方法。

背景技術(shù)：

地表沉降(垂直向地表形變，為便于表述，后續(xù)采用“形變”代表“沉降”)是發(fā)生范圍最廣的地質(zhì)災(zāi)害之一，具有持續(xù)時(shí)間長(zhǎng)的特點(diǎn)，且多發(fā)于城市及其近郊等經(jīng)濟(jì)發(fā)達(dá)和人口聚集區(qū)，對(duì)經(jīng)濟(jì)發(fā)展、城市建設(shè)和人民生活均會(huì)構(gòu)成持久危害。我國(guó)是世界上地表沉降災(zāi)害最為嚴(yán)重的國(guó)家之一，累積沉降大于200毫米的面積超過(guò)15萬(wàn)平方公里，主要集中在華北平原、長(zhǎng)江三角洲和汾渭盆地等經(jīng)濟(jì)發(fā)達(dá)地區(qū)，且出現(xiàn)了嚴(yán)重的沉降漏斗，造成了嚴(yán)重的經(jīng)濟(jì)損失。對(duì)地表沉降(尤其是沉降漏斗)開(kāi)展大范圍精確監(jiān)測(cè)，對(duì)沉降防控及避免相應(yīng)的危害具有十分重要的現(xiàn)實(shí)意義。

目前，時(shí)序差分合成孔徑雷達(dá)(Synthetic Aperture Radar,SAR)干涉(Differential SAR Interferometry,DInSAR)技術(shù)已在地表形變監(jiān)測(cè)中得到廣泛應(yīng)用，且是微波遙感、衛(wèi)星大地測(cè)量以及地球物理學(xué)領(lǐng)域研究和應(yīng)用的熱點(diǎn)之一。時(shí)序DInSAR(Time Series DInSAR,TS-DInSAR)具有覆蓋范圍廣、空間分辨率高、效率高、精度高且不易受云霧和陰雨天氣影響的技術(shù)優(yōu)勢(shì)，非常適用于開(kāi)展大范圍地表形變監(jiān)測(cè)。TS-DInSAR對(duì)緩慢累積性地表形變具有較好的監(jiān)測(cè)效果和精度，但當(dāng)形變較快和形變梯度較大時(shí)易造成解算精度較低或解算失敗(如形變欠估計(jì)和形變漏斗區(qū)的“空洞”現(xiàn)象，即結(jié)果缺失)。

技術(shù)實(shí)現(xiàn)要素：

有鑒于此，本申請(qǐng)針對(duì)TS-DInSAR難以滿足快速或大梯度地表形變監(jiān)測(cè)需求的問(wèn)題，提供了一種適用于大梯度地表沉降監(jiān)測(cè)的時(shí)序差分雷達(dá)干涉方法。

為了解決上述技術(shù)問(wèn)題，本申請(qǐng)公開(kāi)了一種適用于大梯度地表沉降監(jiān)測(cè)的時(shí)序差分雷達(dá)干涉方法，具體按照以下步驟實(shí)施：

步驟1，對(duì)所有篩選后的SAR影像進(jìn)行任意干涉組合并計(jì)算時(shí)間基線(形成一個(gè)干涉對(duì)的兩幅SAR影像的獲取時(shí)間差)和空間基線(兩次成像時(shí)刻SAR傳感器之間的空間距離，一般取該距離垂直于SAR視線方向的分量)；

步驟2，設(shè)置時(shí)空基線閾值進(jìn)行干涉(干涉組合)對(duì)初始篩選，在保證干涉對(duì)數(shù)量的前提下限制時(shí)間失相干和空間失相干，按照干涉組合進(jìn)行差分干涉處理得到差分干涉圖集；

步驟3，永久散射體(Persistent Scatterer,PS)和相干目標(biāo)(Coherent Target,CT)探測(cè)及點(diǎn)集合并，得到相干散射體(Coherent Scatterer,CS)點(diǎn)集，提取CS點(diǎn)集上的差分干涉相位，并以CS為節(jié)點(diǎn)構(gòu)建不規(guī)則三角網(wǎng)(Triangular Irregular Network,TIN)，本發(fā)明技術(shù)方案中采用Delaunay三角網(wǎng)；

步驟4，基于Delaunay三角網(wǎng)和最小費(fèi)用流(Minimum Cost Flow,MCF)方法的CS相位解纏；

步驟5，線性形變速率和高程誤差建模及解算；

步驟6，基于數(shù)據(jù)模擬或地面測(cè)量數(shù)據(jù)的線性形變解算結(jié)果檢驗(yàn)，若結(jié)果不可靠則執(zhí)行步驟7，否則轉(zhuǎn)向步驟8；

步驟7，重新篩選參與計(jì)算的差分干涉圖，并重新執(zhí)行步驟5和步驟6，當(dāng)達(dá)到計(jì)算終止條件時(shí)，執(zhí)行步驟8；

步驟8，從所有原始差分干涉圖中減去線性形變相位分量，對(duì)殘差相位(從差分干涉圖中減去線性形變相位后剩余的相位)進(jìn)行重新解纏，然后再將線性形變相位分量加回重新解纏后的相位中，重新執(zhí)行步驟5，得到更新后的線性形變速率和高程誤差(這一步主要是讓更多干涉對(duì)參與計(jì)算，提高線性形變速率和高程誤差的解算精度)；

步驟9，記錄步驟8中所得到的線性形變速率，并從所有的原始差分干涉圖中減去更新后的線性形變和高程誤差相位分量，重新執(zhí)行基于離散點(diǎn)的相位解纏，然后將步驟8中的線性形變分量重新加回新的解纏相位中，執(zhí)行形變時(shí)間序列建模和解算過(guò)程，最終得到形變時(shí)間序列。

與現(xiàn)有技術(shù)相比，本申請(qǐng)可以獲得包括以下技術(shù)效果：

本發(fā)明技術(shù)方案依據(jù)TS-DInSAR時(shí)序差分干涉圖中形變相位大小以及形變相位梯度大小與時(shí)間間隔(差分干涉圖的時(shí)間基線，即形成差分干涉圖的兩幅SAR影像獲取時(shí)間差)成正相關(guān)關(guān)系這一特點(diǎn)，提出短時(shí)間基線差分干涉圖篩選、離散點(diǎn)相位解纏、基于短時(shí)間基線差分干涉圖的形變分量建模和解算、形變分量可靠性檢驗(yàn)、差分干涉圖相位梯度修正、對(duì)上述過(guò)程進(jìn)行迭代以確保形變分量解算正確，以及基于修正后差分干涉圖的形變時(shí)間序列建模和解算這一整體的技術(shù)方案，解決原有TS-DInSAR技術(shù)中由于形變和相位梯度較大導(dǎo)致形變相位模糊度解算困難或失敗的問(wèn)題，最終達(dá)到正確提取大梯度地表形變速率和形變時(shí)間序列的目的。

采用本發(fā)明技術(shù)方案，只要差分干涉圖序列中存在至少6個(gè)短時(shí)間基線干涉對(duì)(由4幅SAR影像構(gòu)成)即可實(shí)現(xiàn)線性形變分量的解算，然后開(kāi)展長(zhǎng)時(shí)間基線差分干涉圖的形變相位梯度修正，促使更多的差分干涉圖得到正確解纏，并參與形變分量的重新估算以及形變時(shí)間序列的建模和解算。因此，無(wú)需很高的SAR影像使用量便可實(shí)現(xiàn)大梯度形變的有效提取。

當(dāng)然，實(shí)施本申請(qǐng)的任一產(chǎn)品并不一定需要同時(shí)達(dá)到以上所述的所有技術(shù)效果。

附圖說(shuō)明

此處所說(shuō)明的附圖用來(lái)提供對(duì)本申請(qǐng)的進(jìn)一步理解，構(gòu)成本申請(qǐng)的一部分，本申請(qǐng)的示意性實(shí)施例及其說(shuō)明用于解釋本申請(qǐng)，并不構(gòu)成對(duì)本申請(qǐng)的不當(dāng)限定。在附圖中：

圖1是一種適用于大梯度地表沉降監(jiān)測(cè)的時(shí)序差分雷達(dá)干涉方法的實(shí)施流程圖；

圖2是35幅SAR影像任意組合所形成的干涉對(duì)的時(shí)間基線和空間基線分布圖；

圖3是空間基線閾值為30米時(shí)干涉對(duì)的時(shí)空基線分布圖；

圖4是第一次解算所得CS點(diǎn)上線性形變速率結(jié)果圖；

圖5是使用第一次解算所得的線性形變模擬的差分干涉圖與原始差分干涉圖的對(duì)比；其中，a為原始差分干涉圖，b為使用第一次解算所得的線性形變模擬的差分干涉圖；

圖6是經(jīng)3次迭代后解算所得CS點(diǎn)上的線性形變速率結(jié)果圖；

圖7是使用第三次解算所得的線性形變模擬的差分干涉圖與原始差分干涉圖的對(duì)比；其中，a為原始差分干涉圖，b為使用第三次解算所得的線性形變模擬的差分干涉圖；

圖8是2009年6月23日至2010年7月2日期間的累積形變量；

圖9是圖8中所標(biāo)示三個(gè)特征點(diǎn)(A、B和C)的形變時(shí)間序列。

具體實(shí)施方式

以下將配合附圖及實(shí)施例來(lái)詳細(xì)說(shuō)明本申請(qǐng)的實(shí)施方式，藉此對(duì)本申請(qǐng)如何應(yīng)用技術(shù)手段來(lái)解決技術(shù)問(wèn)題并達(dá)成技術(shù)功效的實(shí)現(xiàn)過(guò)程能充分理解并據(jù)以實(shí)施。

一種適用于大梯度地表沉降監(jiān)測(cè)的時(shí)序差分雷達(dá)干涉方法，如附圖1所示，具體按照以下步驟實(shí)施：

步驟1，對(duì)所有篩選后的SAR影像進(jìn)行任意干涉組合并計(jì)算時(shí)間基線和空間基線；

在篩選出合適的SAR影像(排除受雨雪等天氣以及積雪影響的SAR影像)后，進(jìn)行任意干涉組合配對(duì)，假設(shè)有N+1幅SAR影像，通過(guò)任意干涉組合可形成N(N+1)/2個(gè)干涉對(duì)。每個(gè)干涉對(duì)由主、從兩幅SAR影像構(gòu)成。然后，根據(jù)每個(gè)干涉對(duì)主、從SAR影像的獲取時(shí)間計(jì)算該干涉對(duì)的時(shí)間基線(即SAR影像獲取的時(shí)間差)，根據(jù)主、從SAR影像的參數(shù)文件所記錄的SAR傳感器運(yùn)行位置及相應(yīng)的時(shí)間參數(shù)計(jì)算該干涉對(duì)的空間基線(主、從影像成像時(shí)SAR傳感器的空間距離，實(shí)際中一般取該空間距離在垂直于SAR傳感器視線方向上的分量為空間基線，也即垂直基線)。

本實(shí)施例中采用覆蓋天津市精武鎮(zhèn)的35幅SAR影像為數(shù)據(jù)進(jìn)行展示。附圖2所示為35幅SAR影像進(jìn)行任意組合所形成的干涉對(duì)的時(shí)間基線和空間基線分布。圖2中以“年年年年月月日日”的方式標(biāo)注了每幅SAR影像的獲取時(shí)間，如“20090327”代表2009年3月27日。通過(guò)虛線相連的兩幅SAR影像為一個(gè)干涉對(duì)，縱軸代表干涉對(duì)的空間基線大小，橫軸代表時(shí)間，時(shí)間的差值即為時(shí)間基線。

步驟2，設(shè)置時(shí)空基線閾值進(jìn)行干涉對(duì)的初始篩選，在保證干涉對(duì)數(shù)量的前提下限制時(shí)間失相干和空間失相干，按照干涉組合進(jìn)行差分干涉處理得到差分干涉圖集；

在進(jìn)行時(shí)序差分干涉處理時(shí)，為降低時(shí)間失相干和空間失相干的影響，常采用時(shí)間基線閾值和空間基線閾值方法排除時(shí)間基線和空間基線大于某給定閾值的干涉對(duì)(受時(shí)間失相干和空間失相干影響較顯著的干涉對(duì))。大量研究表明，時(shí)間基線不宜超過(guò)2年，空間基線不宜超過(guò)150米。但是，當(dāng)SAR影像獲取時(shí)間整體跨度較小時(shí)(如第一幅SAR影像和最后一幅SAR影像的獲取時(shí)間間隔小于2年)可不考慮對(duì)時(shí)間基線進(jìn)行限制。在本實(shí)施例中，考慮到SAR影像獲取時(shí)間的整體跨度較短(1年零9個(gè)月)，故未設(shè)置時(shí)間基線閾值，僅考慮對(duì)空間基線進(jìn)行限制，設(shè)置空間基線閾值為30米，排除空間基線大于30米的干涉對(duì)。附圖3所示為排除空間基線大于30米的干涉對(duì)后，得到保留的干涉對(duì)的時(shí)間基線和空間基線分布?；诖烁缮娼M合，進(jìn)行差分干涉處理，得到相應(yīng)的差分干涉圖集，處理過(guò)程如下：

以其中一個(gè)干涉對(duì)為例，設(shè)組成此干涉對(duì)的兩幅SAR影像分別為I₁和I₂，分別記為主影像和從影像，對(duì)于影像中所記錄的一點(diǎn)P(即影像像元P)，其在I₁和I₂中所對(duì)應(yīng)的同名像元的值分別為：

Ψ₁＝A₁exp(jφ₁) (1)

Ψ₂＝A₂exp(jφ₂) (2)

其中，exp(·)代表以e為底的指數(shù)函數(shù)；Ψ₁為目標(biāo)P在I₁中的像元值，φ₁為目標(biāo)P在I₁中的相位值；Ψ₂為目標(biāo)P在I₂中的像元值，φ₂為目標(biāo)P在I₂中的相位值；A₁為目標(biāo)P在I₁中信號(hào)的振幅值；A₂為目標(biāo)P在I₂中信號(hào)的振幅值。其中，Ψ₁和Ψ₂為復(fù)數(shù)。由復(fù)數(shù)的共軛相乘可得P在由I₁和I₂形成的干涉圖中的像元值為：

Ψ_int，P＝Ψ₁Ψ₂^*＝A₁A₂exp(j(φ₁-φ₂)) (3)

其中，“*”代表復(fù)數(shù)的共軛，即Ψ₂^*為Ψ₂的共軛復(fù)數(shù)。

由式(3)可得P點(diǎn)所對(duì)應(yīng)的干涉相位(即兩次觀測(cè)的相位差表達(dá)形式)為：

φ_int,P＝φ₁-φ₂ (4)

對(duì)SAR影像中的所有像元進(jìn)行相同的處理，即可得到由主、從SAR影像I₁和I₂所形成的干涉圖。以P點(diǎn)為例，此時(shí)，P點(diǎn)的干涉相位可表達(dá)為：

φ_int,P＝φ_ref,P+φ_top,P+φ_def,P+φ_atm,P+φ_add,P (5)

其中，φ_int,P為P點(diǎn)在由SAR影像I₁和I₂所形成的干涉圖中的干涉相位；φ_ref,P為P點(diǎn)的參考橢球面相位，由SAR傳感器兩次對(duì)同一地區(qū)觀測(cè)時(shí)所處位置不同以及地球曲面引起；φ_top,P為P點(diǎn)的地形相位，由地表真實(shí)高程及其變化引起；φ_def,P為P點(diǎn)所發(fā)生的地表形變對(duì)應(yīng)的相位；φ_atm,P為P點(diǎn)上的大氣延遲相位，由雷達(dá)波在大氣層中傳播發(fā)生折射引起；φ_add,P是P點(diǎn)上的附加相位，主要包括各種隨機(jī)噪聲以及失相干噪聲(兩次成像時(shí)地面目標(biāo)發(fā)生變化，導(dǎo)致雷達(dá)回波信號(hào)產(chǎn)生差異所致)。

所謂差分干涉處理，即在上述干涉處理基礎(chǔ)上，去除地球參考橢球面相位和地形相位，得到以形變相位為主的差分干涉相位，形成差分干涉圖。首先，在不考慮大氣延遲和噪聲相位的情況下，P點(diǎn)上的干涉相位可表達(dá)為：

其中，等式右邊第一項(xiàng)即為參考橢球面相位的表達(dá)式，第二項(xiàng)為地形相位的表達(dá)式，第三項(xiàng)為形變相位的表達(dá)式；B為干涉對(duì)的空間基線；θ為雷達(dá)側(cè)視角；α為空間基線與水平方向的夾角；λ為雷達(dá)波長(zhǎng)；R₁為雷達(dá)傳感器到觀測(cè)目標(biāo)之間的斜距；為空間基線B在垂直于雷達(dá)側(cè)視方向上的投影分量，即垂直基線；h_P為P點(diǎn)的高程；Δr為P點(diǎn)所發(fā)生的地表形變。需要說(shuō)明的是，Bsin(θ-α)即為干涉對(duì)的平行基線(空間基線B在沿雷達(dá)側(cè)視方向上的投影分量)且有：

其中，R₁和R₂分別為雷達(dá)傳感器兩次觀測(cè)成像時(shí)到地面點(diǎn)P的斜距；ΔR＝R₁-R₂為兩次成像雷達(dá)傳感器到P點(diǎn)斜距的差。

差分干涉處理即要從干涉相位中減去公式(6)等號(hào)右邊的前兩項(xiàng)。首先需要計(jì)算并扣除參考橢球面相位。由公式(6)等號(hào)右邊第一項(xiàng)以及公式(7)可知，參考橢球面相位與衛(wèi)星兩次觀測(cè)至參考橢球面間斜距的差ΔR有關(guān)，即：

由公式(7)可知，要獲得ΔR的值，需首先已知R₁和R₂，此二者可通過(guò)兩點(diǎn)間距離公式計(jì)算：

其中，d(·)為兩點(diǎn)間距離方程；M和S分別為對(duì)應(yīng)于主、從影像中P點(diǎn)成像時(shí)刻的衛(wèi)星三維空間坐標(biāo)矢量，由SAR影像參數(shù)文件提供的參數(shù)計(jì)算得到；P為P點(diǎn)在參考橢球面上的空間三維坐標(biāo)矢量，可通過(guò)衛(wèi)星成像參數(shù)、斜距方程、多普勒方程和橢球方程計(jì)算得到。此時(shí)，可通過(guò)式(8)和(9)得到P點(diǎn)上的參考橢球面相位。

對(duì)于地形相位，通過(guò)外部數(shù)字高程模型(Digital Elevation Model,DEM)可獲取地面點(diǎn)P的高程值h_P(美國(guó)地質(zhì)調(diào)查局網(wǎng)站公布了全球免費(fèi)的DEM，本實(shí)施例中即采用該DEM)，R₁為主影像斜距，可由P點(diǎn)在主、從影像上的成像參數(shù)計(jì)算得到。對(duì)于垂直基線有：

其中，空間基線B＝d(M,S)，即主、從影像成像時(shí)雷達(dá)傳感器的空間距離；在得到h_P、R₁、和后，可通過(guò)式(6)等號(hào)右邊第二項(xiàng)計(jì)算P點(diǎn)的地形相位。

在得到P點(diǎn)的參考橢球面相位和地形相位后，從式(5)中扣除這兩部分相位分量即可得到差分干涉相位。對(duì)SAR影像中所有的像元點(diǎn)進(jìn)行上述差分干涉處理，即可得到相應(yīng)的差分干涉圖。

步驟3，PS和CT探測(cè)及點(diǎn)集合并，得到CS點(diǎn)集，提取CS點(diǎn)集上的差分干涉相位，并以CS為節(jié)點(diǎn)構(gòu)建不規(guī)則三角網(wǎng)(Triangular Irregular Network,TIN)，本發(fā)明技術(shù)方案中采用Delaunay三角網(wǎng)；

PS探測(cè)采用振幅閾值和振幅離差指數(shù)(Amplitude Dispersion Index,ADI)閾值雙閾值方法。對(duì)于SAR影像中某一特定的像元，其在N+1幅SAR影像中的振幅值可直接從SAR振幅影像中提取，則其時(shí)序振幅均值和ADI值為：

其中，i為SAR影像索引號(hào)(序號(hào))；a_i和分別為該像元在第i幅SAR影像中的振幅值及在所有影像中的振幅均值；D_a和σ_a分別為該像元的ADI值及時(shí)序振幅標(biāo)準(zhǔn)差。當(dāng)該像元的振幅信息滿足如下兩個(gè)條件時(shí)，認(rèn)為其為PS：

其中，和σ_A分別為N+1幅SAR影像所有像元時(shí)序振幅值的均值和標(biāo)準(zhǔn)差；c和l分別為像元的列號(hào)和行號(hào)，C和L分別為影像列數(shù)和行數(shù)；A_i和a_cl分別為第i幅SAR影像所有像元振幅值的均值和行、列號(hào)分別為c和l的某像元的振幅值。公式(13)中0.25即為ADI閾值。

CT探測(cè)采用相干系數(shù)閾值方法。假設(shè)由N+1幅SAR影像形成了L個(gè)干涉對(duì)，并通過(guò)差分干涉數(shù)據(jù)處理得到L幅差分干涉圖。通過(guò)下式計(jì)算每個(gè)像元在所有干涉圖中相應(yīng)的相干系數(shù)：

其中，為某像元在第l幅干涉圖中的相干系數(shù)值；IM_l和IS_l分別為第l個(gè)干涉對(duì)的主影像和從影像；Z為相干系數(shù)估計(jì)窗口內(nèi)像元數(shù)目，z為窗口內(nèi)像元索引；l為干涉圖索引。當(dāng)某個(gè)像元滿足以下條件時(shí)被識(shí)別為CT：

其中，min(·)表示取變量的最小值；γ_crit為判別某個(gè)像元是否為CT的相干系數(shù)閾值，一般可取0.25至0.3，本實(shí)施例中采用0.25。

當(dāng)探測(cè)出所有的PS和CT后，將二者進(jìn)行合并，并去除重復(fù)點(diǎn)，得到CS點(diǎn)集，最后從所有差分干涉圖中提取CS點(diǎn)上的差分干涉相位值。

步驟4，基于Delaunay三角網(wǎng)和MCF方法的CS相位解纏；

相位解纏一般是基于規(guī)則格網(wǎng)進(jìn)行，常規(guī)的相位解纏方法針對(duì)差分干涉圖進(jìn)行處理，易受失相干區(qū)域的影響，而基于離散點(diǎn)的相位解纏方法可有效避免這一問(wèn)題。本發(fā)明技術(shù)方案中采用基于離散點(diǎn)的相位解纏方法對(duì)CS上的差分干涉相位進(jìn)行解纏處理。對(duì)于某個(gè)CS點(diǎn)P，其在差分干涉圖中的相位和該點(diǎn)上的真實(shí)相位的關(guān)系為：

其中，為P點(diǎn)上的真實(shí)相位(即相位解纏要得到的相位值)；φ_P為P點(diǎn)上的纏繞相位值(即P在差分干涉圖中對(duì)應(yīng)的相位值)，位于[-π,π)區(qū)間內(nèi)，只記錄了不足整周(2π)的小數(shù)部分，存在整周模糊度；n_P為整周模糊數(shù)。

首先，根據(jù)Delaunay剖分法則，以所有的CS為頂點(diǎn)構(gòu)造Delaunay三角網(wǎng)，然后以網(wǎng)絡(luò)邊端點(diǎn)對(duì)應(yīng)的兩個(gè)CS之間的纏繞相位差為觀測(cè)量，估算兩點(diǎn)間的絕對(duì)相位差。通過(guò)最小費(fèi)用流(Minimum Cost Flow,MCF)方法對(duì)與相位不連續(xù)性相關(guān)的網(wǎng)絡(luò)費(fèi)用流進(jìn)行估算，并尋找最小費(fèi)用流對(duì)應(yīng)的積分路徑，進(jìn)行相位積分，完成相位解纏。此過(guò)程也即求解公式(16)中n_P的過(guò)程。

步驟5，線性形變速率和高程誤差建模及解算；

在對(duì)CS上的差分干涉相位進(jìn)行相位解纏后，根據(jù)步驟4中構(gòu)建的Delaunay三角網(wǎng)重新計(jì)算相鄰CS點(diǎn)間的解纏相位差。此處，相鄰的CS點(diǎn)是指Delaunay三角網(wǎng)中三角形邊的端點(diǎn)。以任意一邊的兩端點(diǎn)P和Q為例，二者在第l個(gè)差分干涉圖中對(duì)應(yīng)于的解纏相位為：

其中，和分別為第l幅差分干涉圖的時(shí)間基線和垂直基線；λ為雷達(dá)波波長(zhǎng)；θ_P和θ_Q分別為P和Q點(diǎn)處的雷達(dá)波入射角；R_P和R_Q分別為SAR傳感器到P和Q之間的斜距；和分別為P和Q點(diǎn)上解纏后的差分干涉相位；v_P和v_Q分別為P和Q點(diǎn)的線性形變速率；δh_P和δh_Q分別為P和Q點(diǎn)上的高程誤差；和分別為P和Q點(diǎn)在第l幅差分干涉圖中的非線性形變相位；和分別為P和Q點(diǎn)在第l幅差分干涉圖中的軌道誤差相位；和分別為P和Q點(diǎn)在第l幅差分干涉圖中的大氣延遲相位；和分別為P和Q點(diǎn)在第l幅差分干涉圖中的噪聲相位。

線性形變和高程誤差建模及解算的目的是對(duì)v_P和v_Q及δh_P和δh_Q進(jìn)行估算。在本發(fā)明技術(shù)方案中，以P和Q上的解纏相位為觀測(cè)量進(jìn)行網(wǎng)絡(luò)鄰域差分建模。對(duì)于P和Q，二者之間的網(wǎng)絡(luò)鄰域差分相位增量(解纏相位差)為：

其中，和分別為P和Q點(diǎn)處斜距的平均值和入射角的平均值；為P和Q之間的鄰域差分相位增量；Δv_PQ為P和Q之間的線性形變速率增量(即線性形變速率差)；Δδh_PQ為P和Q之間的高程誤差增量(即高程誤差之差)；為第l幅差分干涉圖中P和Q之間的相位殘差增量(即相位殘差之差)，即P和Q點(diǎn)上非線性形變、軌道誤差、大氣延遲和噪聲相位之間的差值之和。

以L幅差分干涉圖為例，對(duì)于P和Q所對(duì)應(yīng)的網(wǎng)絡(luò)邊而言，可列出L個(gè)與式(19)相同的方程，組成相應(yīng)的線性方程組。將其表達(dá)為矩陣的形式有：

ΔΨ＝AX+W (20)

其中，

A＝[κ,η] (24)

其中，

此時(shí)，方程組中僅有Δv_PQ和Δδh_PQ兩個(gè)未知數(shù)，觀測(cè)值數(shù)量為L(zhǎng)，通過(guò)最小二乘可得未知參數(shù)的估值為：

對(duì)于所有通過(guò)網(wǎng)絡(luò)邊連接的CS點(diǎn)，均通過(guò)該過(guò)程解算相鄰兩點(diǎn)間的線性形變速率增量和高程誤差增量。在得到所有CS點(diǎn)間的線性形變速率增量和高程誤差增量后，以此作為觀測(cè)量，以每個(gè)CS點(diǎn)的線性形變速率和高程誤差為待估參數(shù)，在空間域進(jìn)行最小二乘網(wǎng)絡(luò)平差計(jì)算，解算所有CS點(diǎn)上的線性形變速率和高程誤差。對(duì)于線性形變速率有：

其中，和分別為點(diǎn)P和Q的線性形變速率估計(jì)值(待估參數(shù))；為點(diǎn)P和Q之間的線性形變速率增量估值；r_PQ為最小二乘殘差。對(duì)于所有的CS點(diǎn)(設(shè)數(shù)量為G)和網(wǎng)絡(luò)邊(設(shè)數(shù)量為H)而言，根據(jù)公式(28)得到相應(yīng)的觀測(cè)方程，表達(dá)為矩陣的形式為：

其中，B為系數(shù)矩陣，其元素由弧段所對(duì)應(yīng)的方程式確定，由1、-1和0組成；矩陣中為待估參數(shù)，即每個(gè)CS點(diǎn)的線性形變速率；中為所有CS點(diǎn)線性形變速率增量的估值；r為殘差矩陣。由最小二乘平差可得：

其中，P為線性速率增量的權(quán)陣，可通過(guò)弧段最小二乘殘差進(jìn)行估算得到。

如附圖4所示為使用附圖3中所展示的干涉對(duì)所解算得出的CS點(diǎn)上的線性形變速率。圖4中以分級(jí)設(shè)色的方式表達(dá)線性形變速率的量值，LOS DR表示雷達(dá)視線向形變速率，mm/yr為形變速率單位，即毫米/年。在步驟6中，將對(duì)此結(jié)果進(jìn)行檢驗(yàn)。

步驟6，基于數(shù)據(jù)模擬或地面測(cè)量數(shù)據(jù)的線性形變解算結(jié)果檢驗(yàn)，若結(jié)果不可靠則執(zhí)行步驟7，否則轉(zhuǎn)向步驟8；

一般地，對(duì)TS-DInSAR形變結(jié)果進(jìn)行驗(yàn)證主要有兩種途徑：基于外部測(cè)量數(shù)據(jù)的驗(yàn)證方法和基于差分干涉圖模擬的驗(yàn)證方法。由于外部數(shù)據(jù)往往不容易獲得，因此本發(fā)明技術(shù)方案中采用基于差分干涉圖模擬的驗(yàn)證方法。

該方法是使用估計(jì)所得到的年形變速率或累積形變量和相應(yīng)的干涉基線參數(shù)模擬差分干涉圖。即將形變量轉(zhuǎn)換為形變相位，形變到相位的轉(zhuǎn)換公式為：

其中，為線性形變?cè)诘趌幅差分干涉圖中對(duì)應(yīng)的相位；v為通過(guò)步驟5中的數(shù)據(jù)處理所得的線性形變速率。

在大形變梯度區(qū)域(如沉降漏斗區(qū))，差分干涉圖中的形變相位梯度隨時(shí)間基線的增大而顯著提升，在差分干涉圖中表現(xiàn)為密集的條紋變化，若估計(jì)所得形變量符合或接近真實(shí)情況，使用形變結(jié)果模擬的差分干涉圖與原始差分干涉圖應(yīng)該相似，尤其是形變條紋數(shù)目應(yīng)相等或差別較小。若估計(jì)所得形變量不符合真實(shí)情況(如欠估計(jì))，則模擬所得條紋數(shù)目與原始差分干涉圖會(huì)存在明顯區(qū)別。據(jù)此，可簡(jiǎn)單有效地對(duì)形變解算結(jié)果的有效性進(jìn)行驗(yàn)證，且這種驗(yàn)證方法能夠有效識(shí)別形變欠估計(jì)現(xiàn)象。

如附圖5所示為由步驟5解算所得線性形變模擬所得的差分干涉圖與原始差分干涉圖的對(duì)比。其中，左側(cè)(圖5a)為原始差分干涉圖，右側(cè)(圖5b)為模擬的差分干涉圖。二者條紋數(shù)目差別明顯，模擬所得差分干涉圖中條紋數(shù)目低于原始差分干涉圖，表明發(fā)生了形變欠估計(jì)現(xiàn)象。因此，第一次解算的線性形變并不可靠。此時(shí)，需要執(zhí)行步驟7。

步驟7，重新篩選參與計(jì)算的差分干涉圖，并重新執(zhí)行步驟5和步驟6，當(dāng)達(dá)到計(jì)算終止條件時(shí)，執(zhí)行步驟8；

該步驟是實(shí)現(xiàn)正確解算線性形變的關(guān)鍵。經(jīng)過(guò)步驟6的檢驗(yàn)，確定所解算結(jié)果不可靠時(shí)，則通過(guò)設(shè)置更為嚴(yán)格的時(shí)間基線重新篩選出短時(shí)間基線差分干涉圖，重新執(zhí)行步驟5和步驟6，并進(jìn)行檢驗(yàn)。具體地，重新篩選短時(shí)間基線差分干涉圖是指設(shè)置更小的時(shí)間基線，篩選出比上一次解算所用干涉圖時(shí)間基線更短的干涉圖。這是由于在一般的地表沉降區(qū)，時(shí)間基線越短則形變?cè)叫?，且形變梯度越小，差分干涉圖中的形變條紋數(shù)目就越少，更容易對(duì)差分干涉圖進(jìn)行正確的相位解纏，為形變的建模和解算提供正確的觀測(cè)量。例如，本實(shí)施例中第一次解算時(shí)未對(duì)時(shí)間基線進(jìn)行限制，解算所得到的線性形變不可靠，則設(shè)置較短的時(shí)間基線值(本實(shí)施例中采用180天)，篩選出時(shí)間基線小于180天的干涉圖重新解算線性形變，然后再次進(jìn)行結(jié)果檢驗(yàn)，若不可靠，則進(jìn)行第三次解算，本實(shí)施例中第三次解算時(shí)選取90天作為時(shí)間基線閾值，即篩選出時(shí)間基線小于90天的干涉圖重新解算線性形變。時(shí)間基線逐漸縮短，即所謂的設(shè)置更為嚴(yán)格的時(shí)間基線限制。

本實(shí)施例中，經(jīng)過(guò)3次迭代計(jì)算可獲取正確的線性形變結(jié)果，如附圖6所示。附圖7所示為利用第三次解算所得線性形變模擬的差分干涉圖與原始差分干涉圖的對(duì)比，其中左側(cè)(圖7a)為原始差分干涉圖，右側(cè)(圖7b)為模擬的差分干涉圖。此時(shí)，二者條紋數(shù)目幾乎一致。得到可靠的線性形變后，可執(zhí)行步驟8。

需要指出的是，在本實(shí)施例中，通過(guò)3次迭代計(jì)算獲取了正確的線性形變速率，但對(duì)于不同的研究區(qū)域，迭代次數(shù)是可變的，如增加或減少，以能夠獲取正確的線性形變?yōu)闇?zhǔn)。

步驟8，從所有原始差分干涉圖中減去線性形變相位分量，對(duì)殘差相位(從差分干涉圖中減去線性形變相位后剩余的相位)進(jìn)行重新解纏，然后再將線性形變相位分量加回重新解纏后的相位中，重新執(zhí)行步驟5，得到最終的線性形變速率和高程誤差(這一步主要是讓更多干涉對(duì)參與計(jì)算，提高線性形變速率和高程誤差的解算精度)；

在通過(guò)第7步驟得到正確的線性形變后，從所有的原始差分干涉圖中減去線性形變所對(duì)應(yīng)的相位分量。為便于闡述，設(shè)此時(shí)所得到的正確的線性形變?yōu)槿砸缘趌幅差分干涉圖為例，其所對(duì)應(yīng)的相位分量為：

此時(shí)，通過(guò)下式獲取扣除線性形變分量后的相位：

其中，為扣除線性形變分量后的纏繞相位；φ_diff,l為原始的差分干涉相位；conj(·)為求復(fù)數(shù)共軛。對(duì)所有L幅差分干涉圖執(zhí)行此步運(yùn)算，得L幅修正時(shí)空相位梯度后的差分干涉圖。此時(shí)，差分干涉圖中的相位梯度將顯著降低，條紋數(shù)目減少，有利于相位解纏的正確執(zhí)行。

在得到L幅經(jīng)過(guò)時(shí)空相位梯度修正的差分干涉圖后，重新對(duì)干涉圖集執(zhí)行步驟4中所述的相位解纏過(guò)程，得到L幅相應(yīng)的解纏相位圖，新的解纏相位表達(dá)為：

其中，為扣除線性形變分量后的解纏相位；為高程誤差相位分量；為非線性形變相位分量；為軌道誤差相位分量；為大氣延遲相位分量；為噪聲相位分量。在得到修正后差分干涉相位的解纏值后，將步驟7中所得線性形變分量(即所對(duì)應(yīng)的形變相位)與其相加，得到修正后的解纏相位：

在得到修正的解纏相位后，以其為觀測(cè)值重新執(zhí)行步驟5中的解算流程，得到最終的線性形變速率和高程誤差。為便于后續(xù)闡述，此處設(shè)最終的線性形變速率和高程誤差分別為和則二者所對(duì)應(yīng)的相位為：

現(xiàn)有研究表明，參與運(yùn)算的干涉對(duì)數(shù)目越多，線性形變速率和高程誤差的解算精度越高，這也是對(duì)相位梯度進(jìn)行修正使更多干涉對(duì)參與運(yùn)算的重要目的之一。此外，在后續(xù)的數(shù)據(jù)處理中，最終的線性形變和高程誤差將用于最終的相位梯度修正。

步驟9，記錄步驟8中所得到的線性形變速率，并從所有的原始差分干涉圖中減去更新后的線性形變和高程誤差相位分量，重新執(zhí)行基于離散點(diǎn)的相位解纏，然后將步驟8中的線性形變分量重新加回新的解纏相位中，執(zhí)行形變時(shí)間序列建模和解算過(guò)程，最終得到形變時(shí)間序列。

在得到最終的線性形變速率和高程誤差結(jié)果后，首先保留并輸出線性形變速率結(jié)果，然后可根據(jù)下式將線性形變和高程誤差相位(可通過(guò)公式(36)和(37)計(jì)算)從原始差分干涉相位中扣除，得到新的差分干涉相位值：

其中，為扣除線性形變和高程誤差相位后的差分干涉相位，且有：

在得到最終經(jīng)線性形變和高程誤差修正的相位后，重新執(zhí)行第4步驟中的相位解纏過(guò)程。在得到新的解纏相位后，重新將步驟8中得到的線性形變相位加回，得到扣除高程誤差后的解纏相位：

其中，為經(jīng)高程誤差校正后的解纏相位；為由公式(36)計(jì)算得到的線性形變相位。

對(duì)于某個(gè)CS而言，其在第l幅扣除了高程誤差的解纏相位圖中的相位表達(dá)為：

其中，為形變相位(包括線性形變和非線性形變)；為軌道誤差相位；為大氣延遲相位；為噪聲相位。對(duì)于L個(gè)差分干涉圖，經(jīng)高程誤差改正后的相位為：

仍以上述N+1幅SAR影像為例，設(shè)其獲取時(shí)刻為T＝[t₀,t₁,…,t_n,…,t_N]，由其所形成的L個(gè)干涉對(duì)的主(IM)、從(IS)SAR影像集分別為：

IM＝[IM₁,IM₂,…,IM_l,…,IM_L]；IS＝[IS₁,IS₂,…,IS_l,…,IS_L] (43)

此時(shí)，以t₀時(shí)刻為參考時(shí)刻，其余任意時(shí)刻t_n(n＝1,2,…,N)相對(duì)于t₀時(shí)刻的相位為待估參數(shù)，以扣除高程誤差后的解纏相位為觀測(cè)值，建立觀測(cè)方程并恢復(fù)每個(gè)時(shí)刻對(duì)應(yīng)的相位(即相位時(shí)間序列)。為得到符合實(shí)際物理意義的形變估計(jì)參數(shù)，假設(shè)相鄰兩幅SAR影像獲取時(shí)間間隔內(nèi)相位的變化符合線性累積趨勢(shì)，將對(duì)相位時(shí)間序列的求解轉(zhuǎn)化為對(duì)相位變化速率(相鄰時(shí)間段內(nèi)的平均相位變化速率)v_ph(注意與步驟5、6、7和8中形變速率v的區(qū)別)的求解，此時(shí)，待估參數(shù)為：

其中，為t_n時(shí)刻相對(duì)于t₀時(shí)刻的累積相位，且有式(44)實(shí)際的物理意義是相鄰兩幅SAR影像獲取時(shí)間間隔內(nèi)的相位變化速率，也稱為分段相位速率。相應(yīng)地，該模型可被稱為分段線性模型。此處的“段”指相鄰兩幅SAR影像之間的時(shí)間段。此時(shí)可得觀測(cè)方程為：

將式(45)表達(dá)為矩陣的形式有：

其中，B為L(zhǎng)×N的矩陣，矩陣元素B[l,n]＝t_n-t_n-1(IS_l+1≤n≤IM_l)，其它元素值為零。由于干涉圖集可能是不連續(xù)的，即在某個(gè)SAR影像獲取時(shí)刻斷開(kāi)，進(jìn)而造成系數(shù)矩陣B發(fā)生奇異(此時(shí)方程組為欠定狀態(tài)，不存在唯一解)。因此，實(shí)際處理中首先對(duì)B進(jìn)行奇異值分解(Singular Value Decomposition,SVD)處理，然后估算出最小范數(shù)意義下各時(shí)間段的平均相位速率，通過(guò)在時(shí)間域上進(jìn)行積分即可恢復(fù)與SAR影像獲取時(shí)刻相對(duì)應(yīng)的相位時(shí)間序列v_ph。B的奇異值分解為：

B＝USW^T (47)

其中，U為L(zhǎng)×L階正交矩陣，被稱作B的左奇異向量，其前N列是BB^T的特征向量；W是N×L階正交矩陣，被稱作B的右奇異向量。S是半正定L×L階對(duì)角矩陣，其元素為相應(yīng)于BB^T特征向量的均方根，也即奇異值σ_i，且有，

S＝diag(σ₁,…,σ_N-r+1,0,…,0) (48)

其中，diag(·)代表對(duì)角矩陣，除對(duì)角線元素為σ_i外，其余元素值均為零；N-r+1為矩陣B的秩，r代表差分干涉圖子集的數(shù)量(即由于設(shè)置了時(shí)間基線和空間基線閾值限制，導(dǎo)致由N+1幅SAR影像所形成的干涉圖集合被分成了多個(gè)子集，造成干涉圖集不連續(xù))。最終，可對(duì)公式(46)進(jìn)行解算得：

其中，S⁺＝diag(1/σ₁,…,1/σ_N-r+1,…,0,…,0)。

在解算出相鄰時(shí)間段內(nèi)的平均相位變化速率后，根據(jù)相位變化速率與時(shí)間之積并求和得到對(duì)應(yīng)于SAR影像獲取時(shí)刻的相位序列，即的值。然后，從相位序列中扣除由步驟8所得到的線性形變分量，并在空間域進(jìn)行低通濾波去除噪聲相位，在時(shí)間域進(jìn)行高通濾波得到大氣延遲和軌道誤差相位序列。最后，從相位序列中扣除大氣延遲和軌道誤差相位序列即可得到地表形變所對(duì)應(yīng)的相位序列此時(shí)可得形變時(shí)間序列為：

其中，D為與每幅SAR影像獲取時(shí)刻相對(duì)應(yīng)的累積形變量所組成的向量；為與每幅SAR影像獲取時(shí)刻相對(duì)應(yīng)的累積形變相位。即：

如附圖8所示為2009年6月23日至2010年7月2日期間的累積形變量。圖8中標(biāo)記了A、B和C三個(gè)特征點(diǎn)，附圖9所示為三個(gè)特征點(diǎn)的形變時(shí)間序列。至此，便完成了大梯度形變區(qū)域(沉降漏斗區(qū)域)形變速率和形變時(shí)間序列的解算。

本發(fā)明技術(shù)方案帶來(lái)的有益效果：

(1)解決原有TS-DInSAR技術(shù)難以進(jìn)行大梯度地表形變建模和解算的問(wèn)題

原有TS-DInSAR技術(shù)難以適用于大梯度地表形變的建模和解算，本發(fā)明技術(shù)方案依據(jù)TS-DInSAR時(shí)序差分干涉圖中形變相位大小以及形變相位梯度大小與時(shí)間間隔(差分干涉圖的時(shí)間基線，即形成差分干涉圖的兩幅SAR影像獲取時(shí)間差)成正相關(guān)關(guān)系這一特點(diǎn)，提出短時(shí)間基線差分干涉圖篩選、離散點(diǎn)相位解纏、基于短時(shí)間基線差分干涉圖的形變分量建模和解算、形變分量可靠性檢驗(yàn)、差分干涉圖相位梯度修正、對(duì)上述過(guò)程進(jìn)行迭代以確保形變分量解算正確，以及基于修正后差分干涉圖的形變時(shí)間序列建模和解算這一整體的技術(shù)方案，解決原有TS-DInSAR技術(shù)中由于形變和相位梯度較大導(dǎo)致形變相位模糊度解算困難或失敗的問(wèn)題，最終達(dá)到正確提取大梯度地表形變速率和形變時(shí)間序列的目的。

針對(duì)本發(fā)明技術(shù)方案中具體的思路進(jìn)行有益效果分析如下：

(a)短時(shí)間基線差分干涉圖篩選。短時(shí)間基線差分干涉圖中形變相位梯度較小，較容易實(shí)現(xiàn)正確的相位模糊度解算(即相位解纏)，得到正確的解纏相位后，可用于后續(xù)形變分量建模和解算。此外，短時(shí)間基線差分干涉圖篩選是一個(gè)可以重復(fù)執(zhí)行的過(guò)程。

(b)離散點(diǎn)相位解纏。相位解纏也即解算相位模糊度，離散點(diǎn)相位解纏僅針對(duì)高質(zhì)量的觀測(cè)目標(biāo)進(jìn)行處理，避免低質(zhì)量目標(biāo)的影響。解纏后的相位將被用于形變分量建模。相比于原有技術(shù)基于非解纏相位進(jìn)行形變分量解算而言，基于解纏相位進(jìn)行形變分量建?？杀苊鈱?duì)SAR影像時(shí)間采樣率的依賴。

(c)基于短時(shí)間基線差分干涉圖的形變分量建模和解算，并檢驗(yàn)形變分量的可靠性。在(a)和(b)基礎(chǔ)上進(jìn)行形變分量(線性形變速率)建模和解算(方法見(jiàn)技術(shù)方案的詳細(xì)闡述部分)，并采用差分干涉圖模擬方法檢驗(yàn)形變分量解算結(jié)果的可靠性。先保證形變分量解算正確，后續(xù)將被應(yīng)用于對(duì)差分干涉圖的形變相位梯度進(jìn)行修正，即從差分干涉圖中減去形變分量，起到降低形變相位梯度的作用，進(jìn)而實(shí)現(xiàn)所有差分干涉圖的正確解纏。

(d)迭代計(jì)算。在對(duì)形變分量進(jìn)行可靠性檢驗(yàn)時(shí)，若發(fā)現(xiàn)結(jié)果不可靠，則重新進(jìn)行(a)、(b)和(c)，確保形變分量解算結(jié)果可靠。迭代計(jì)算是連接各個(gè)步驟的關(guān)鍵。

(e)采用正確的形變分量結(jié)果重新對(duì)差分干涉圖進(jìn)行相位梯度修正，然后重新進(jìn)行相位解纏，并基于此對(duì)形變時(shí)間序列進(jìn)行建模和解算。

綜上所述，以上各個(gè)過(guò)程緊密結(jié)合，共同形成了本發(fā)明的整體技術(shù)方案，最終實(shí)現(xiàn)大梯度形變區(qū)域形變分量和形變時(shí)間序列的精確解算。

(2)可降低SAR影像使用量，降低TS-DInSAR技術(shù)的應(yīng)用經(jīng)濟(jì)成本

在大梯度形變解算方面，現(xiàn)有研究指出SAR時(shí)間采樣率的提高可有效緩解梯度較大造成的問(wèn)題，但這往往受現(xiàn)有SAR系統(tǒng)觀測(cè)周期(對(duì)同一地區(qū)進(jìn)行重復(fù)觀測(cè)的周期)以及研究區(qū)域可用影像數(shù)量的限制，并且數(shù)據(jù)成本會(huì)顯著提高(如針對(duì)大梯度形變，以往的TS-DInSAR技術(shù)往往需要很高的SAR影像時(shí)間采樣率，即較高的數(shù)據(jù)量需求)。尤其是高分辨率SAR影像，價(jià)格較高，數(shù)據(jù)成本將十分顯著。

采用本發(fā)明技術(shù)方案，理論上只要差分干涉圖序列中存在至少6個(gè)短時(shí)間基線干涉對(duì)(由4幅SAR影像構(gòu)成)即可實(shí)現(xiàn)線性形變分量的解算，然后開(kāi)展長(zhǎng)時(shí)間基線差分干涉圖的形變相位梯度修正，促使更多的差分干涉圖(尤其是長(zhǎng)時(shí)間基線差分干涉圖)得到正確解纏，并參與形變分量的重新估算以及形變時(shí)間序列的建模和解算。因此，無(wú)需很高的SAR影像時(shí)間采樣率便可實(shí)現(xiàn)大梯度形變的有效提取。

上述說(shuō)明示出并描述了發(fā)明的若干優(yōu)選實(shí)施例，但如前所述，應(yīng)當(dāng)理解發(fā)明并非局限于本文所披露的形式，不應(yīng)看作是對(duì)其他實(shí)施例的排除，而可用于各種其他組合、修改和環(huán)境，并能夠在本文所述發(fā)明構(gòu)想范圍內(nèi)，通過(guò)上述教導(dǎo)或相關(guān)領(lǐng)域的技術(shù)或知識(shí)進(jìn)行改動(dòng)。而本領(lǐng)域人員所進(jìn)行的改動(dòng)和變化不脫離發(fā)明的精神和范圍，則都應(yīng)在發(fā)明所附權(quán)利要求的保護(hù)范圍內(nèi)。