本發(fā)明涉及一種基于兩景SAR強度影像的礦區(qū)地表大量級三維形變估計方法。

背景技術(shù)：

礦區(qū)地表在東西、南北和垂直方向的三維形變監(jiān)測對于提前評估和控制地下開采導致的潛在地質(zhì)災(zāi)害和建構(gòu)筑物損壞有著重要作用。

合成孔徑雷達干涉測量(interferometric synthetic aperture radar,InSAR)技術(shù)是一種新興的雷達遙感技術(shù)。該技術(shù)基于兩景SAR單視復數(shù)(single look complex,SLC)影像生成的單個InSAR干涉對，實現(xiàn)了地表雷達視線方向(line-of-sight,LOS)形變的大范圍、高精度、低成本、非接觸監(jiān)測。然而，由于InSAR技術(shù)僅能從單個InSAR干涉對中獲取一維LOS向形變，而非真正的地表三維形變。從而大大地制約了InSAR技術(shù)在礦區(qū)的應(yīng)用前景。為了獲取礦區(qū)地表三維形變，傳統(tǒng)的方法至少需要三個來自于不同雷達成像幾何學SAR數(shù)據(jù)的InSAR干涉對，即至少需要六景SLC影像。受限于目前較少的可用SAR衛(wèi)星數(shù)量，該方法在實際應(yīng)用中很難滿足。

發(fā)明專利ZL201210440875.4提出了一種利用單個InSAR干涉對獲取礦區(qū)地表三維形變場的方法，該發(fā)明有效地克服了傳統(tǒng)方法無法從單個InSAR干涉對中獲取礦區(qū)地表三維形變的局限，大大降低了傳統(tǒng)方法對于數(shù)據(jù)的苛刻要求，為礦區(qū)地質(zhì)災(zāi)害和建構(gòu)筑物風險評估提供了重要的參考依據(jù)。由于InSAR技術(shù)是基于兩景SLC影像的干涉相位的，所以當干涉相位失相關(guān)時，InSAR技術(shù)則無法準確地獲取地表雷達視線向形變。然而，大量級形變常常導致InSAR干涉相位失相關(guān)，從而致使發(fā)明專利ZL201210440875.4無法獲取礦區(qū)大量級三維形變。

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

本發(fā)明提供了一種基于兩景SAR強度影像的礦區(qū)地表大量級三維形變估計方法，其目的在于克服傳統(tǒng)基于單個InSAR干涉對的方法無法獲取礦區(qū)地表大量級三維形變的局限，利用礦區(qū)地表水平移動與下沉梯度的比例關(guān)系從SAR強度偏移量追蹤算法得到的礦區(qū)地表視線向和方位向形變中估計礦區(qū)地表沿著垂直、東西和南北方向的三維形變，大大地拓寬了SAR技術(shù)在礦區(qū)的應(yīng)用前景。

一種基于兩景SAR強度影像的礦區(qū)地表大量級三維形變估計方法，包括以下幾個步驟：

步驟1：利用SAR強度偏移量追蹤算法從兩景SAR強度影像中生成待監(jiān)測礦區(qū)地表雷達視線方向d_LOS和方位方向形變d_AZI；

待監(jiān)測的礦區(qū)的形變圖的總行數(shù)和總列數(shù)分別為n和m；

步驟2：基于礦區(qū)地表水平移動與下沉梯度的比例關(guān)系建立東西、南北方向的水平移動U_E、U_N與下沉W的函數(shù)關(guān)系；

步驟3：分別根據(jù)視線向形變d_LOS和方位向形變d_AZI與地表真實三維形變的投影關(guān)系建立礦區(qū)地表下沉與利用偏移量追蹤算法獲取的視線向與方位向形變之間的投影方程；

步驟4：按照步驟3構(gòu)建的投影方程建立待監(jiān)測礦區(qū)地表各像素點的下沉值與視線向形變的觀測方程；

由于礦區(qū)邊緣受地下開采影響較小，在待監(jiān)測礦區(qū)最后一行和最后一列的水平移動設(shè)置為0；

步驟5：求解待監(jiān)測礦區(qū)各像素點的下沉值；

步驟6：基于求解的下沉值并利用步驟2中所述的東西、南北方向的水平移動與下沉梯度的關(guān)系估計待監(jiān)測礦區(qū)地表在東西、南北方向的水平移動。

所述利用偏移量追蹤算法獲取的視線向與方位向形變中的偏移量追蹤算法屬于現(xiàn)有技術(shù)。

進一步地，所述基于礦區(qū)地表水平移動與下沉梯度的比例關(guān)系建立的東西、南北方向的水平移動U_E、U_N與下沉W的函數(shù)關(guān)系如下：

其中，i和j為地表點的像素坐標；W(i,j)表示地表點(i,j)處的下沉值，C_E和C_N分別示地表點(i,j)處在東西、南北方向的下沉梯度比例系數(shù)，b、H和β分別表示礦區(qū)水平移動系數(shù)、采深和主要影響角，b、H和β從監(jiān)測礦區(qū)獲??；R_E和R_N分別表示視線向和方位向形變圖在東西和南北方向的空間分辨率，R_E和R_N從視線向和方位向形變圖中獲取。

進一步地，所述礦區(qū)地表下沉與利用偏移量追蹤算法獲取的視線向與方位向形變之間的投影方程如下：

和/或

其中，為投影系數(shù)矩陣；

A₁₁(i,j)＝-A₁₂(i,j)-A₁₃(i,j)，A₂₁(i,j)＝cosθ-A₂₂(i,j)-A₂₃(i,j)；

A₁₂(i,j)＝sin(α-3π/2)C_E(i,j)，A₂₂(i,j)＝-b·C_E(i,j)·sin(α-3π/2)；

A₁₃(i,j)＝cos(α-3π/2)C_N(i,j)，A₂₃(i,j)＝-b·C_N(i,j)·cos(α-3π/2)；

α和θ表示雷達飛行方位角和入射角，從SAR強度影像的頭文件中獲得。

用方位向或者視線向監(jiān)測三維形變，即用任意一個投影方程來估計三維形變時，采用線性方程求解算法直接求解礦區(qū)下沉值；

兩者聯(lián)合使用時，采用最小二乘求解算法直接求解礦區(qū)下沉值，聯(lián)合使用時增加了多余觀測，使得下沉值的求解結(jié)果精度更高。

進一步地，所述待監(jiān)測礦區(qū)地表各像素點的下沉值與視線向形變的觀測方程如下：

其中，W_i和W_i+1分別表示待監(jiān)測的礦區(qū)形變圖中第i行和第i+1行像素中各像素的下沉值；

W_i＝[W(i,1),W(i,2),…,W(i,m)]^T，W_i+1＝[W(i+1,1),W(i+1,2),…,W(i+1,m)]^T；

d_LOSi和d_AZIi分別表示待監(jiān)測的礦區(qū)形變圖中第i行像素中各像素的視線方向和方位方向形變觀測值；

d_LOSi＝[d_LOS(i,1),d_LOS(i,2),…,d_LOS(i,m)]^T，d_AZIi＝[d_AZI(i,1),d_AZI(i,2),…,d_AZI(i,m)]^T；

C_i、C_i′、D_i及D_i′表示觀測系數(shù)矩陣；

C_i′＝diag[A₁₃(i,1) A₁₃(i,2) … A₁₃(i,m-1) 0]_m×m,

D_i′＝diag[A₂₃(i,1) A₂₃(i,2) … A₂₃(i,m-1) 0]_m×m

其中，待監(jiān)測的礦區(qū)形變圖中最后一行像素中各像素的視線方向和方位方向形變觀測值為D″＝diag(cosθ cosθ … cosθ)_m×m，diag表示對角陣。

進一步地，按照下述公式求解待監(jiān)測礦區(qū)各像素點的下沉值：

利用最小二乘法求解：

其中，E_i和L_i表示中間變量矩陣，為E_i的轉(zhuǎn)置矩陣，表示的逆矩陣；

其中，最后一行各像素的下沉值按照公式計算獲得：W_n＝(D″)^-1·d_LOSn。

有益效果

本發(fā)明公開了一種基于兩景SAR強度影像的礦區(qū)地表大量級三維形變估計方法，利用SAR強度偏移量追蹤算法從兩景SAR強度影像中生成礦區(qū)地表雷達視線向和方位向形變；基于礦區(qū)地表水平移動與下沉梯度的比例關(guān)系建立東西、南北方向水平移動與下沉之間的函數(shù)模型；分別建立礦區(qū)地表下沉與視線向和方位向形變的觀測方程組；利用最小二乘法估計求解礦區(qū)地表下沉；基于求解的地表下沉，使用礦區(qū)地表水平移動與下沉梯度的比例關(guān)系估計礦區(qū)地表在東西、南北方向的水平移動。本發(fā)明實現(xiàn)了僅利用兩景SAR強度影像估計礦區(qū)地表大量級三維形變，有效地克服了傳統(tǒng)方法無法從單個InSAR干涉對中獲取大量級三維形變的局限。此外，該方法數(shù)據(jù)制約少，監(jiān)測結(jié)果準確有效，大大拓寬了InSAR技術(shù)的應(yīng)用前景。本方法其對于拓寬SAR技術(shù)的應(yīng)用空間，降低礦區(qū)三維形變監(jiān)測成本有著重要意義。此外，其對于指導礦區(qū)安全生產(chǎn)、分析礦區(qū)沉降機理、預警礦區(qū)地表地質(zhì)災(zāi)害以及生態(tài)環(huán)境保護也起著重要作用。

附圖說明

圖1為本發(fā)明的所述方法的流程示意圖。

圖2為應(yīng)用本發(fā)明所述方法得到的形變圖，其中，圖(a)為視線向形變圖，圖(b)為方位向形變圖；

圖3為應(yīng)用本發(fā)明所述方法獲得的模擬的三維形變與實際監(jiān)測獲得的三維形變對比示意圖，其中，圖(a)為應(yīng)用本發(fā)明所述方法模擬獲得的垂直方向形變圖，圖(b)為圖(a)對應(yīng)的監(jiān)測獲得的形變圖，圖(c)為AA’剖面的垂直模擬形變值和監(jiān)測值的對比圖，圖(d)為應(yīng)用本發(fā)明所述方法模擬獲得的東西方向形變圖，圖(e)為圖(d)對應(yīng)的監(jiān)測獲得的形變圖，圖(f)為AA’剖面在東西方向上的模擬形變值和監(jiān)測值的對比圖，圖(g)為應(yīng)用本發(fā)明所述方法模擬獲得的南北方向形變圖，圖(h)為圖(g)對應(yīng)的監(jiān)測獲得的形變圖，圖(i)為BB’剖面在南北方向上的模擬形變值和監(jiān)測值得對比圖。

具體實施方式

下面將結(jié)合附圖1-3對本發(fā)明做進一步的說明。

如圖1所示，一種基于兩景SAR強度影像的礦區(qū)地表大量級三維形變估計方法，包括以下步驟：

步驟1：實施例利用快速拉格朗日分析軟件FLAC3D模擬生成采深為550m，煤層傾角為0°，采厚為5m，采空區(qū)尺寸為700m×150m的礦區(qū)地表三維形變場，如圖3中圖(a)、圖(d)、圖(g)所示。之后將生成的三維形變按照入射角為38°、飛行方位角為350°投影到雷達視線向和方位向(即用式(2)計算)；為了模擬偏移量追蹤技術(shù)導致的形變誤差，實施例在投影后的視線向和方位向形變中加入-0.1到0.1之間服從均勻分布的誤差，并以此模擬利用SAR強度偏移量追蹤算法從兩景SAR強度影像中獲取礦區(qū)地表沿著視線方向d_LOS和方位方向的形變d_AZI，其結(jié)果分別如圖2中(a)和(b)所示；

步驟2：基于礦區(qū)地表水平移動與下沉梯度的比例關(guān)系建立東西、南北方向的水平移動U_E、U_N與下沉W的函數(shù)關(guān)系

其中，i和j為地表點的像素坐標；W(i,j)表示地表點(i,j)處的下沉值，C_E和C_N分別示地表點(i,j)處在東西、南北方向的下沉梯度比例系數(shù)，b、H和tanβ分別表示礦區(qū)水平移動系數(shù)、采深和主要影響角正切，其中，b＝0.35、H＝550m和tanβ＝2.0；R_E和R_N分別表示視線向和方位向形變圖在東西和南北方向的空間分辨率；其中，R_E＝3m和R_N＝3m；

步驟3：視線向d_LOS和方位向形變d_AZI是地表真實三維形變按照雷達成像幾何學的投影：

其中，θ和α為雷達的入射角和飛行方位角，其中，θ＝38°，α＝350°。

將式(1)代入式(2)可得：

其中，為投影系數(shù)矩陣；

A₁₁(i,j)＝-A₁₂(i,j)-A₁₃(i,j)，A₂₁(i,j)＝cosθ-A₂₂(i,j)-A₂₃(i,j)；

A₁₂(i,j)＝sin(α-3π/2)C_E(i,j)，A₂₂(i,j)＝-b·C_E(i,j)·sin(α-3π/2)；

A₁₃(i,j)＝cos(α-3π/2)C_N(i,j)，A₂₃(i,j)＝-b·C_N(i,j)·cos(α-3π/2)；

步驟4：由于礦區(qū)邊緣受地下開采影響較小，因此假設(shè)在待監(jiān)測礦區(qū)最后一行和最后一列的水平移動為0，即：

其中，n和m為待監(jiān)測礦區(qū)的像素總行數(shù)和總列數(shù)，在實施例中，n＝301，m＝467。

根據(jù)式(3)和(4)建立礦區(qū)各像素下沉值與視線向和方位形變之間的觀測方程組，該觀測方程組以矩陣形式表示為：

其中，W_i和W_i+1分別表示待監(jiān)測的礦區(qū)形變圖中第i行和第i+1行像素中各像素的下沉值；

W_i＝[W(i,1),W(i,2),…,W(i,m)]^T，W_i+1＝[W(i+1,1),W(i+1,2),…,W(i+1,m)]^T；

d_LOSi和d_AZIi分別表示待監(jiān)測的礦區(qū)形變圖中第i行像素中各像素的視線方向和方位方向形變觀測值；

d_LOSi＝[d_LOS(i,1),d_LOS(i,2),…,d_LOS(i,m)]^T，d_AZIi＝[d_AZI(i,1),d_AZI(i,2),…,d_AZI(i,m)]^T；

C_i、C_i′、D_i及D_i′表示觀測系數(shù)矩陣；

C_i′＝diag[A₁₃(i,1) A₁₃(i,2) … A₁₃(i,m-1) 0]_m×m,

D_i′＝diag[A₂₃(i,1) A₂₃(i,2) … A₂₃(i,m-1) 0]_m×m

D″＝diag(cosθ cosθ … cosθ)_m×m，diag表示對角陣

步驟5：式(6)中有n未知數(shù)W_n和n個線性無關(guān)的方程，求解式(6)可得：

W_n＝(D″)^-1·d_LOSn (7)

將解出的W_n代入式(5)后，式(5)中有n未知數(shù)W_n-1和2n個方程，所以W_n-1的最小二乘為：

其中，E_i和L_i表示中間變量矩陣，為E_i的轉(zhuǎn)置矩陣，將解出的W_n-1代入式(5)，即可利用式(8)求解出W_n-2，重復該過程直至解出W₁，如圖3(b)所示；

步驟6：基于求解的下沉值利用式(1)估計待監(jiān)測礦區(qū)地表在東西、南北方向的水平移動，如圖3(e)和3(h)所示。

通過對比模擬的三維形變與監(jiān)測的三維形變(如圖3所示)發(fā)現(xiàn)，兩者吻合較好，其在垂直、東西和南北方向的均方根誤差分別為0.003、0.03和0.06m。該結(jié)果表明本發(fā)明是可行和可靠的。為了更進一步驗證本發(fā)明的應(yīng)用效果，實施例中在三維形變場中分別選取了AA’剖面的下沉和東西方向水平移動以及BB’剖面在南北方向的水平移動對比分析，結(jié)果亦表明：兩個剖面上的模擬值和計算值吻合非常好，從而說明本發(fā)明是可行的和可靠的。

以上所述僅為本發(fā)明的優(yōu)選實施例而已，并不用于限制本發(fā)明，對于本領(lǐng)域的技術(shù)人員來說，本發(fā)明可以有各種更改和變化。凡在本發(fā)明的精神和原則之內(nèi)，所作的任何修改、等同替換、改進等，均應(yīng)包含在本發(fā)明的保護范圍之內(nèi)。