融合像元偏移跟蹤和短基線集的礦區地表形變解算方法
【專利摘要】一種融合像元偏移跟蹤和短基線集的礦區地表形變解算方法。適用于地表礦區地表形變研究使用。其步驟為：對礦區地表形變區域進行預估；選取時間序列上相干圖中的高相干點；利用時序偏移量跟蹤算法恢復大形變區域高相干點的相位整周數N；大形變區域和小形變區域的高相干點聯合建立短基線集解算模型；礦區地表高程誤差及形變速率的解算。其監測精度高、范圍大，實現操作過程簡單，費用低，克服了傳統時序InSAR方法無法正確獲取礦區大形變梯度下的地表沉降問題，也解決了基于SAR幅度信息的時序像元偏移跟蹤算法和基于差分相位的短基線集技術難以聯合解算地表形變速率和高程誤差的問題，具有廣泛的實用性。
【專利說明】
融合像元偏移跟蹤和短基線集的礦區地表形變解算方法
技術領域
[0001] 本發明設及一種礦區地表形變解算方法，尤其適用于一種融合像元偏移跟蹤和短 基線集的礦區地表形變解算方法。 技術背景
[0002] 合成孔徑雷達干設測量技術（InSAR，Interferometric Synthetic Aperture Radar)具有全天候、覆蓋范圍大、精度高、具有穿透力等優點，在地形測繪、形變監測、地表 參數反演等領域的應用越來越廣。DInSAR(DifferentialInSAR)主要用于地表形變監測， 在火山監測、地面沉降、地震形變場獲取、滑坡監測等領域具有極大的應用潛能。但該技術 受到時間基線、空間基線、大氣延遲等因素的影響，僅對部分高相干點或永久散射點進行處 理分析的時序111541?方法應運而生。該類方法主要包括:？5-1]1541?(化1'1]1日]1日]118。日1:1日的1·- InSAR)、SBAS(Small Baseline Subsets)、IPTA(Interferometric Point Target Analysis)等。
[0003] 時序InSAR處理技術中，因主影像不唯一，時空基線較短，因此相干性較高的短基 線集技術(SBAS)應用較廣。但該方法在解算時首先要對差分相位進行解纏，而礦區開采沉 陷速度快、形變量大，相位解纏只能解算相鄰點形變梯度在(-η，π]間的沉降值，大形變區域 的地表下沉難W通過相位解纏方法正確得到，導致運些區域監測量級較小，難W滿足礦區 大形變監測要求。同時，像元偏移跟蹤算法能夠解算地表大形變，目前主要應用于冰川移 動、地震形變等，但利用SAR幅度信息得到的地表形變精度較低，并且非常耗時，解算結果中 的地形誤差與基于差分相位建立的短基線集技術中的地形誤差并不等價，無法將其帶入短 基線集模型進行直接解算。因此，目前沒有將像元偏移跟蹤和短基線集融合解算地表形變 和局程誤差的方法。
【發明內容】
[0004] 針對上述技術問題，提供一種克服了時序InSAR方法缺點，計算量小，精度高的融 合像元偏移跟蹤和短基線集的礦區地表形變解算方法。
[0005] 為實現上述技術目的，本發明的融合像元偏移跟蹤和短基線集的礦區地表形變解 算方法步驟如下：
[0006] 步驟1:預估礦區地表形變區域：
[0007] 選取被測因煤炭開采導致的礦區，利用衛星雷達獲取選取礦區的SAR影像，利用 DInSAR或開采沉陷預計得到的地表沉降結果，將SAR影像區域中相鄰像元形變差大于SAR影 像波長一半的下沉盆地初步劃分成大形變區域，將像元形變量小于10mm設定為無形變區 域，其余區域設定為小形變區域；
[0008] 步驟2:選取時間序列上相干圖中的高相干點：
[0009] 設定SAR影像的時間和空間基線闊值，選取小于闊值的SAR影像進行干設處理，再 依據相干性獲取整個或大部分時間序列上具有高相干性的點位；
[0010] 步驟3:利用時序偏移量跟蹤算法恢復大形變區域高相干點的相位整周數N:
[0011] 在劃分出的無形變區域內選取相干性大于0.3的高相干點計算整體偏移量，利用 最小二乘原理建立偏移量的雙線性多項式內插函數，將形變區域內相干點的影像坐標帶入 內插函數計算形變區域內相干點的整體偏移量；
[0012] 利用小窗口最大頻譜配準對大形變區域的高相干點進行精確配準，獲取局部區域 偏移量；
[0013] 從局部精確偏移量中去除整體偏移量，并去除未下沉區域因地形起伏引起的殘 差，從而獲取最終的視線向像兀偏移量off set;
[0014] 利用視線向像元偏移量offset得到大形變區域高相干點的相位整周數N，
[0015] 利用相位整周數N和大形變區域內高相干點的原差分纏繞相位聯合恢復大形變區 域高相干點的真實相位差S Ob;
[0016] 步驟4:利用大形變區域和小形變區域的高相干點聯合建立短基線集解算模型：
[0017] 在隨機時間點ti、隨機時間點t2分別采集了兩幅被測礦區的SAR影像信息，判斷步 驟2選取的相干點位于SAR影像的小形變和無形變區域時，通過傳統最小費用流方法解纏差 分干設圖中高相干點（1，m)的真實相位差δ Φ S;
[0018] 判斷高相干點（l，m)位于大形變區域時，不能采用傳統最小費用流方法解纏差分 干設圖中高相干點（l，m)真實相位差，此時則由步驟3獲取高相干點（l，m)的真實相位差δ 巫b ;
[0019]最終，利用δΦ3、δ(1Π 聯合構建短基線集的解算模型；
[0020] 步驟5,解算礦區地表高程誤差及形變速率：
[0021] 利用短基線集處理模型，采用最小二乘解算礦區高相干點處地表高程誤差及形變 速率，在獲取礦區地表形變后，通過反演礦區概率積分法預計參數，最終獲取礦區地表沉陷 預計及分析模型，為礦區沉陷工程應用提供基礎數據。
[002^ 所述未下沉區域因地形起伏引起的殘差利用公式:峨妊,。,。=。。巧+。2(//-豆V 計算得到；式中，offsettopo為地形起伏引起的偏移量誤差；曰日，曰1，曰2為最小二乘的擬合系 數;Η為對應位置的高程;忌為研究區域的平均高程；
[0023] 利用公式:N = ceil(2offset · Δ/λ)，計算得到大形變區域高相干點的相位整周 數Ν，式中，ceil()為向下取整運算，Δ為斜距向像元尺寸，λ為SAR影像波長;從而克服像元 偏移跟蹤算法求解的偏移量僅用到SAR影像幅度信息，不能同傳統短基線技術得到的差分 相位進行最小二乘建模的問題。
[0024] 利用公式：S〇b = 2炯+δφ<Η??，計算大形變區域高相干點的真實相位差，式中，δ Φ diff為纏繞的差分相位；
[0025] 利用真實相位差δ Ob建立短基線集處理模型A方法為:首先利用真實相位差δ Ob構建解 纏相位矩陣δ Φ = [ δ Φ sS Φ b]，再利用公式：
建立短基線集處理模型;式中：A hi,m為相干點（l，m)的高程誤差;V為相干點的下沉速率;W 為殘余相位;λ為衛星雷達波長;B丄,i,m為兩幅SAR影像相干點處的垂直基線；0(l，m)為相干 點處的雷達入射角;ri,m為衛星距離點目標(1，m)的斜距；
[0026] 局部區域偏移量還可W通過最小二乘匹配方法對大形變區域的高相干點進行精 確配準獲取；
[0027] 所述礦區高相干點處地表高程誤差及形變速率利用公式：
計算，式中P為權陣，可取單位陣。
[0028] 有益效果:本發明將采用幅度獲取地表大形變的時序偏移量跟蹤算法和利用相位 獲取地表微小形變的時序短基線算法有機融合，將兩者納入到一個模型中統一計算，利用 兩種方法的互補性，通過利用真實相位差建立短基線集處理模型解決了兩者無法統一建模 計算和傳統時序InSAR方法無法正確獲取礦區大形變梯度下的地表沉降問題，改變了原有 短基線集算法中只能解算地表微小形變的方式，現有時序偏移量方法不能有效解算地形誤 差的問題也得到了的解決。方法簡單，監測精度高、范圍大，費用低，對指導礦區生產、預警 地質災害、治理生態環境、維護礦群關系等具有重要的實際意義和應用價值。
【附圖說明】
[0029] 圖1為本發明實施例一種融合像元偏移跟蹤和短基線集的礦區地表形變解算方法 流程圖
[0030] 圖2為本發明實施例采用相干性闊值選取的高相干點分布圖 [0031 ]圖3為本發明實施例解算的相干點處地表年形變速率圖
[0032] 圖4為本發明實施例估計的相干點處DEM誤差分布圖
【具體實施方式】
[0033] 下面將結合圖和具體實施過程對本發明做進一步詳細說明：
[0034] 如圖1所示，本發明的融合像元偏移跟蹤和短基線集的礦區地表形變解算方法，包 括如下步驟：
[0035] 步驟1:礦區地表形變區域預估：
[0036] 本實例選取被測因煤炭開采導致的礦區，采用的是裁切后的利用衛星獲取選取礦 區的11景Im分辨率的Terrasar-X衛星影像，范圍為2000 X 1600像元。利用DInSAR或開采沉 陷預計得到的地表沉降結果，將SAR影像區域中相鄰像元形變差大于SAR影像波長一半的下 沉盆地初步劃分成大形變區域，將像元形變量小于10mm設定為無形變區域，其余區域設定 為小形變區域；
[0037] 步驟2:選取時間序列上相干圖中的高相干點：
[003引設定SAR影像的時間和空間基線闊值，選取小于闊值的SAR影像進行干設處理，再 依據相干性獲取整個或大部分時間序列上具有高相干性的點位;設定時間和空間基線闊值 分別為11天和160m，選取小于闊值的SAR影像進行干設處理，W相干性闊值為0.3選取時間 序列上的高相干性點位，本發明實施例采用相干性闊值選取的高相干點分布圖如圖2所示， 共選出885997個點，中部大形變區域內因失相干影響，只選出了部分點位；
[0039] 步驟3:利用時序偏移量跟蹤算法恢復大形變區域高相干點的相位整周數N:
[0040] 在劃分出的無形變區域內選取相干性大于0.3的高相干點計算整體偏移量，基于 最小二乘原理建立偏移量的雙線性多項式內插函數，將形變區域內相干點的影像坐標帶入 內插函數計算形變區域內相干點的整體偏移量；
[0041] 利用小窗口最大頻譜配準對大形變區域的高相干點進行精確配準，獲取局部區域偏移 量;從局部精確偏移量中去除整體偏移量，并去除利用公式:<斯'啤I aW W I 3:W //)2 得到的未下沉區域因地形起伏引起的殘差，從而獲取最終的視線向像元偏移量offset,式 中，offsettopo為地形起伏引起的偏移量誤差;曰日，曰1，恥為最小二乘的擬合系數;Η為對應位 置的高程;77為研究區域的平均高程；
[0042] 由于利用上述像元偏移跟蹤算法求解的偏移量僅用到SAR影像幅度信息，不能同 傳統短基線技術得到的差分相位進行最小二乘建模。利用公式:N = ceil(2offset . Δ/λ) 得到大形變區域高相干點的相位整周數Ν，式中，ceil()為向下取整運算，Δ為斜距向像元 尺寸，λ為SAR影像波長。由相位整周數Ν和大形變區域內高相干點的原差分纏繞相位聯合恢 復大運些點的真實相位差δΦ6 = 2侃+Sd)diff，式中，Sd)diff為纏繞的差分相位；
[0043] 步驟4:大形變區域和小形變區域的高相干點聯合建立短基線集解算模型：
[0044] 在隨機兩個時間點ti，t2時刻分別獲取了兩幅被測礦區的SAR影像，判讀相干點位 于SAR影像的小形變和無形變區域時，通過傳統最小費用流方法解纏差分干設圖中高相干 點（l，m)的相位δΦ3相位;判讀相干點位于大形變區域時，相位解纏結果則不可靠，需利用 步驟3得到的5〇b構建解纏相位矩陣：δΦ = [S0sS0b]，然后直接利用公式：
i立模型，改變了原有短基線集算 法中只能解算地表微小形變的方式；
[0045] 式中：Ahi,m為相干點（l，m)的高程誤差;V為相干點的下沉速率;W為殘余相位;λ為 雷達波長;為兩幅SAR影像相干點處的垂直基線；0(l，m)為相干點處的雷達入射角； ri,m為衛星距離點目標（1，m)的斜距；
[0046] 步驟5,礦區地表高程誤差及形變速率的解算
[0047] 利用步驟4中的短基線集處理模型，采用最小二乘解算礦區高相干點處地表高程 誤差及形變速率
其中P為權陣，可取單位陣。
[004引本發明實施例解算的相干點處地表年形變速率如圖3所示，從中可W看出，每年下 沉大于1000mm的部分地表點已經可W得到，解決了短基線集技術無法監測地表大變形的問 題。圖4為本發明實施例估計的相干點處DEM誤差分布圖。因此，本發明解決了基于SAR幅度 信息的時序像元偏移跟蹤算法和基于差分相位的短基線集技術難W聯合解算地表形變速 率和高程誤差的問題，為礦區地表沉降時序監測提供了新的處理方法。
【主權項】
1. 一種融合像元偏移跟蹤和短基線集的礦區地表形變解算方法，其特征在于步驟如 下： 步驟1:預估礦區地表形變區域： 選取被測因煤炭開采導致的礦區，利用衛星雷達獲取選取礦區的SAR影像，利用DInSAR 或開采沉陷預計得到的地表沉降結果，將SAR影像區域中相鄰像元形變差大于SAR影像波長 一半的下沉盆地初步劃分成大形變區域，將像元形變量小于l〇mm設定為無形變區域，其余 區域設定為小形變區域； 步驟2:選取時間序列上相干圖中的高相干點： 設定SAR影像的時間和空間基線閾值，選取小于閾值的SAR影像進行干涉處理，再依據 相干性獲取整個或大部分時間序列上具有高相干性的點位； 步驟3:利用時序偏移量跟蹤算法恢復大形變區域高相干點的相位整周數N: 在劃分出的無形變區域內選取相干性大于0.3的高相干點計算整體偏移量，利用最小 二乘原理建立偏移量的雙線性多項式內插函數，將形變區域內相干點的影像坐標帶入內插 函數計算形變區域內相干點的整體偏移量； 利用小窗口最大頻譜配準對大形變區域的高相干點進行精確配準，獲取局部區域偏移 量； 從局部精確偏移量中去除整體偏移量，并去除未下沉區域因地形起伏引起的殘差，從 而獲取最終的視線向像元偏移量offset; 利用視線向像元偏移量off set得到大形變區域高相干點的相位整周數N， 利用相位整周數N和大形變區域內高相干點的原差分纏繞相位聯合恢復大形變區域高 相干點的真實相位差 步驟4:利用大形變區域和小形變區域的高相干點聯合建立短基線集解算模型： 在隨機時間點t、隨機時間點t2分別采集了兩幅被測礦區的SAR影像信息，判斷步驟2選 取的相干點位于SAR影像的小形變和無形變區域時，通過傳統最小費用流方法解纏差分干 涉圖中高相干點（1，m)的真實相位差δ Φ s; 判斷高相干點（l，m)位于大形變區域時，不能采用傳統最小費用流方法解纏差分干涉 圖中高相干點（1，m)真實相位差，此時則由步驟3獲取高相干點（1，m)的真實相位差δ Φ b; 最終，利用S Φ s、S Φ b聯合構建短基線集的解算模型； 步驟5,解算礦區地表高程誤差及形變速率： 利用短基線集處理模型，采用最小二乘解算礦區高相干點處地表高程誤差及形變速 率，在獲取礦區地表形變后，通過反演礦區概率積分法預計參數，最終獲取礦區地表沉陷預 計及分析模型，為礦區沉陷工程應用提供基礎數據。2. 根據權利要求1所述的融合像元偏移跟蹤和短基線集的礦區地表形變解算方法，其特征 在于:所述未下沉區域因地形起伏引起的殘差利用公式:喻… 計算得到；式中，off sett%。為地形起伏引起的偏移量誤差;a〇，ai，a2為最小二乘的擬合系 數;Η為對應位置的高程;互為研究區域的平均高程。3. 根據權利要求1所述的融合像元偏移跟蹤和短基線集的礦區地表形變解算方法，其 特征在于:利用公式:N=ceil(2offset · △ /λ)，計算得到大形變區域高相干點的相位整周 數N，式中，ceil〇為向下取整運算，Λ為斜距向像元尺寸，λ為SAR影像波長;從而克服像元 偏移跟蹤算法求解的偏移量僅用到SAR影像幅度信息，不能同傳統短基線技術得到的差分 相位進行最小二乘建模的問題。4. 根據權利要求1所述的融合像元偏移跟蹤和短基線集的礦區地表形變解算方法，其 特征在于：利用公式：δ Φb = 2Νπ+δ φ diff，計算大形變區域高相干點的真實相位差，式中，δ Φ diff為纏繞的差分相位。5. 根據權利要求1所述的融合像元偏移跟蹤和短基線集的礦區地表形變解算方法，其特征在 于:利用真實相位差S Φ b建立短基線集處理模型A方法為:首先利用真實相位差δ Φ b構建解纏相 位矩陣δΦ = [δΦ3 δΦ」，再利用公式：建立短基線集處理模型;式中：Ain,?為相干點（l，m)的高程誤差;V為相干點的下沉速率;W 為殘余相位;λ為衛星雷達波長;為兩幅SAR影像相干點處的垂直基線；0(l，m)為相干 點處的雷達入射角;n, m為衛星距離點目標(1，m)的斜距。6. 根據權利要求1所述的融合像元偏移跟蹤和短基線集的礦區地表形變解算方法，其 特征在于:局部區域偏移量還可以通過最小二乘匹配方法對大形變區域的高相干點進行精 確配準獲取。7. 根據權利要求1或6所述的融合像元偏移跟蹤和短基線集的礦區地表形變解算方法，其特 - Λ - 征在于:所述礦區高相干點處地表高程誤差及形變速率利用公式： _ V _ ; 計算，式中Ρ為權陣，可取單位陣。
【文檔編號】G01S13/90GK106066478SQ201610362009
【公開日】2016年11月2日
【申請日】2016年5月27日 公開號201610362009.6, CN 106066478 A, CN 106066478A, CN 201610362009, CN-A-106066478, CN106066478 A, CN106066478A, CN201610362009, CN201610362009.6
【發明人】范洪冬, 杜森, 黃繼磊, 閆世勇, 鄧喀中, 汪云甲
【申請人】中國礦業大學