本發(fā)明涉及基于遙感影像的大地測量領域中合成孔徑雷達干涉測量(insar，synthetic?aperture?radar?interferometry)技術和全球導航衛(wèi)星系統(tǒng)(gnss，globalnavigation?satellite?system)技術，特別涉及一種時空kalman模型融合insar和gnss數(shù)據(jù)重建礦區(qū)地表形變監(jiān)測方法。

背景技術：

1、高強度的煤礦開采給礦區(qū)地表、地質結構和周圍環(huán)境帶來巨大影響，甚至引發(fā)一系列如滑坡、泥石流以及崩塌等地質災害問題，嚴重影響了當?shù)厝嗣竦纳敭a(chǎn)安全。形變是地質災害發(fā)生最直觀的前兆反映，因此采用科學的方法對采空區(qū)的地表形變進行監(jiān)測，對準確評估和預防礦區(qū)地質災害極為重要。在過去的20年，合成孔徑雷達干涉測量(insar，interferometric?synthetic?aperture?radar)技術以高空間分辨率、高精度和大范圍快速監(jiān)測的獨特優(yōu)勢，迅速在礦區(qū)、滑坡、地震、大氣、凍土和城市基礎設施等各個領域取得了巨大的突破和成就，成為目前地質災害監(jiān)測的最有效手段之一。然而insar形變監(jiān)測的時間分辨率卻受限于sar衛(wèi)星傳感器的重訪周期，一般為6-14天，導致其監(jiān)測的形變具有較大的滯后性。在面對礦區(qū)、滑坡等突發(fā)性強、變形梯度大的形變時無法及時反映形變的態(tài)勢，這極大限制了insar技術的發(fā)展和應用。而全球導航衛(wèi)星系統(tǒng)(gnss,global?navigationsatellite?system)具有實時、高精度和長時間連續(xù)觀測等優(yōu)點，已被廣泛應用于各類地質災害和工程應用的安全監(jiān)測。但是gnss的成本較高，只能在個別重要位置監(jiān)測，無法進行大面積高密度的布設，因此其監(jiān)測結果的空間分辨率較低，難以反映監(jiān)測區(qū)域形變的空間變化，尤其在形變梯度變化大的礦區(qū)。

2、時空kalman濾波是一種充分顧及了時空相關性的kalman濾波方法，在時空形變分析領域展示了極大的優(yōu)越性，能對形變監(jiān)測數(shù)據(jù)進行時空動態(tài)濾波、時空插值以及時空預測，尤其時空kalman中的時空隨機效應模型(stre，spatio-temporalrandom?effects)還可以對數(shù)據(jù)進行降維，大大提高了時空大數(shù)據(jù)的計算效率。因此有學者提出結合insar形變高空間分辨率和gnss形變高時間分辨率的優(yōu)勢，采用時空kalman濾波方法融合insar和gnss的監(jiān)測結果重建地表高時空分辨率的形變，并將其應用于南加州等斷層形變的恢復中，取得了較好的效果。然而這些研究中依舊存在兩個問題：一是在利用stre模型在融合過程中通過布設空間基來降維，以此達到快速解算的目的，且越少的空間基降維效果越好。而目前空間基的布設主要在多個尺度均勻布設，且依靠人為的經(jīng)驗指導，雖然可以在斷層這樣的均勻形變中滿足使用需求，卻不適用于形變不均勻且梯度大的礦區(qū)。二是在insar和gnss融合生成高時空分辨率的地表形變結果時，stre模型融合的形變精度與gnss的數(shù)量緊密相關，gnss數(shù)量越多則融合的形變精度越高，gnss數(shù)量越少則融合結果主要由insar形變結果主導。但是礦區(qū)通常無法大量布設gnss站點，導致已有的時空kalman融合方法不適用于時空形變變化復雜的礦區(qū)。因此面對目前礦區(qū)形變監(jiān)測的需求，考慮如何優(yōu)化時空kalman融合過程中礦區(qū)形變的空間建模方法以及解決礦區(qū)gnss站點數(shù)量少導致融合形變精度低的問題是目前急需解決的難題。

技術實現(xiàn)思路

1、本發(fā)明的目的在于克服目前基于時空kalman方法融合insar和gnss形變重建地表高時空分辨率的形變過程中難以適用于礦區(qū)形變梯度大、gnss數(shù)量少的缺陷，提出一種時空kalman模型融合insar和gnss數(shù)據(jù)重建礦區(qū)地表形變監(jiān)測方法。該方法根據(jù)礦區(qū)空間形變特征自適應構建空間模型，提高了算法的自動性和效率。同時附加礦區(qū)先驗形變的時間特征約束，不僅減少了傳統(tǒng)stre融合方法對gnss數(shù)據(jù)的依賴，而且進一步提高了形變?nèi)诤系木取?/p>

2、本發(fā)明首先提供了如下的技術方案：

3、一種時空kalman模型融合insar和gnss數(shù)據(jù)重建礦區(qū)地表形變監(jiān)測方法，其包括以下步驟：

4、s1、insar和gnss數(shù)據(jù)處理：對獲取的insar干涉圖進行濾波、解纏、大范圍的趨勢去除、大氣校正和時空濾波之后進行地理編碼，然后基于sbas-insar技術解算研究區(qū)域的平均形變速率和時間序列形變結果；將gnss三維形變監(jiān)測結果轉換到insar視線向；最后將gnss形變與insar形變進行時空基準統(tǒng)一。

5、s2、建立融合insar和gnss形變的時空kalman模型。首先建立stre模型，然后根據(jù)獲取的insar和gnss數(shù)據(jù)建立融合的空間模型、時間模型，形成時空相關的kalman融合模型。

6、s3、根據(jù)礦區(qū)的空間形變特征針對大尺度主要形變構建自適應的空間模型：以insar平均形變速率為基準，首先在全圖布設一層均勻空間基來捕捉全局的光滑形變；然后計算第一層空間基建模后的殘差，識別該殘差中值較大的部分，并對其進行形態(tài)學的膨脹操作，對獲得的不規(guī)則區(qū)域進行第二層空間基布設用于捕捉局部形變。接著計算第一層和第二層總的空間基建模后的殘差，識別該殘差中值較大的部分，并對其進行形態(tài)學的膨脹操作，對獲得的不規(guī)則區(qū)域進行第三層空間基布設用于捕捉局部小形變。最后合并insar三層空間基和gnss空間基，完成大尺度主要形變的自適應空間建模。

7、s4、針對礦區(qū)的時間形變特征，引入礦區(qū)先驗信息，為狀態(tài)量附加logistic模型約束，提高沒有insar觀測時間的狀態(tài)量的估計精度，同時減少對gnss數(shù)據(jù)的依賴。最后構建形變時空特征約束的融合insar和gnss形變的時空kalman模型。

8、s5、求解融合過程中的模型參數(shù)：首先通過半變異函數(shù)擬合、方差分量估計確定insar小尺度的細微形變的方差、insar觀測噪聲的方差以及gnss觀測噪聲的方差；然后通過最小二乘擬合、向前濾波和向后平滑確定融合過程中的模型參數(shù)；最后對整個融合過程用em迭代估計提高模型參數(shù)的估計精度。

9、s6、重建地表高時空分辨率的形變結果：根據(jù)上述步驟求解的空間基、狀態(tài)量和小尺度的細微形變重建恢復地表高時空分辨率的形變結果，并對形變結果進行精度評估。

10、根據(jù)本發(fā)明的一些優(yōu)選實施方式，所述s1包括：

11、s11、對獲取的雷達影像，經(jīng)過配準、小基線構網(wǎng)、干涉、去平地和去地形、相位濾波、相位解纏、大范圍的趨勢去除、大氣校正和時空濾波之后進行地理編碼后形成干涉圖集。值得注意的是由于后續(xù)沒有對形變中的趨勢進行建模，因此這一步必須要去除干涉圖中的趨勢項。

12、s12、假設研究區(qū)域一共有n景sar影像，根據(jù)時空基線閾值生成m幅干涉圖，且只有一個小基線集，則由短基線集(small?baseline?subsets,sbas)方法可以構建如下方程：

13、

14、對上式由最小二乘可得：

15、

16、其中是平差的最終形變序列結果；a是上式中的系數(shù)矩陣；p是權矩陣，可根據(jù)干涉圖的相干性等確定；l是干涉圖相位。

17、獲得形變序列之后可以采用最小二乘擬合求得整個時間段的平均形變速率。

18、s13、由于gnss有三個方向上的形變，而insar一般只能獲取到雷達視線向(los，line?ofsight)的形變，因此需要將gnss的三維方向轉換到雷達視線向：

19、

20、其中表示gnss三維形變轉換到los向后的形變，de、dn和du分別表示gnss點的東西、南北和垂直向的形變，θ表示衛(wèi)星的雷達局部入射角，α表示雷達方位角。

21、s14、將insar點和gnss點轉換到同一個坐標框架。

22、在空間上，找到距離gnss點最近的insar點作為gnss的同名點。

23、之后選擇其中一個穩(wěn)定的gnss點對其余gnss點和insar的時間維和空間維進行校正，從而實現(xiàn)兩者的時空基準統(tǒng)一。

24、根據(jù)本發(fā)明的一些優(yōu)選實施方式，所述s2包括：

25、s21、對于一系列時空形變數(shù)據(jù)z(s，t)，則可將其表示為：

26、z(s，t)＝y(tǒng)(s，t)+ω(s，t)

27、其中s＝1,2，..，n和t＝1,2，...，t分別表示區(qū)域d中點的空間位置和觀測時間，y(s，t)表示時空觀測數(shù)據(jù)的真值，ω(s，t)～n(0，σn2)表示時空觀測數(shù)據(jù)的噪聲。

28、進一步時空觀測數(shù)據(jù)的真值y(s，t)可以表示為如下形式：

29、y(s，t)＝vt(s)+ξt(s)

30、其中vt(s)表示空間上大尺度的主要形變，ξt(s)表示空間上小尺度的細微形變。注意如果形變中含有趨勢項，也可以對趨勢項進行空間建模。由于礦區(qū)形變監(jiān)測中一般會提前去除趨勢，因此不存在趨勢項形變。

31、進一步空間上大尺度的主要形變vt(s)一般建模為空間場st(s)和時變狀態(tài)量xt的線性組合進行描述，即：

32、

33、其中st(s)＝[s1，t(s)，s2，t(s),...，sr,t(s)]t，表示建立的r×n維的空間基函數(shù)，xt＝[x1，t(·)，x2，t(·),...，xr,t(·)]t表示r×m的零均值高斯隨機狀態(tài)量。注意基函數(shù)可以隨時間變化，但是對于一個礦區(qū)而言可以通過平均形變速率非常輕易捕捉形變區(qū)域，同時也為了減小模型的復雜度，因此將其建模為不隨時間變化的函數(shù)，即s(s)＝st(s)＝[s1(s)，s2(s),...，sr(s)]t。

34、基函數(shù)選擇在時空分析中常用的bisquare函數(shù)，其表達如下：

35、

36、其中，sl，j(s)表示第l層空間基中第j個點在空間s上的空間基函數(shù)值，j＝1,2,...，r；pl,j表示第l層空間基中第j個點的空間基位置；dl為第l層空間中空間變異的變程(一般取第l層尺度下最短空間基距離的1.5倍)。

37、s22、現(xiàn)在將狀態(tài)量xt在時間上進行演化構成狀態(tài)轉移方程：

38、xt＝hxt-1+et

39、其中，h是r×r的狀態(tài)轉移矩陣，是狀態(tài)轉移過程中的噪聲。通過以上所述可以建立stre模型的最終形式：

40、

41、從以上分析可知，stre的模型參數(shù)為小尺度的細微形變ξt(s)、狀態(tài)量xt、狀態(tài)轉移矩陣h和狀態(tài)轉移過程中的方差

42、s23、假設在研究區(qū)域d內(nèi)，共獲得了n1個insar點，其有m1個觀測日期t1；同時共獲得了n2個gnss點，其有m2個觀測日期假設待重建的所有時刻為t＝t2，insar與gnss的時間對應關系為id，則有t1＝t2(id)＝t(id)。令n＝n1+n2，m＝m2。

43、則根據(jù)以上s21所述stre建立的融合insar和gnss的時空隨機效應模型：

44、

45、其中各個變量的含義分別如下表1：

46、表1變量的含義

47、

48、注意在融合過程中insar和gnss的小尺度細微形變的方差是相同的。

49、最后同步驟s22，將狀態(tài)量xt在時間上進行演化構成狀態(tài)轉移方程，形成最終融合的時空kalman模型。

50、根據(jù)本發(fā)明的一些優(yōu)選實施方式，所述s3包括：

51、s31、在步驟s23中一般為了充分捕捉形變的空間變化，會在多個尺度上布設多層空間基，這相當于多個變異函數(shù)的套合。但是形變在空間上并不是均勻變化的，因此需要對空間基的位置布設做一些自適應的調整。

52、首先根據(jù)insar平均形變速率在整個形變場布設一層均勻的空間基，然后計算這層均勻空間基建模剩余的殘余形變：

53、

54、其中，v1表示布設第一層空間基建模后的殘差，za表示insar平均形變速率，s1表示第一層布設的r1個空間基。

55、s32、識別v1中值較大的部分，在該過程中會出現(xiàn)一些孤島和面積較小的區(qū)域。因此需要先刪除面積小的區(qū)域，同時進行形態(tài)學的膨脹操作，填補其中空洞的同時對邊界進行一定擴充，獲得第一層空間基未能充分捕捉形變的區(qū)域。

56、s33、對獲得的不規(guī)則區(qū)域采用更高密度布設第二層均勻空間基(假設為r2個)，并與第一次的空間基合并計算新的殘差：

57、

58、用s1,2表示第一次和第二次布設的r1+r2個空間基，則s1,2＝[s1,s2]t。

59、s34、由于礦區(qū)的形變梯度比較大，因此兩層空間基依然不能完全捕捉所有形變，因此還需要在此基礎上布設間距更小的第三層空間基。與步驟s32相似，識別v2中值較大的部分，在該過程中會出現(xiàn)一些孤島和面積較小的區(qū)域。因此需要先刪除面積小的區(qū)域，同時進行形態(tài)學的膨脹操作，填補其中空洞的同時對邊界進行一定擴充，獲得第一層和第二層空間基未能充分捕捉的區(qū)域。

60、s35、對獲得的不規(guī)則區(qū)域布設第三層層均勻空間基，并與前兩次的空間基合并計算新的殘差：

61、

62、用s3表示第三次布設的r3個空間基，則sins＝[s1,s2,s3]t，共r＝r1+r2+r3個空間基。

63、在上述過程中，只需指定每層空間基的密度，即可獲取整個區(qū)域自適應的空間基位置和數(shù)量。這不僅減少了人為的干預，提高了算法的自動性，還降低了空間基的數(shù)量，提升了算法的效率。特別是在進行大范圍形變計算時，可大幅度節(jié)省時間和空間成本。

64、值得注意的是一般情況下三次布設之后，殘差都會非常小，可以滿足后續(xù)計算。但是當形變梯度非常大時三層尺度的空間基可能依舊無法充分捕捉所有形變，這時可以在此三層的基礎上布設更多尺度的空間基以進一步提高空間建模的精度。但是過多的空間基會耗費大量的計算時間，且可能造成過擬合現(xiàn)象，因此需要折中選擇。

65、s36、用所有gnss點來計算gnss的空間基sgnss。至此分別獲得了insar和gnss的空間基。

66、根據(jù)本發(fā)明的一些優(yōu)選實施方式，所述s4包括：

67、s41、由于礦區(qū)形變具有明顯的時間特征，其一般符合logistic模型，因此可以對每一個insar時間序列采用非線性最小二乘擬合以下三參數(shù)的logistic模型：

68、

69、其中f(tins)表示tins時的insar形變量，a,b,c為待估計的模型參數(shù)。

70、求得參數(shù)之后，則對于任意時刻t可以得到其模型預計的形變：

71、

72、其中為非線性最小二乘估計得到的參數(shù)。

73、s42、此時，由s22步驟中的觀測方程可以建立模型預計值與狀態(tài)量的關系：

74、

75、一般情況下細微形變的量級都比較小，因此在添加約束方程時可以忽略小尺度的細微形變，則可以得到狀態(tài)量的約束條件：

76、gxt＝fs(t)

77、其中表示評估的空間基。

78、普通的stre并沒有包含系統(tǒng)的額外信息，但是對于礦區(qū)這樣的形變，狀態(tài)量應該滿足這些約束條件，因此利用這些額外信息來獲得比kalman濾波器更好的濾波性能。

79、根據(jù)以上所述可以得到附有狀態(tài)約束條件的stre模型，簡稱為cstre：

80、

81、其中fs(t,θ)代表建立的空間不相關，時間相關的函數(shù)約束模型，θ表示模型參數(shù)。

82、s43、最初的stre模型是基于線性系統(tǒng)構建的，然而礦區(qū)形變具有高度非線性特性，其狀態(tài)量也表現(xiàn)為高度非線性。因此原本的stre模型中存在超出系統(tǒng)之外的誤差，即完整的系統(tǒng)描述與標準時空kalman濾波所建立的描述不同。此時附加的狀態(tài)約束可以提高原本stre的性能。

83、kalman濾波是一種狀態(tài)最優(yōu)估計方法，在添加模型約束之后，應該保持這中最優(yōu)估計特性，因此需要求解最優(yōu)約束狀態(tài)。

84、針對這種附加先驗約束的kalman濾波問題，已經(jīng)有學者做了大量研究，一般是將標準kalman濾波的求解轉換為標準kalman濾波和二次規(guī)劃問題的組合，并取得了較好的性能。

85、用表示t時刻的約束狀態(tài)量，用和表示無約束kalman估計的最優(yōu)狀態(tài)量，用無約束kalman估計后驗方差，則根據(jù)已有研究，一種求解方法是將無約束的狀態(tài)估計投影到約束面上，因此根據(jù)最小方差準則建立如下目標函數(shù)：

86、

87、其中w是正定權矩陣，可以將其設置為或者單位矩陣i。

88、則該約束kalman的最優(yōu)解可以表示為：

89、

90、由于引入了具有不確定性的模型，因此還需要重新估計附有模型約束的狀態(tài)量方差，根據(jù)誤差傳播率可得：

91、

92、其中j＝wg′(gwg′)-1。

93、以上約束表明當kalman評估的狀態(tài)量方差較小時，可以認為其更可靠，此時給模型約束賦予較小的權重；當kalman評估的狀態(tài)量方差較大時，認為其需要糾正，此時給模型約束賦予較大的權重。通過這種方式，可以保證kalman濾波結果不偏離模型擬合結果。在設計時應當注意距離gnss點較近的狀態(tài)量盡量保持不變，以免破壞gnss的精度。

94、至此便得到了附加約束的狀態(tài)量及其方差的估計，當進行下一輪的kalman濾波時，只需要在狀態(tài)轉移預測中用代替用代替然后進行新的迭代即可。

95、通過以上過程可以獲得形變時空特征約束的insar和gnss融合的時空kalman模型，即cstre模型的構建過程。

96、根據(jù)本發(fā)明的一些優(yōu)選實施方式，所述s5包括：

97、s51、insar空間建模后的殘差信息中包含了小尺度的細微形變和噪聲，需要分離這兩種信息的方差。對每一景insar形變空間建模后的殘差采用半變異函數(shù)進行擬合，獲得不同距離上的方差值：

98、

99、式中，h為形變點之的間距；n(h)為以h為間距的所有觀測點的成對數(shù)目；z(xi)和z(xi+h)分別表示相對距離為h的insar形變觀測值；為間距為h的方差，即變異值，在一定范圍內(nèi)隨h的增大而增大。當實測點距離大于最大相關距離時，該值趨于穩(wěn)定。通過本步驟可以獲得insar殘余形變的方差與間距h的樣本值。

100、s52、采用經(jīng)典球狀函數(shù)模型對每個時間的insar數(shù)據(jù)樣本進行擬合，獲得塊金、變程、基臺和偏基臺值等模型參數(shù)。其中塊金體現(xiàn)的就是insar空間不相關的噪聲ω(s,t)的方差而基臺體現(xiàn)的則是空間相關的小尺度的細微形變ξt(s)的方差通過以上操作分離得到了insar空間相關的小尺度細微形變的方差和觀測噪聲的方差。

101、s53、由于insar和gnss的觀測精度不同，因此可以根據(jù)s52步驟得到的insar小尺度細微形變的方差和觀測噪聲的方差采用方差分量估計來確定這兩種觀測方法的精度比k，然后根據(jù)下式計算gnss的觀測噪聲：

102、

103、可得到gnss的觀測噪聲為：

104、

105、s54、由以上分析可知cstre模型的未知參數(shù)為：小尺度的細微形變ξt(s)、附加約束的狀態(tài)量xt、狀態(tài)轉移矩陣h和狀態(tài)轉移過程中的方差

106、在求解之前其余參數(shù)可以設定初始經(jīng)驗值，而狀態(tài)轉移矩陣h的初始值，一般狀態(tài)轉移矩陣的初始值設置為：

107、h＝ρe

108、其中0<ρ<1是乘數(shù)因子，e是單位矩陣。

109、將以上步驟求得的參數(shù)帶入融合模型中進行求解。該融合模型一般采用向前濾波和向后平滑進行求解。即依次計算：狀態(tài)量預測、狀態(tài)量先驗方差估計、計算時空kalman增益、狀態(tài)量及其方差的后驗最優(yōu)估計、約束狀態(tài)量及其方差估計。然后對所有時間重復以上過程，最后對整個融合過程用em迭代進一步提高參數(shù)估計的精度。

110、根據(jù)本發(fā)明的一些優(yōu)選實施方式，所述s6包括：

111、s61、根據(jù)上述步驟評估獲得模型參數(shù)，重建地表高時空分辨率形變：

112、

113、其中，s表示重建后空間點的位置，一般與insar點的分布一致，t表示重建點的時間，一般與gnss觀測時間一致。

114、s62、針對重建的形變監(jiān)測結果可以分別在時間和空間上進行精度驗證。

115、在時間上，利用一部分未參與融合的地表真實gnss數(shù)據(jù)，計算其與融合結果的均方根誤差。

116、在空間上，為了方便精度驗證，可以在重建前將一部分insar數(shù)據(jù)刪除，使其不參與融合，然后用融合結果和這部分insar數(shù)據(jù)檢校，計算其均方根誤差。

117、通過以上兩個過程可以分別驗證本發(fā)明提出的方法在時間和空間上形變的重建精度。

118、根據(jù)上述時空kalman模型融合insar和gnss數(shù)據(jù)重建礦區(qū)地表形變監(jiān)測方法，可以得到一種基于先驗形變特征約束的時空kalman模型融合insar和gnss數(shù)據(jù)重建礦區(qū)地表高時空分辨率的形變監(jiān)測系統(tǒng)。

119、上述融合方法或系統(tǒng)可用于礦區(qū)及其它場景的形變監(jiān)測中。

120、本發(fā)明針對形變空間變化不均勻且梯度大、時間變化高度非線性的礦區(qū)，提出的時空kalman模型融合insar和gnss數(shù)據(jù)重建礦區(qū)地表形變監(jiān)測方法充分利用了insar技術高空間分辨率和gnss技術高時間分辨率的優(yōu)勢。根據(jù)礦區(qū)形變的空間特征自適應構建空間模型，提高了算法的自動性和效率。同時附加礦區(qū)先驗形變的時間特征約束，不僅減少了傳統(tǒng)stre融合方法對gnss數(shù)據(jù)的依賴，而且進一步提高了形變?nèi)诤系木取Ｔ摲椒軌蛲瑫r顧及礦區(qū)的insar形變監(jiān)測結果、gnss形變監(jiān)測結果和先驗形變特征，有效地將這三種信息最大程度地融合，實現(xiàn)了對礦區(qū)高時空分辨率的形變監(jiān)測。