技術特征：

1.一種適用于大梯度地表沉降監(jiān)測的時序差分雷達干涉方法，其特征在于，具體按照以下步驟實施：

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

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

步驟3，永久散射體PS和相干目標CT探測及點集合并，得到相干散射體CS點集，提取CS點集上的差分干涉相位，并以CS為節(jié)點構建不規(guī)則Delaunay三角網(wǎng)TIN；

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

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

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

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

步驟8，從所有原始差分干涉圖中減去線性形變相位分量，對殘差相位進行重新解纏，然后再將線性形變相位分量加回重新解纏后的相位中，重新執(zhí)行步驟5，得到更新后的線性形變速率和高程誤差；

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

2.如權利要求1所述的一種適用于大梯度地表沉降監(jiān)測的時序差分雷達干涉方法，其特征在于，步驟1具體實施方式為：排除受雨雪天氣以及積雪影響的SAR影像后進行任意干涉組合配對，假設有N+1幅SAR影像，通過任意干涉組合形成N(N+1)/2個干涉對，每個干涉對由主、從兩幅SAR影像構成，然后，根據(jù)每個干涉對主、從SAR影像的獲取時間計算該干涉對的時間基線即SAR影像獲取的時間差；根據(jù)主、從SAR影像的參數(shù)文件所記錄的SAR傳感器運行位置及相應的時間參數(shù)計算該干涉對的空間基線，主、從影像成像時SAR傳感器的空間距離，取該空間距離在垂直于SAR傳感器視線方向上的分量為空間基線，也即垂直基線。

3.如權利要求1所述的一種適用于大梯度地表沉降監(jiān)測的時序差分雷達干涉方法，其特征在于，步驟2中在進行時序差分干涉處理時，為降低時間失相干和空間失相干的影響，采用時間基線閾值和空間基線閾值方法排除時間基線和空間基線大于某給定的閾值的干涉對，當SAR影像獲取時間整體跨度小于2年時不考慮對時間基線進行限制。

4.如權利要求1所述的一種適用于大梯度地表沉降監(jiān)測的時序差分雷達干涉方法，其特征在于，步驟3具體實施方式為：PS探測采用振幅閾值和振幅離差指數(shù)ADI閾值雙閾值方法，對于SAR影像中某一特定的像元，其在N+1幅SAR影像中的振幅值直接從SAR振幅影像中提取，則其時序振幅均值和ADI值為：

$\stackrel{\‾}{a}=\frac{\underset{i=1}{\overset{N+1}{\Σ}}{a}_i}{N+1}---\left(1\right)$

$D_a=\frac{{\mathrm{\σ}}_a}{\stackrel{\‾}{a}},{\mathrm{\σ}}_a=\sqrt{\frac{1}{N+1}\underset{i=1}{\overset{N+1}{\Σ}}{(a_i-\stackrel{\‾}{a})}^2}---\left(2\right)$

其中，i為SAR影像索引號(序號)；a_i和分別為該像元在第i幅SAR影像中的振幅值及在所有影像中的振幅值的均值；D_a和σ_a分別為該像元的ADI值及時序振幅標準差。當該像元的振幅信息滿足如下條件時，認為其為PS：

$\left\{\begin{array}{c}{D}_a\≤0.25\\ \stackrel{\‾}{a}\≥\stackrel{\‾}{A}\±2{\mathrm{\σ}}_A,\stackrel{\‾}{A}=\frac{\underset{i=1}{\overset{N+1}{\mathrm{\Σ}}}{A}_i}{N+1},A_i=\frac{\underset{c=1}{\overset{C}{\mathrm{\Σ}}}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{a}_{cl}}{CL},{\mathrm{\σ}}_A=\sqrt{\frac{1}{N+1}\underset{i=1}{\overset{N+1}{\mathrm{\Σ}}}{\left(A_i-\stackrel{\‾}{A}\right)}^2}\end{array}\right.---\left(3\right)$

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

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

${\widehat{\mathrm{\γ}}}_l=\frac{\left|\underset{z=1}{\overset{Z}{\mathrm{\Σ}}}{i}_{{\mathrm{IM}}_l}^{\left(z\right)}{i_{{\mathrm{IS}}_l}^{*}}^{\left(z\right)}\right|}{\sqrt{\underset{z=1}{\overset{Z}{\mathrm{\Σ}}}|i_{{\mathrm{IM}}_l}^{\left(z\right)}{|}^2+\underset{z=1}{\overset{Z}{\mathrm{\Σ}}}|i_{{\mathrm{IS}}_l}^{\left(z\right)}{|}^2}},l=1,2,\mathrm{...},L---\left(4\right)$

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

$min\left({\widehat{\mathrm{\γ}}}_l\right)\>{\mathrm{\γ}}_{crit},l=1,2,...,L---\left(5\right)$

其中，min(·)表示取變量的最小值；γ_crit為判別某個像元是否為CT的相干系數(shù)閾值，取0.25至0.3；

當探測出所有的PS和CT后，將二者進行合并，并去除重復點，得到CS點集，最后從所有差分干涉圖中提取CS點上的相位值。

5.如權利要求1所述的一種適用于大梯度地表沉降監(jiān)測的時序差分雷達干涉方法，其特征在于，步驟4中相位解纏采用基于離散點的相位解纏方法對CS上的差分干涉相位進行解纏處理；對于某個CS點P，其在差分干涉圖中的相位和該點上的真實相位的關系為：

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

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

6.如權利要求1所述的一種適用于大梯度地表沉降監(jiān)測的時序差分雷達干涉方法，其特征在于，步驟5中在對CS上的差分干涉相位進行相位解纏后，根據(jù)步驟4中構建的Delaunay三角網(wǎng)重新計算相鄰CS點間的解纏相位差，此處，相鄰的CS點是指Delaunay三角網(wǎng)中三角形邊的端點，以任意一邊的兩端點P和Q為例，二者在第l個差分干涉圖中對應于的解纏相位為：

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

線性形變和高程誤差建模及解算的目的是對v_P和v_Q及δh_P和δh_Q進行估算，以P和Q上的解纏相位為觀測量進行網(wǎng)絡鄰域差分建模，對于P和Q，二者之間的網(wǎng)絡鄰域差分相位增量為：

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

以L幅差分干涉圖為例，對于P和Q所對應的網(wǎng)絡邊而言，列出L個與式(9)相同的方程，組成相應的線性方程組，將其表達為矩陣的形式有：

ΔΨ＝AX+W (10)

其中，

$X=\left[\begin{array}{c}{\mathrm{\Δv}}_{PQ}\\ {\mathrm{\Δ\δh}}_{PQ}\end{array}\right]---\left(12\right)$

A＝[κ,η] (14)

其中，

$\mathrm{\κ}=-\frac{4\mathrm{\π}}{\mathrm{\λ}}{\[B_1^T,B_2^T,...,B_l^T,...,B_L^T\]}^T---\left(15\right)$

$\mathrm{\η}=-\frac{4\mathrm{\π}}{\mathrm{\λ}}{\[\frac{B_1^{\⊥}}{\stackrel{\‾}{R}\sin \stackrel{\‾}{\mathrm{\θ}}},\frac{B_2^{\⊥}}{\stackrel{\‾}{R}\sin \stackrel{\‾}{\mathrm{\θ}}},\mathrm{...},\frac{B_l^{\⊥}}{\stackrel{\‾}{R}\sin \stackrel{\‾}{\mathrm{\θ}}},\mathrm{...},\frac{B_L^{\⊥}}{\stackrel{\‾}{R}\sin \stackrel{\‾}{\mathrm{\θ}}}\]}^T---\left(16\right)$

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

$\left[\begin{array}{c}\mathrm{\Δ}{\hat{v}}_{PQ}\\ \mathrm{\Δ}\mathrm{\δ}{\hat{h}}_{PQ}\end{array}\right]={\left(A^{T}A\right)}^{-1}{A}^T\mathrm{\Δ}\mathrm{\Ψ}---\left(17\right)$

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

${\hat{v}}_P-{\hat{v}}_Q=\mathrm{\Δ}{\hat{v}}_{PQ}+r_{PQ}---\left(18\right)$

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

$B\hat{V}=\mathrm{\Δ}\hat{V}+r---\left(19\right)$

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

$\hat{V}={\left(B^{T}PB\right)}^{-1}{B}^{T}P\mathrm{\Δ}\hat{V}---\left(20\right)$

其中，P為線性速率增量的權陣，通過弧段最小二乘殘差進行估算得到。

7.如權利要求1所述的一種適用于大梯度地表沉降監(jiān)測的時序差分雷達干涉方法，其特征在于，步驟6中采用基于差分干涉圖模擬的驗證方法；

使用估計所得到的年形變速率或累積形變量和相應的干涉基線參數(shù)模擬差分干涉圖，即將形變量轉換為形變相位，形變速率到相位的轉換公式為：

其中，為線性形變在第l幅差分干涉圖中對應的相位；v為通過步驟5中的數(shù)據(jù)處理所得的線性形變速率；

在大形變梯度區(qū)域，差分干涉圖中的形變相位梯度隨時間基線的增大而顯著提升，在差分干涉圖中表現(xiàn)為密集的條紋變化，若估計所得形變量符合或接近真實情況，使用形變結果模擬的差分干涉圖與原始差分干涉圖相似，則執(zhí)行步驟8；若估計所得形變量不符合真實情況，則模擬所得條紋數(shù)目與原始差分干涉圖會存在明顯區(qū)別，則執(zhí)行步驟7。

8.如權利要求1所述的一種適用于大梯度地表沉降監(jiān)測的時序差分雷達干涉方法，其特征在于，步驟8具體為：

在得到正確的線性形變后，從所有的原始差分干涉圖中減去線性形變所對應的相位分量，設此時所得到的正確的線性形變?yōu)?img id="icf0047" file="FDA0001207696740000075.GIF" wi="67" he="62" img-content="drawing" img-format="GIF" orientation="portrait" inline="no" />仍以第l幅差分干涉圖為例，其所對應的相位分量為：

此時，通過下式獲取扣除線性形變分量后的相位：

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

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

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

在得到修正的解纏相位后，以其為觀測值重新執(zhí)行步驟5中的解算流程，得到最終的線性形變速率和高程誤差，設最終的線性形變速率和高程誤差分別為和則二者所對應的相位為：

9.如權利要求1所述的一種適用于大梯度地表沉降監(jiān)測的時序差分雷達干涉方法，其特征在于，步驟9具體為：

在步驟8得到最終的線性形變速率和高程誤差結果后，首先保留并輸出線性形變速率結果，然后根據(jù)下式將線性形變和高程誤差相位從原始差分干涉相位中扣除，得到新的差分干涉相位值：

且有：

${\widehat{\mathrm{\φ}}}_{linear\_top\_crct,l}={\mathrm{\φ}}_{def\_n\ln r,l}+{\mathrm{\φ}}_{orb\_res,l}+{\mathrm{\φ}}_{atm,l}+{\mathrm{\φ}}_{noi,l}---\left(29\right)$

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

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

對于某個CS而言，其在第l幅扣除了高程誤差的解纏相位圖中的相位表達為：

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

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

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

此時，以t₀時刻為參考時刻，其余任意時刻t_n(n＝1,2,…,N)相對于t₀時刻的相位為待估參數(shù)，以扣除高程誤差后的解纏相位為觀測值，建立觀測方程并恢復每個時刻的累積形變相位。為得到符合實際物理意義的形變估計參數(shù)，假設相鄰兩幅SAR影像獲取時間間隔內(nèi)相位的變化符合線性累積趨勢，將對相位時間序列的求解轉化為對相位變化速率v_ph的求解，此時，待估參數(shù)為：

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

將式(35)表達為矩陣的形式有：

其中，B為L×N的矩陣，矩陣元素B[l,n]＝t_n-t_n-1(IS_l+1≤n≤IM_l)，其它元素值為零；

首先對B進行奇異值分解SVD處理，然后估算出最小范數(shù)意義下各時間段的平均相位速率，通過在時間域上進行積分即可恢復與SAR影像獲取時刻相對應的相位時間序列v_ph；B的奇異值分解為：

B＝USW^T (37)

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

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

其中，diag(·)代表對角矩陣，除對角線元素為σ_i外，其余元素值均為零；N-r+1為矩陣B的秩，r代表差分干涉圖子集的數(shù)量，對公式(36)進行解算得：

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

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

最后，從相位序列中扣除大氣延遲和軌道誤差相位序列，得到地表形變所對應的相位序列此時得形變時間序列為：

其中，D為與每幅SAR影像獲取時刻相對應的累積形變量所組成的向量；為與每幅SAR影像獲取時刻相對應的累積形變相位，即：

$D=\[d_{t_1},d_{t_2},...,d_{t_N}\]---\left(41\right)$