本發(fā)明涉及一種基于sar強(qiáng)度影像的礦區(qū)大量級三維時(shí)序形變估計(jì)方法及裝置。

背景技術(shù)：

礦區(qū)地表三維時(shí)序形變監(jiān)測對于理解礦區(qū)開采沉陷動態(tài)機(jī)理以及評估礦區(qū)潛在地質(zhì)災(zāi)害起著重要作用。偏移量追蹤(offsettracking，ot)技術(shù)能夠從兩景配準(zhǔn)的sar強(qiáng)度影像中獲取地表沿著雷達(dá)視線方向和方位方向的二維大量級形變(比如幾米或者幾十米)。2011年，casu等提出了一種名為ot-sbas(smallbaselinesubset)方法(參見文獻(xiàn)1casuf,manconia,pepea,etal.deformationtime-seriesgenerationinareascharacterizedbylargedisplacementdynamics:thesaramplitudepixel-offsetsbastechnique[j].ieeetransactionsongeoscienceandremotesensing,2011,49(7):2752-2763.)，該方法利用sar偏移量追蹤算法，實(shí)現(xiàn)了單個(gè)雷達(dá)成像幾何學(xué)sar影像集的地表二維(沿著視線向和方位向)大量級時(shí)序形變監(jiān)測。然而，該方法無法獲取傳統(tǒng)意義上的沿著垂直、東西和南北方向的三維大量級時(shí)序形變。為了克服該局限，raucoules等于2013年提出基于ot-sbas方法(參見文獻(xiàn)2raucoulesd,demichelem,maletjp,etal.time-variable3dgrounddisplacementsfromhigh-resolutionsyntheticapertureradar(sar).applicationtolavalettelandslide(southfrenchalps)[j].remotesensingofenvironment,2013,139:198-204.)和兩個(gè)顯著不同的雷達(dá)成像幾何學(xué)sar強(qiáng)度影像集，實(shí)現(xiàn)了地表大量級三維時(shí)序形變獲取方法。

然而，raucoules的方法存在著兩個(gè)明顯的局限：1)對于像礦區(qū)地表形變這樣的大量級、高度非線性的時(shí)序形變而言，兩個(gè)顯著不同雷達(dá)成像幾何學(xué)sar強(qiáng)度影像集必須是時(shí)間同步才能保證監(jiān)測的結(jié)果可靠。遺憾地是，受限于當(dāng)前較少的可用衛(wèi)星數(shù)量和太陽同步的軌道配置，該要求在實(shí)際生產(chǎn)中幾乎不可能滿足。2)偏移量追蹤算法獲取的形變觀測值常常含有粗差(通常定義為超過3倍中誤差的誤差)，而raucoules的方法使用傳統(tǒng)的等權(quán)最小二乘方法估計(jì)大量級三維時(shí)序形變。由于等權(quán)最小二乘方法穩(wěn)健性較差，所以，其估計(jì)的時(shí)序形變精度受粗差影響較大，從而大大地削弱了三維時(shí)序形變的精度。

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

本發(fā)明提供了一種基于sar強(qiáng)度影像的礦區(qū)大量級三維時(shí)序形變估計(jì)方法及裝置，其目的在于克服上述現(xiàn)有技術(shù)中的局限，利用實(shí)際應(yīng)用過程中容易獲取的單一雷達(dá)成像幾何學(xué)sar強(qiáng)度影像集，并結(jié)合aot-sap方法和三維形變速率估計(jì)值獲取礦區(qū)地表大量級三維時(shí)序形變。

一種基于sar強(qiáng)度影像的礦區(qū)大量級三維時(shí)序形變估計(jì)方法，包括以下步驟：

步驟1：根據(jù)sar數(shù)據(jù)特征和待監(jiān)測礦區(qū)地形設(shè)定sar強(qiáng)度影像對的時(shí)空基線閾值，對待監(jiān)測礦區(qū)單個(gè)雷達(dá)成像幾何學(xué)sar數(shù)據(jù)集生成時(shí)空基線小于時(shí)空基線閾值的sar強(qiáng)度影像對；

已有的方法至少需要兩個(gè)同步的不同雷達(dá)成像幾何學(xué)sar數(shù)據(jù)，本方案只需要一個(gè)雷達(dá)成像幾何學(xué)sar數(shù)據(jù)即可。比如，監(jiān)測某礦區(qū)從2016年1月到2017年1月的時(shí)間間隔為1個(gè)月的大量級時(shí)序形變，傳統(tǒng)方法需要24景(每月兩景)來自于兩個(gè)不同雷達(dá)成像幾何學(xué)的sar影像，而本發(fā)明只需要12景(每月1景)，成本大約下降了50％；

待監(jiān)測礦區(qū)單個(gè)雷達(dá)成像幾何學(xué)sar數(shù)據(jù)集的數(shù)量為m+1，生成的sar強(qiáng)度影像對的數(shù)量為g；

步驟2：利用aot-sap方法分別處理步驟1中生成的sar強(qiáng)度影像對，獲得待監(jiān)測礦區(qū)在垂直方向lw＝[lw1,lw2,…,lwg]、東西方向le＝[le1,le2,…,leg]和南北方向ln＝[ln1,ln2,…,lng]的三維多時(shí)域觀測值；

所述aot-sap方法是指一種基于兩景sar強(qiáng)度影像的礦區(qū)地表大量級三維形變估計(jì)方法，該方法是申請?zhí)枮?01710038785.5的中國專利，公開日為2017-05-31；

步驟3：設(shè)置相鄰sar影像期間礦區(qū)地表在垂直、東西和南北方向的三維形變速率，構(gòu)建三維多時(shí)域觀測值與三維形變速率之間的函數(shù)方程組；

其中，vw＝[vw1,vw2,…,vwm],ve＝[ve1,ve2,…,vem]和vn＝[vn1,vn2,…,vnm]分別為相鄰sar影像期間礦區(qū)地表在垂直方向、東西方向和南北方向的三維形變速率；

b為sar強(qiáng)度影像對輔影像和主影像獲取時(shí)間差的系數(shù)矩陣，其維度為g×m；

對于b的任意第k行，第imk個(gè)元素前的所有元素均為0，從第imk到第isk-1個(gè)元素，依次為：第isk-1個(gè)元素以后的所有元素均為0；

其中，imk和isk分別為生成第k個(gè)sar強(qiáng)度影像對的主、輔影像獲取時(shí)間索引，根據(jù)sar強(qiáng)度影像對的組成情況獲得；

步驟4：求解步驟3中函數(shù)方程組中的三維形變速率估計(jì)值：

和

步驟5：利用步驟4獲得的三維形變率估計(jì)值估計(jì)待監(jiān)測礦區(qū)在垂直方向、東西方向和南北方向上相對于第一景sar影像獲取時(shí)刻的三維大量級時(shí)序形變(w(tk)，e(tk)，n(tk))；

其中，tk表示sar影像的獲取時(shí)間，k＝1,2,…,m+1，w(t0)＝e(t0)＝n(t0)≡0。

進(jìn)一步地，采用最小二乘法求解步驟3中函數(shù)方程組中的三維形變速率估計(jì)值。

進(jìn)一步地，采用穩(wěn)健估計(jì)法求解步驟3中函數(shù)方程組中的三維形變速率估計(jì)值。

一種基于sar強(qiáng)度影像的礦區(qū)大量級三維時(shí)序形變估計(jì)裝置，包括：

sar強(qiáng)度影像對生成單元，用于根據(jù)sar數(shù)據(jù)特征和待監(jiān)測礦區(qū)地形設(shè)定sar強(qiáng)度影像對的時(shí)空基線閾值，對待監(jiān)測礦區(qū)單個(gè)雷達(dá)成像幾何學(xué)sar數(shù)據(jù)集生成時(shí)空基線小于時(shí)空基線閾值的sar強(qiáng)度影像對；

三維多時(shí)域觀測值獲取單元，利用aot-sap方法分別處理生成的sar強(qiáng)度影像對，獲得待監(jiān)測礦區(qū)在垂直方向lw＝[lw1,lw2,…,lwg]、東西方向le＝[le1,le2,…,leg]和南北方向ln＝[ln1,ln2,…,lng]的三維多時(shí)域觀測值；

函數(shù)方程組構(gòu)建單元，通過設(shè)置相鄰sar影像期間礦區(qū)地表在垂直、東西和南北方向的三維形變速率，構(gòu)建三維多時(shí)域觀測值與三維形變速率之間的函數(shù)方程組；

三維形變速率估計(jì)值計(jì)算單元，計(jì)算函數(shù)方程組中的三維形變速率估計(jì)值；

三維大量級時(shí)序形變求解單元，利用三維形變率估計(jì)值估計(jì)待監(jiān)測礦區(qū)相對于第一景sar影像獲取時(shí)刻的三維大量級時(shí)序形變(w(tk)，e(tk)，n(tk))；

所述三維大量級時(shí)序形變求解公式如下：

其中，tk表示sar影像的獲取時(shí)間，k＝1,2,…,m+1，w(t0)＝e(t0)＝n(t0)≡0。

進(jìn)一步地，所述三維形變速率估計(jì)值計(jì)算單元采用穩(wěn)健估計(jì)法進(jìn)行計(jì)算。

有益效果

本發(fā)明提供了一種基于sar強(qiáng)度影像的礦區(qū)大量級三維時(shí)序形變估計(jì)方法及裝置，該方法首先基于設(shè)定的時(shí)空基線閾值從單個(gè)雷達(dá)成像幾何學(xué)sar強(qiáng)度影像集中生成sar強(qiáng)度影像對；利用已有的aot-sap方法分別處理各sar強(qiáng)度影像對，獲得礦區(qū)地表多時(shí)域三維形變觀測值；分別建立相鄰sar影像期間礦區(qū)地表三維形變速率與多時(shí)域三維形變觀測值之間的函數(shù)模型；利用穩(wěn)健估計(jì)求解礦區(qū)地表三維大量級時(shí)序形變。該裝置結(jié)構(gòu)簡單，實(shí)現(xiàn)方便；相較于現(xiàn)有技術(shù)本方案的優(yōu)點(diǎn)主要體現(xiàn)在以下幾點(diǎn)：

第一、該方法首次實(shí)現(xiàn)了僅利用單個(gè)雷達(dá)成像幾何學(xué)sar強(qiáng)度影像集的礦區(qū)地表大量級三維時(shí)序形變獲取，有效地克服了傳統(tǒng)方法對于sar數(shù)據(jù)的苛刻要求，節(jié)約了礦區(qū)大量級三維形變sar監(jiān)測的成本(約50％)，提高了三維時(shí)序形變的精度。

第二、不同于利用合成孔徑雷達(dá)干涉測量(interferometricsyntheticapertureradar，insar)估計(jì)礦區(qū)地表三維時(shí)序形變(比如專利：“基于單個(gè)雷達(dá)成像幾何學(xué)sar影像的礦區(qū)地表三維形變估計(jì)”(申請?zhí)枺?01610546270.1))，本發(fā)明基于偏移量追蹤技術(shù)，而非insar技術(shù)，所以其能基于單個(gè)雷達(dá)成像幾何學(xué)sar數(shù)據(jù)估計(jì)礦區(qū)地表大量級(幾米甚至幾十米)三維形變，而基于insar技術(shù)的則不能。

第三、本方法其對于拓寬sar應(yīng)用空間，指導(dǎo)礦區(qū)安全生產(chǎn)、分析礦區(qū)沉降機(jī)理、預(yù)警礦區(qū)地表地質(zhì)災(zāi)害以及礦區(qū)生態(tài)環(huán)境保護(hù)和治理也起著重要作用。

附圖說明

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

圖2模擬的三維時(shí)序形變圖；

圖3本實(shí)施例的設(shè)計(jì)矩陣；

圖4本發(fā)明估計(jì)的三維時(shí)序形變與模擬值之間差值的統(tǒng)計(jì)直方圖，其中，(a)為垂直方向，(b)為東西方向，(c)為南北方向。

具體實(shí)施方式

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

本實(shí)施例利用概率積分法(一種開采沉陷預(yù)計(jì)數(shù)學(xué)模型)模擬了地下一條工作面開采導(dǎo)致的礦區(qū)地表在八個(gè)時(shí)間點(diǎn)的三維大量級時(shí)序形變(時(shí)間間隔30天)，如圖2所示。將任意兩個(gè)時(shí)間點(diǎn)的時(shí)序形變相減，從而生成了28個(gè)多時(shí)域三維形變觀測值。為了使模擬的多時(shí)域三維形變觀測值更接近真實(shí)情況，本發(fā)明在每個(gè)觀測值中加入了均值為0，均方根誤差為0.1m的高斯誤差。之后，對每一個(gè)像素的三維多時(shí)域觀測值隨機(jī)加入或者減去一個(gè)10到20倍均方根誤差(即從1到2m)的粗差，從而模擬了含粗差的多時(shí)域三維形變觀測值。

如圖1所示，一種基于sar強(qiáng)度影像的礦區(qū)地表大量級三維時(shí)序形變估計(jì)方法，包括以下步驟：

步驟1：令覆蓋待監(jiān)測礦區(qū)的單個(gè)雷達(dá)成像幾何學(xué)sar強(qiáng)度影像的數(shù)量為m+1(本實(shí)施例中m＝7)，其時(shí)間按照先后順序可表示為[t0,t1,…,tm]；根據(jù)影像特征(比如波長、分辨率)和研究區(qū)域地形等因素設(shè)定時(shí)空基線閾值，并從m+1景sar強(qiáng)度影像集中生成時(shí)空基線小于設(shè)定時(shí)空基線閾值的可用sar強(qiáng)度影像對(其數(shù)量表示為g，本實(shí)施例中g(shù)＝28)；

步驟2：利用aot-sap方法(alternativeoffsettracking-singleamplitudepair，專利申請?zhí)枺?01710038785.5由楊澤發(fā)等提出的基于兩景sar強(qiáng)度影像數(shù)據(jù)和偏移量追蹤技術(shù)的礦區(qū)地表大量級三維形變獲取的方法)分別處理步驟1中生成的sar強(qiáng)度影像對，從而獲得待監(jiān)測礦區(qū)在垂直lw＝[lw1,lw2,…,lwg]、東西le＝[le1,le2,…,leg]和南北方向ln＝[ln1,ln2,…,lng]的三維多時(shí)域觀測值；

步驟3：令vw＝[vw1,vw2,…,vwm],ve＝[ve1,ve2,…,vem]和vn＝[vn1,vn2,…,vnm]分別為相鄰sar影像期間礦區(qū)地表在垂直、東西和南北方向的三維形變速率，所以其于三維多時(shí)域觀測值lw、le和ln之間的函數(shù)方程組可表示為：

式中，b為方程組設(shè)計(jì)矩陣，其取決于sar強(qiáng)度對的分布，在本實(shí)施例中，系數(shù)矩陣b如圖3所示；

式(1)中每一個(gè)方程組均含有g(shù)個(gè)觀測方程和m個(gè)未知數(shù)(即平均形變速率)，所以若g≥m，則式(1)中每個(gè)方程組均為正定或超定方程組，其未知數(shù)能則較為容易地被求解。在實(shí)際中，該條件是很容易滿足的，主要原因是偏移量追蹤算法對于干涉相位不敏感，所以其受相位失相關(guān)的影響較小，從而可形成大量的sar強(qiáng)度影像對；換句話說，在實(shí)際中，觀測方程個(gè)數(shù)g通常大于未知數(shù)個(gè)數(shù)m(即g>m)，且式(1)中方程組為超定方程組；

步驟4：采用穩(wěn)健估計(jì)法求解步驟3中函數(shù)方程組中的三維形變速率估計(jì)值；

對于一個(gè)超定方程組，最常用的求解方法是等權(quán)最小二乘。然而，由于等權(quán)最小二乘的影響函數(shù)是無界的，所以其對于粗差的抵抗能力較弱。換句話說，最小二乘解對于粗差較為敏感，觀測值中的一個(gè)粗差可能大大地降低求解參數(shù)的精度，甚至得出錯(cuò)誤的解。然而，式(1)中的三維形變觀測值是基于偏移量追蹤算法估計(jì)的二維視線向和方位向形變，并利用aot-sap方法獲得，其含有粗差的概率較大。主要因?yàn)橐韵聝蓚€(gè)因素：1)偏移量追蹤算法精度主要取決于影像分辨率和設(shè)置的搜索窗口，因此，偏移量追蹤算法估計(jì)的二維視線向和方位向形變中極有可能含有粗差；2)在aot-sap估計(jì)三維形變時(shí)，偏移量追蹤算法獲取的二維形變誤差通常被放大到二維水平移動觀測值le和ln中，從而增大了三維形變觀測值中含有粗差的可能性。鑒于穩(wěn)健估計(jì)具有較強(qiáng)的抵御粗差的能力，因此，本實(shí)例使用穩(wěn)健估計(jì)方法估計(jì)礦區(qū)地表三維形變速率vw，ve和vn。

本發(fā)明以穩(wěn)健估計(jì)中的m估計(jì)為例詳細(xì)說明求解過程。m估計(jì)是一種選權(quán)迭代法，通過對含有粗差的觀測值加較小的權(quán)值(甚至零權(quán)值)，從而抑制粗差對于方程組解的貢獻(xiàn)。利用m估計(jì)，式(1)中的三個(gè)方程組的穩(wěn)健解可通過迭代式(2)求解；

直到相鄰兩次迭代解的最大差值小于一個(gè)很小的容差(比如10^-6)。在式(2)中，和為m估計(jì)的第i+1次迭代解；和為m估計(jì)中第i+1次迭代的權(quán)對角矩陣；本實(shí)施例通過igg-i加權(quán)函數(shù)確定權(quán)對角矩陣的權(quán)值。以第k個(gè)(k＝1,2,…,g)垂直下沉觀測值的權(quán)為例：

式中，a和b為設(shè)定閾值，其通常取a＝1.5和b＝2.5；其中為式(1)中垂直方向方程組的殘差，即為第i次單位權(quán)中誤差，med表示求中位數(shù)；

步驟5：估計(jì)礦區(qū)地表在垂直w＝[w(t1),w(t2),…,w(tm)]，東西e＝[e(t1),e(t2),…,e(tm)]和南北方向n＝[n(t1),n(t2),…,n(tm)]的三維時(shí)序形變：

式中tk(k＝1,2,…,m+1)表示sar影像(其總數(shù)量為m+1)的獲取時(shí)間,和為步驟4中估計(jì)的三維時(shí)序形變速率的穩(wěn)健估計(jì)值。

通過對比模擬的三維大量級時(shí)序形變與估計(jì)的三維時(shí)序形變發(fā)現(xiàn)(兩者的差值統(tǒng)計(jì)圖見圖4)，兩者吻合較好，其在垂直、東西和南北方向的均方根誤差均為0.062m(略小于模擬的均方根誤差為0.1m的高斯誤差)，該精度能夠滿足礦區(qū)大量級時(shí)序形變監(jiān)測精度需求。為了更進(jìn)一步驗(yàn)證本發(fā)明的應(yīng)用效果，實(shí)施例中還利用傳統(tǒng)的等權(quán)最小二乘估計(jì)了時(shí)序三維形變，結(jié)果表明：穩(wěn)健估計(jì)的結(jié)果精度相對于等權(quán)最小二乘估計(jì)的結(jié)果(其均方根誤差大約為0.144m)提高了約56.9％。通過本實(shí)施例說明：本發(fā)明不僅實(shí)現(xiàn)了基于單個(gè)雷達(dá)成像幾何學(xué)sar強(qiáng)度數(shù)據(jù)集的礦區(qū)地表大量級三維時(shí)序形變獲取，同時(shí)還大大地提高了獲取形變的精度。

一種基于sar強(qiáng)度影像的礦區(qū)大量級三維時(shí)序形變估計(jì)裝置，包括：

sar強(qiáng)度影像對生成單元，用于根據(jù)sar數(shù)據(jù)特征和待監(jiān)測礦區(qū)地形設(shè)定sar強(qiáng)度影像對的時(shí)空基線閾值，對待監(jiān)測礦區(qū)單個(gè)雷達(dá)成像幾何學(xué)sar數(shù)據(jù)集生成時(shí)空基線小于時(shí)空基線閾值的sar強(qiáng)度影像對；

三維多時(shí)域觀測值獲取單元，利用aot-sap方法分別處理生成的sar強(qiáng)度影像對，獲得待監(jiān)測礦區(qū)在垂直方向lw＝[lw1,lw2,…,lwg]、東西方向le＝[le1,le2,…,leg]和南北方向ln＝[ln1,ln2,…,lng]的三維多時(shí)域觀測值；

函數(shù)方程組構(gòu)建單元，通過設(shè)置相鄰sar影像期間礦區(qū)地表在垂直、東西和南北方向的三維形變速率，構(gòu)建三維多時(shí)域觀測值與三維形變速率之間的函數(shù)方程組；

三維形變速率估計(jì)值計(jì)算單元，計(jì)算函數(shù)方程組中的三維形變速率估計(jì)值；

三維大量級時(shí)序形變求解單元，利用三維形變率估計(jì)值估計(jì)待監(jiān)測礦區(qū)相對于第一景sar影像獲取時(shí)刻的三維大量級時(shí)序形變(w(tk)，e(tk)，n(tk))；

所述三維大量級時(shí)序形變求解公式如下：

其中，tk表示sar影像的獲取時(shí)間，k＝1,2,…,m+1，w(t0)＝e(t0)＝n(t0)≡0。

在本實(shí)例中三維形變速率估計(jì)值計(jì)算單元采用穩(wěn)健估計(jì)法進(jìn)行計(jì)算。

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