結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法
【專利摘要】本發(fā)明公開了一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，充分融合了InSAR監(jiān)測手段分辨率高、相對精度高、覆蓋范圍大和地面控制點(diǎn)絕對精度高的特點(diǎn)，既解決InSAR獲取DEM絕對精度不高，又克服了傳統(tǒng)DEM獲取方法例如星載和機(jī)載光學(xué)及激光測高手段獲取高分辨測圖效費(fèi)比低和受云層影響，地面測圖費(fèi)時(shí)費(fèi)力且無法得到整個(gè)野外山區(qū)高分辨地貌等缺點(diǎn)。整個(gè)流程結(jié)構(gòu)清晰，具有實(shí)現(xiàn)簡單、費(fèi)用低、監(jiān)測精度高、監(jiān)測范圍大、自動(dòng)化程度高等優(yōu)點(diǎn)。
【專利說明】
結(jié)合外控點(diǎn)的I nSAR高精度高分辨率DEM獲取方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明屬于基于遙感影像的大地測量領(lǐng)域，是一種利用InSAR(合成孔徑干涉雷 達(dá))技術(shù)和外部控制點(diǎn)獲取高分辨率高精度DEM的方法。
【背景技術(shù)】
[0002] 由于具備全天候，大范圍，高精度的觀測能力，合成孔徑干涉雷達(dá)測量 (Interferometric SyntheticAperture radar，InSAR)已經(jīng)成為國內(nèi)外熱點(diǎn)關(guān)注技術(shù)。隨 著各種 SAR 傳感器的升空，例如 JERS、ERS-l/2、ENVISAT、AL0S/PALSAR、AL0S-2/PALSAR-2、 1^(^"3丨-1/2、161瓜541?1、0)5]\?)-51^]^(1、5611乜1161-1等，可用541?數(shù)據(jù)的時(shí)間和空間分辨 率不斷提高，該技術(shù)在監(jiān)測由于地下水抽取、煤礦開采、石油開采、地下鐵路修建、填海、地 震、冰川等引起的地表形變方面得到了空前的發(fā)展和應(yīng)用。InSAR監(jiān)測地表形變主要是借助 于覆蓋同一個(gè)地區(qū)的兩幅或兩幅以上SAR圖像，利用包含在SAR圖像中的差分相位信息獲取 來實(shí)現(xiàn)。令A(yù) (i>de5f為地表形變引起的相位變化，A為雷達(dá)波長，則形變值A(chǔ) D可通過下面的模 型得到：
(1)
[0004]當(dāng)a (})def = 231時(shí)，對應(yīng)的形變大小，即形變模糊度為：
(2)
[0006]除監(jiān)測形變以外，InSAR技術(shù)另一個(gè)基本應(yīng)用是通過在同一地區(qū)不同角度拍攝的 SAR圖像的相位差別來提取地形起伏。令A(yù) (}> top為地表起伏A h引起的相位變化，R為雷達(dá)視 距，#為雷達(dá)入射角，為垂直基線，則相位差到地形起伏的轉(zhuǎn)換模型為：
⑶
[0008]干涉相位圖上231相位變化引起的高程變化，即模糊高為
[0010] A-般為幾厘米級到幾分米，B±-般為幾十米到上千米，而R為幾百千米，顯然，地 形起伏對相位變化的敏感度要遠(yuǎn)大于形變對相位變化的敏感度。常規(guī)的重復(fù)軌道觀測相位 變化往往包含了兩次觀測間大氣變化引起的相關(guān)部分和去相干引起的噪聲。這些相位誤差 會(huì)極大地影響InSAR高程提取結(jié)果。因此，除美國航空局NASA的SRTM項(xiàng)目（以獲取DEM為目標(biāo) 的雙天線單次觀測）以外，其他單天線重復(fù)觀測SAR項(xiàng)目長期以來以形變監(jiān)測應(yīng)用為主。DEM 即為數(shù)字高程模型，鑒于全球多個(gè)領(lǐng)域?qū)Ω呔雀叻直媛蔇EM的需求，德國宇航局DLR開發(fā) 了具有創(chuàng)新性的一發(fā)雙收觀測模式，即雙基(Bistatic)模式。由串聯(lián)飛行的TerraSAR-X或 TanDEM-X衛(wèi)星發(fā)射雷達(dá)信號(hào)，而TerraSAR-X與TanDEM-X衛(wèi)星同時(shí)接收回波信號(hào)。兩個(gè)衛(wèi)星 接收回波信號(hào)的時(shí)間差極短(一分鐘內(nèi)），這樣保證了兩次成像的大氣狀態(tài)和地表散射特性 的高度相似性，進(jìn)而保證的相位觀測的精度。雙基模式下相位差到地形起伏的轉(zhuǎn)換模型為：
〇) L〇〇12」 目前TerraSAR-X和TanDEM-X獲取的影像對（即TanDEM-X CoSSC影像對)像元大小 可以達(dá)到3m以內(nèi)，通過迭代差分干涉的方式能滿足高分辨DEM生產(chǎn)像元大小這一塊的要求。 但是通過差分干涉方法獲取的DEM是相對于外部DEM的高程值，其絕對精度與解纏參考點(diǎn)的 選取有很大關(guān)系，需要外部控制點(diǎn)來進(jìn)行絕對高程校正;而且，通過InSAR手段獲取的初始 DEM，其平面位置精度主要由地理編碼決定，具體來說是通過真實(shí)SAR影像與外部DEM模擬的 SAR影像之間的強(qiáng)度匹配來實(shí)現(xiàn)(配準(zhǔn)法則是窗口內(nèi)強(qiáng)度相關(guān)系數(shù)達(dá)到最大）。匹配的結(jié)果 一般是通過擬合得到的表征系統(tǒng)偏移量的二次多項(xiàng)式來檢核，也就是基于內(nèi)符合精度。盡 管目前SAR影像星歷參數(shù)和軌道數(shù)據(jù)精度均較高，但山區(qū)坡度起伏大的特點(diǎn)決定即使細(xì)小 的平面偏差也會(huì)引起較大的高程誤差。因此需要一種匹配法則不同于強(qiáng)度相關(guān)最大的配準(zhǔn) 方法來檢核和校正初始DEM的平面和垂直精度。此外，目前通過軌道來估計(jì)基線一般都包含 殘差?；€殘差影響傳遞到DEM產(chǎn)品表現(xiàn)為平面上的系統(tǒng)誤差趨勢；SAR的斜距成像使得影 像不可避免地出現(xiàn)透視收縮。透視收縮的強(qiáng)度與局部坡度和坡向密切相關(guān)。因此還要考慮 初始DEM產(chǎn)品與坡度和坡向相關(guān)的系統(tǒng)誤差。上述系統(tǒng)誤差都要借助外部控制點(diǎn)來進(jìn)行最 小二乘平差來削弱或者消除。
[0013] 綜上所述，目前采用TanDEM-X CoSSC數(shù)據(jù)和InSAR差分干涉手段可以得到高分辨 的DEM產(chǎn)品，但其高程是相對的，而且受各種現(xiàn)有的InSAR干涉處理過程中無法完全去除的 誤差影響，需要外部控制點(diǎn)來進(jìn)行校正以提高絕對的平面和垂直精度。

【發(fā)明內(nèi)容】
：
[0014] 為了克服上述【背景技術(shù)】的缺陷，本發(fā)明提供一種結(jié)合外控點(diǎn)的InSAR高精度高分 辨率DEM獲取方法，具有實(shí)現(xiàn)簡單、監(jiān)測精度高范圍大、自動(dòng)化程度高的優(yōu)點(diǎn)。
[0015] 為了解決上述技術(shù)問題本發(fā)明的所采用的技術(shù)方案為：
[0016] 一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，包括：
[0017] 步驟1，通過雙基InSAR差分干涉處理單軌影像對，獲取單軌DEM;
[0018]步驟2,采取不同軌道的影像對分別獲取DEM，并進(jìn)行逐點(diǎn)融合；
[0019]步驟3,將經(jīng)融合的DEM與地面控制點(diǎn)進(jìn)行配準(zhǔn)；
[0020]步驟4,校正DEM與坡度、平面位置和坡向相關(guān)的系統(tǒng)誤差。
[0021]較佳地，步驟1包括：
[0022]步驟1.1，根據(jù)DEM需求區(qū)域坐標(biāo)范圍定制影像對，影像對包括主影像和輔影像；
[0023] 步驟1.2，將主影像和輔影像共輒相乘生成初始干涉圖；
[0024] 步驟1.3通過軌道數(shù)據(jù)估計(jì)基線，根據(jù)基線由外部DEM模擬包含平地相位的地形相 位圖；
[0025] 步驟1.4將地形相位圖與初始干涉相位圖相減得到差分干涉圖；
[0026] 步驟1.5對差分干涉圖進(jìn)行相位濾波和相位解纏，將解纏相位轉(zhuǎn)化為高程差；
[0027] 步驟1.6將高程差與外部DEM疊加得到新DEM;
[0028] 步驟1.7通過地理編碼將新DEM從SAR坐標(biāo)系轉(zhuǎn)換到外部DEM的地理坐標(biāo)系得到單 軌 DEM。
[0029] 較佳地，還包括:將步驟1.7所得單軌DEM作為外部DEM，對步驟1.1至步驟1.7進(jìn)行 迭代，直至獲取原始影像像元大小的DEM作為單軌DEM。
[0030]較佳地，步驟1.5對差分干涉圖進(jìn)行相位濾波采用的濾波方式包括Nonlocal均值 濾波方式，濾波模型包括：
[0032] s為濾波窗口內(nèi)的中心像元，t為濾波窗口內(nèi)的其他任意像元，真實(shí)的干涉相位為 U，含噪聲的原始干涉相位為巾，相位估計(jì)值U(S)為濾波結(jié)果，W(S，t)為由像元t的鄰域像元 與像元S的鄰域像元之間的SAR影像強(qiáng)度和相干性的綜合相似度而決定的權(quán)重。
[0033]較佳地，步驟2的具體步驟包括：
[0034] 步驟21，將不同軌道的影像對統(tǒng)一配準(zhǔn)到相同坐標(biāo)系之下，不同軌道的影像對包 括升軌和降軌的影像對；
[0035] 步驟22，采用融合模型進(jìn)行融合，融合模型包括 其中，Hfuslcin為融合后高程，83、13和仏分別為升軌的垂直基線值，疊影或陰影掩膜值及高程 值，Bd、Wd和Hd分別為降軌的垂直基線值，疊影或陰影掩膜值及高程值。
[0036]較佳地，步驟3包括：
[0037]步驟3.1，將外部控制點(diǎn)的空間直角坐標(biāo)轉(zhuǎn)為大地經(jīng)煒度坐標(biāo)；
[0038]步驟3.2,采用雙線性插值提取外部控制點(diǎn)所在DEM高程值；
[0039]步驟3.3,將經(jīng)融合的DEM與數(shù)個(gè)外部控制點(diǎn)進(jìn)行統(tǒng)計(jì)配準(zhǔn)，配準(zhǔn)模型包括
[0040] 其中，0為控制點(diǎn)位置的局部坡向，a為控制點(diǎn)位置的局部坡度，5為所有控制點(diǎn)位 置坡度的平均值，a為DEM平面偏移矢量的模，b為DEM偏移矢量的方向角，A Z為DEM垂直偏 移，A h為DEM任一點(diǎn)的高程值與控制點(diǎn)高程值的差。
[0041] 較佳地，對配準(zhǔn)過程進(jìn)行迭代，直至高程差標(biāo)準(zhǔn)差改善率低于2 %。
[0042] 較佳地，步驟4對DEM結(jié)果與平面位置相關(guān)的誤差進(jìn)行校正，是通過二次曲面擬合 來完成的，二次曲面擬合模型包括
[0043] A h(x,y) =aix2+a2y2+a3xy+a4x+a5y+a6
[0044] 其中x和y表示某個(gè)像元的橫坐標(biāo)和縱坐標(biāo)，A h(x,y)為該點(diǎn)的DEM值與控制點(diǎn)高程 差，ai……a 6為待估模型參數(shù)。采用解算控制點(diǎn)位的位置與高程差作為觀測值，通過最小二 乘平差方法來求解待估參數(shù)ai……a 6,然后將模型求取的誤差趨勢值與經(jīng)步驟3改正后DEM 進(jìn)行相減，得到進(jìn)一步改正后的DEM。
[0045]較佳地，步驟4對DEM結(jié)果與坡度相關(guān)系統(tǒng)誤差進(jìn)行校正是通過二次或三次多項(xiàng)式 模型來完成的，
[0046]二次或三次多項(xiàng)式包括&]1(31。1)(5) = 131(12+匕2(1+匕3或&]1(31。1)(5) = 131(13+匕2(12+匕3(1+匕4，其 中a表示某個(gè)像元點(diǎn)的坡度值，A h(sic^為該點(diǎn)的DEM值與控制點(diǎn)高程差h……b4為待估模 型參數(shù)，
[0047]采用解算控制點(diǎn)位的坡度值與高程差作為觀測值，通過最小二乘平差方法來求解 待估參數(shù)匕……b4,然后將模型求取的誤差趨勢與經(jīng)過平面位置誤差改正后的DEM進(jìn)行相 減，得到進(jìn)一步改正后的DEM。
[0048]較佳地，步驟4對DEM結(jié)果與坡向相關(guān)系統(tǒng)誤差進(jìn)行校正是通過正弦模型來完成 的，正弦模型為 Ahaspect=Aisin( 〇0+巾）+Ci
[0049] 其中0為某個(gè)像元點(diǎn)的坡向值，A haspec；t為該點(diǎn)的DEM值與控制點(diǎn)高程差，Ai、《、小 和&為待估模型參數(shù)。
[0050] 采用解算控制點(diǎn)位的坡向值與高程差作為觀測值，通過最小二乘平差方法來求解 待估參數(shù)Ai、co、巾和心，然后將模型求取的誤差趨勢與經(jīng)過坡度誤差改正后的DEM進(jìn)行相 減，得到進(jìn)一步改正后的DEM。
[0051] 本發(fā)明的有益效果在于:本發(fā)明所述方法在差分干涉InSAR獲取大范圍高分辨率 DEM的基礎(chǔ)上，針對差分干涉獲取的高程是相對值，SAR基線估計(jì)和地理編碼包含誤差，SAR 成像受透視收縮影響等主要問題，采用外部控制點(diǎn)來進(jìn)行約束和改善，提高InSAR差分干涉 技術(shù)獲取的高分辨率DEM的精度。充分融合了 InSAR監(jiān)測手段分辨率高、相對精度高、覆蓋范 圍大和地面控制點(diǎn)絕對精度高的特點(diǎn)，既解決InSAR獲取DEM絕對精度不高，又克服了傳統(tǒng) DEM獲取方法例如星載和機(jī)載光學(xué)及激光測高手段獲取高分辨測圖效費(fèi)比低和受云層影 響，地面測圖費(fèi)時(shí)費(fèi)力且無法得到整個(gè)野外山區(qū)高分辨地貌等缺點(diǎn)。整個(gè)流程結(jié)構(gòu)清晰，具 有實(shí)現(xiàn)簡單、費(fèi)用低、監(jiān)測精度高、監(jiān)測范圍大、自動(dòng)化程度高等優(yōu)點(diǎn)。
【附圖說明】
[0052]圖1是本發(fā)明實(shí)施例的方法流程圖。
[0053]圖2是本發(fā)明實(shí)施例統(tǒng)計(jì)配準(zhǔn)時(shí)平面偏移矢量大小在坡向?yàn)?的山坡處投影Sxy的 示意圖；
[0054]圖3是本發(fā)明實(shí)施例統(tǒng)計(jì)配準(zhǔn)時(shí)DEM某點(diǎn)位的高程與真實(shí)高程之間之差A(yù) h的示意 圖；
[0055]圖4是本發(fā)明實(shí)施例獲取的DEM與地面水準(zhǔn)檢查點(diǎn)的對比結(jié)果圖。
【具體實(shí)施方式】
[0056] 下面結(jié)合附圖和實(shí)施例對本發(fā)明做進(jìn)一步的說明。
[0057] -種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，所述外控點(diǎn)即為外部控制 點(diǎn)，如圖1所示，包括：
[0058] 步驟1，通過雙基InSAR差分干涉處理單軌影像對，獲取單軌DEM;
[0059]步驟1.1，根據(jù)DEM需求區(qū)域坐標(biāo)范圍定制影像對，影像對包括主影像和輔影像； [0060]步驟1.2，將主影像和輔影像共輒相乘生成初始干涉圖；
[0061 ]步驟1.3，通過軌道數(shù)據(jù)估計(jì)基線，根據(jù)基線由外部DEM模擬包含平地相位的地形 相位圖；
[0062] 步驟1.4,將地形相位圖與初始干涉相位圖相減得到差分干涉圖；
[0063] 步驟1.5,對差分干涉圖進(jìn)行相位濾波和相位解纏，將解纏相位轉(zhuǎn)化為高程差；
[0064]對差分干涉圖進(jìn)行相位濾波采用的濾波方式包括Nonlocal均值濾波方式，所述的 Nonlocal均值濾波方式即為非局部均值濾波方式，濾波模型包括：
(6)
[0066] s為濾波窗口內(nèi)的中心像元，t為濾波窗口內(nèi)的其他任意像元，真實(shí)的干涉相位為 U，含噪聲的原始干涉相位為巾，相位估計(jì)值U(S)為濾波結(jié)果，W(S，t)為由像元t的鄰域像元 與像元S的鄰域像元之間的SAR影像強(qiáng)度和相干性的綜合相似度而決定的權(quán)重。定義I z I、 z' |分別為兩幅SAR影像的強(qiáng)度。8氺、〖?xì)旆謩e表示兩個(gè)以像元8和〖為中心的小窗口中處于 相同位置的第k個(gè)像元。令P(k)為像元s，k與t，k的相似度。而P(k)的計(jì)算公式為：
(7)
[0068] A=(|zs,k|2+|zt,k|2+|z，s,k| 2+|z，t,k|2) (8)
[0069] B = 4(|zs，k|2|z's，k|2+|zt，k| 2|z't，k| 2+2 |zs，k| |z's，k| |zt，k| |z't，k|cos(<i)s，k- 4>t,k)) (9)
[0070] C= | zs，k | | z's，k | | zt，k | | z't，k | (10)
[0071] Vf(.v,〇 = PJ/;(^)! ；， (11) k
[0072] 其中h為平滑參數(shù)，h越大權(quán)值的差異性越小，濾波的強(qiáng)度則越大。
[0073] 步驟1.6,將高程差與外部DEM疊加得到新DEM;
[0074] 步驟1.7，通過地理編碼將新DEM從SAR坐標(biāo)系轉(zhuǎn)換到外部DEM的地理坐標(biāo)系得到單 軌 DEM。
[0075] 步驟1.8,將步驟1.7所得單軌DEM作為外部DEM，對步驟1.1至步驟1.7進(jìn)行迭代，直 至獲取原始影像像元大小的DEM作為單軌DEM。
[0076] 一次標(biāo)準(zhǔn)的差分干涉InSAR提取DEM的過程即完成。因?yàn)門anDEM-X影像初始分辨率 (3m左右）與外部DEM(例如30m分辨率的SRTM數(shù)據(jù))相差較大，在利用外部DEM進(jìn)行地理編碼 和相位模擬時(shí)，如果采用較大過采樣倍數(shù)會(huì)引入較大誤差，所以初次生成DEM時(shí)需對 TanDEM-X CoSSC影像多視以縮小外部DEM的過采樣倍數(shù)。通過迭代的方式來逐級提高DEM產(chǎn) 品的分辨率，即將上一次生成的新DEM作為外部DEM重復(fù)上述InSAR差分干涉過程，直到獲取 TanDEM-X CoSSC單視下的DEM產(chǎn)品。
[0077]在上述過程中，相位濾波可以極大地影響相位圖的質(zhì)量，進(jìn)而影響生成DEM的精 度，因此至關(guān)重要。本套方法采用的是考慮濾波窗口內(nèi)像元異質(zhì)性的Non-Local均值濾波， 該方法可以有效抑制相位噪聲而不損害相位的原本真實(shí)信息。Non-Local均值濾波方法的 核心是對濾波窗口內(nèi)的像素值進(jìn)行加權(quán)平均后賦予中心像元點(diǎn)。給予與中心像元相似度高 的像元更高的權(quán)值，反之與中心像元差異性越大更低的權(quán)值。因此權(quán)值的計(jì)算是否合理直 接關(guān)系到該估計(jì)方法的可靠性和精度。
[0078]步驟2,采取不同軌道的影像對分別獲取DEM，并進(jìn)行逐點(diǎn)融合以降低SAR影像透視 收縮和陰影的影響。
[0079] 步驟21，將不同軌道的影像對統(tǒng)一配準(zhǔn)到相同坐標(biāo)系之下，不同軌道的影像對主 要是時(shí)間間隔較短的升軌和降軌影像對；
[0080] 步驟22,采用融合模型進(jìn)行融合，融合模型包括
(12)
[0082]其中，Hfuslcin為融合后高程兒1和仏分別為升軌的垂直基線值，疊影或陰影掩膜 值(疊影或陰影為真則該值為〇,否則為1)及高程值，Bd、Wd和Hd分別為降軌的垂直基線值，疊 影或陰影掩膜值及高程值。
[0083]步驟3，將經(jīng)融合的DEM與地面控制點(diǎn)進(jìn)行配準(zhǔn)；
[0084] 步驟3.1，將外部控制點(diǎn)的空間直角坐標(biāo)轉(zhuǎn)為大地經(jīng)煒度坐標(biāo)；
[0085] 國內(nèi)常見的控制點(diǎn)平面基準(zhǔn)為北京54平面高斯投影坐標(biāo)系和1985黃海高程基準(zhǔn)， 而InSAR生成的DEM為WGS84大地經(jīng)煒度坐標(biāo)系和大地高(WGS84橢球高）。
[0086]首先將外部控制點(diǎn)的北京54平面高斯投影坐標(biāo)轉(zhuǎn)換到北京54大地經(jīng)煒度坐標(biāo)，即 (X，y，h)beijing544(B，L，H)beijing54。轉(zhuǎn)換過程中需要已知高斯投影參數(shù)如中央子午線。
[0087]其次將外部控制點(diǎn)的北京54大地經(jīng)煒度坐標(biāo)轉(zhuǎn)換到北京54空間直角坐標(biāo)，即（B， L , H) bei jing54~^ ( X , Y , Z ) bei jing54 〇
[0088]然后再將外部控制點(diǎn)的北京54空間直角坐標(biāo)轉(zhuǎn)換到WGS84空間直角坐標(biāo)，即(X，Y， Z)beijing54-(X，Y，Z)WGS84 具體轉(zhuǎn)換過程如下：
[0089]采用布爾莎七參數(shù)模型將外部控制點(diǎn)的北京54空間直角坐標(biāo)轉(zhuǎn)換到WGS84空間直 角坐標(biāo)，公式如下： 'X] 「〇: pr] 「瑪-
[0090] Y + .F + 0 +￡x ? ￥ + AF0 (13) ^warn - ，啤54. 54
[0091]其中AXo, A Yo, AZo是平移參數(shù);m是尺度比參數(shù)，ex,ey,ez是旋轉(zhuǎn)參數(shù)。
[0092]對其進(jìn)行變換得到如下式子： 爲(wèi)- 輯 -尤] m [i00 xbcipn^ 0： az0
[0093] r = ￥ +010 iL""c54 ， +A ? W (14) -Z_irGSS4. _〇 0 1 〇 _ & 'ey. 5 -
[0094]解算這七參數(shù)，至少需要三個(gè)已知的公共控制點(diǎn)的北京54空間直角坐標(biāo)（X，Y， Z)beijing54和WGS-84空間直角坐標(biāo)(X，Y，Z)WGS84,采用間接平差模型最小二乘準(zhǔn)則進(jìn)行解算， 用矩陣形式表示為：
[0095] V = Bx-l (15) 1 0 0 XbeijwgS4 0 +ez -er
[0096] 5=0 1 0 Y heijhig54 -sz 0 +￡?, (16) _〇 0 1 A -A 〇 _ 「x-
[0097] 1= f - y (17) 7 7 _ 」沢巧84 L 」k"",,g54
[0098] 其中V為觀測值改正數(shù)向量，《為待估參數(shù)殘差值，采用最小二乘迭代方法進(jìn)行參 數(shù)殘差求解，
[0099] x^(BrB) 'Br! (18)
[0100] 解得七參數(shù)，除了本身解算時(shí)的精度評定，還可以采用多余的已知控制點(diǎn)對七參 數(shù)準(zhǔn)確性進(jìn)行檢核。根據(jù)解得的七參數(shù)，就可以求出每個(gè)外部控制點(diǎn)的WGS84空間直角坐 標(biāo)。
[0101] 最后外部控制點(diǎn)的WGS84空間直角坐標(biāo)轉(zhuǎn)為WGS84大地經(jīng)煒度坐標(biāo)，即（X，Y，Z)WGS84 4(B，L，H)WGS84〇
[0102] 步驟3.2,采用雙線性插值提取外部控制點(diǎn)所在DEM高程值；
[0103] 外部控制點(diǎn)在DEM格網(wǎng)內(nèi)對應(yīng)P點(diǎn)位置，Q12，Q22，Q11和Q21為P點(diǎn)四周的DEM格網(wǎng) 點(diǎn)。欲得到P點(diǎn)的高程值，采用P點(diǎn)四周的四個(gè)格網(wǎng)點(diǎn)的高程值進(jìn)行雙線性插值。首先在x軸 方向上，對R1和R2兩個(gè)點(diǎn)進(jìn)行插值。然后根據(jù)R1和R2對P點(diǎn)進(jìn)行插值。
[0104] -般來講，欲求函數(shù)f在點(diǎn)P=(x，y)的值，假設(shè)我們已知函數(shù)f在Qn = (xi，yi)，Qi2 =(xi，y2)，Q21 = (X2，yi)及Q22 = (X2，y2)四個(gè)點(diǎn)的值，貝首先在x方向進(jìn)行線性插值，得到
[0106]然后在y方向進(jìn)行線性插值，得到
[0108]這樣就得到所要的結(jié)果f(x，y)，
[0110] 步驟3.3,將經(jīng)融合的DEM與數(shù)個(gè)外部控制點(diǎn)進(jìn)行統(tǒng)計(jì)配準(zhǔn)，所述配準(zhǔn)模型包括
[0111] 其中，0為控制點(diǎn)位置的局部坡向，a為控制點(diǎn)位置的局部坡度，盡為所有控制點(diǎn)位 置坡度的平均值，a為DEM平面偏移矢量的模，b為DEM偏移矢量的方向角，A Z為DEM垂直偏 移，A h為DEM任一點(diǎn)的高程值與控制點(diǎn)高程值的差。
[0112] 上述配準(zhǔn)模型的推導(dǎo)如下：
[0113] 假設(shè)InSAR生成的DEM的平面偏移矢量，AZs為平面偏移矢量的方位向，則平面偏移 矢量大小在坡向?yàn)?的山坡處的投影S xy的大小為(如圖2所示）：
[0114] Sxy = |S|cos(AZs-P) (22)
[0115] DEM某點(diǎn)位的高程與真實(shí)高程之間的差A(yù)h為(如圖3所示）：
[0116] A h = | Sxy | tan〇+ | Sz (23)
[0117] a為局部坡度，Sz為DEM在垂直向上的偏移矢量，上式可以進(jìn)一步轉(zhuǎn)為
(24)
[0119] 將上式簡化為
(25)
[0121] c等于
，則公式(24)成為一個(gè)干凈的余弦模型，Ah，a，P為觀測值，a、b、c為待 估模型參數(shù)。
[0122] 在所有控制點(diǎn)中選取一部分在坡度、坡向和空間上較為均勻分布的控制點(diǎn)作為參 與模型求解的解算點(diǎn)，剩余的點(diǎn)作為誤差改進(jìn)效果的檢查點(diǎn)。
[0123] 采用解算控制點(diǎn)高程差、坡度和坡向作為觀測值來進(jìn)行最小二乘平差求解配準(zhǔn)模 型參數(shù)，然后根據(jù)配準(zhǔn)模型參數(shù)對DEM進(jìn)行平移達(dá)到配準(zhǔn)目的，完成一次統(tǒng)計(jì)配準(zhǔn)后求取在 檢查控制點(diǎn)位高程差值的標(biāo)準(zhǔn)方差改善幅度。
[0124] 步驟3.4，對步驟3.1至3.3的配準(zhǔn)過程進(jìn)行迭代，直至高程差標(biāo)準(zhǔn)差改善率低于 2%〇
[0125] 因?yàn)椴襟E3.3的統(tǒng)計(jì)配準(zhǔn)中用Um泛替代tana會(huì)引入跟局部坡度相關(guān)的誤差。坡度 分布在0~90之內(nèi)，所以坡度小的區(qū)域配準(zhǔn)改善度偏小，而坡度大的區(qū)域配準(zhǔn)改善度過大。 需用一個(gè)與局部坡度相關(guān)的正切模型來消除這個(gè)引入的誤差：
[0126] A h = a7tan(a)+b7 (26)
[0127] 采用解算控制點(diǎn)位的坡度值與高程差作為觀測值，通過最小二乘平差方法來求解 待估參數(shù)a7和b7,然后將(26)所述的模型求取的誤差趨勢值與經(jīng)過步驟3.3改正后的DEM進(jìn) 行相減，得到進(jìn)一步改正后的DEM。
[0128] 步驟4,校正DEM與坡度、平面位置和坡向相關(guān)的系統(tǒng)誤差。
[0129] 步驟4.1，對DEM結(jié)果與平面位置相關(guān)的誤差進(jìn)行校正，通過二次曲面擬合來完成。 [0130]在InSAR差分干涉處理過程中的基線估計(jì)往往會(huì)有誤差，基線殘差在DEM產(chǎn)品的體 現(xiàn)為平面趨勢。因此可以通過對DEM結(jié)果進(jìn)行二次曲面擬合來消除這部分與平面位置相關(guān) 的誤差，二次曲面擬合模型包括
[0131] A h(x,y) =aix2+a2y2+a3xy+a4x+a5y+a6 (27)
[0132] 其中x和y表示某個(gè)像元的橫坐標(biāo)和縱坐標(biāo)，A hu,y)為該點(diǎn)的DEM值與控制點(diǎn)高程 差，ai……a6為待估模型參數(shù)。采用解算控制點(diǎn)位的位置與高程差作為觀測值，通過最小二 乘平差方法來求解待估參數(shù)ai……a 6,然后將公式(27)所述的模型求取的誤差趨勢值與經(jīng) 過步驟3改正后DEM進(jìn)行相減，得到進(jìn)一步改正后的DEM。
[0133] 步驟4.2,對DEM結(jié)果與坡度相關(guān)系統(tǒng)誤差進(jìn)行校正，通過二次或三次多項(xiàng)式模型 來完成。
[0134] 由于SAR斜距成像特點(diǎn)，使得SAR影像不可避免的出現(xiàn)透視收縮。不同平臺(tái)下得到 的DEM融合可以消除一部分透視收縮影響，但仍有可能具有不可忽視的殘留誤差。因?yàn)橥敢?收縮與局部坡度大小有直接關(guān)系，所以設(shè)定一個(gè)二次或三次多項(xiàng)式模型來消除，具體選擇 哪個(gè)模型視高程差的標(biāo)準(zhǔn)差改善效果而定，
[0135] 二次多項(xiàng)式為
[0136] A h(si〇Pe) = bia2+b2a+b3 (28)
[0137] 三次多項(xiàng)式為
[0138] A h(si〇Pe) = bia3+b2a2+b3a+b4 (29)
[0139] 其中a表示某個(gè)像元點(diǎn)的坡度值，A h^ic^)為該點(diǎn)的DEM值與控制點(diǎn)高程差h…… b4為待估t吳型參數(shù)，
[0140]采用解算控制點(diǎn)位的坡度值與高程差作為觀測值，通過最小二乘平差方法來求解 待估參數(shù)匕……b4,然后將公式(28)或(29)所述的模型求取的誤差趨勢值與經(jīng)過平面位置 誤差改正后的DEM進(jìn)行相減，得到進(jìn)一步改正后的DEM。
[0141] 步驟4.3,對DEM結(jié)果與坡向相關(guān)系統(tǒng)誤差進(jìn)行校正，通過正弦模型來完成的。
[0142] 除局部坡度以外，SAR影像透視收縮還與局部坡向有直接關(guān)系，所以設(shè)定一個(gè)正弦 模型來消除與坡向相關(guān)的系統(tǒng)誤差，所述正弦模型為
[0143] A haspect=Aisin( co0+<i) )+Ci (30)
[0144] 其中0為某個(gè)像元點(diǎn)的坡向值，A haspec；t為該點(diǎn)的DEM值與控制點(diǎn)高程差，Ai、《、<i> 和(^為待估模型參數(shù)。為方便最小二乘平差解算，該正弦模型可進(jìn)一步變換為一階傅立葉 模型
[0145] A hasPect=Aisin ?0+Bicos ?0+Ci (31)
[0146] 其中為某個(gè)像元點(diǎn)的坡向值，Ahaspec；t為該點(diǎn)的DEM值與控制點(diǎn)高程差，? 和&為待估模型參數(shù)，
[0147] 采用解算控制點(diǎn)位的坡向值與高程差作為觀測值，通過最小二乘平差方法來求解 待估參數(shù)和d，然后將公式(31)所述的模型求取的誤差趨勢與經(jīng)過坡度誤差改正 后的DEM進(jìn)行相減，得到進(jìn)一步改正后的DEM。
[0148] 本實(shí)施例所述方法充分融合了 InSAR監(jiān)測手段分辨率高、相對精度高、覆蓋范圍大 和地面控制點(diǎn)絕對精度高的特點(diǎn)，既解決InSAR獲取DEM絕對精度不高，又克服了傳統(tǒng)DEM獲 取方法例如星載和機(jī)載光學(xué)及激光測高手段獲取高分辨測圖效費(fèi)比低和受云層影響，地面 測圖費(fèi)時(shí)費(fèi)力且無法得到整個(gè)野外山區(qū)高分辨地貌等缺點(diǎn)。整個(gè)流程結(jié)構(gòu)清晰，具有實(shí)現(xiàn) 簡單、費(fèi)用低、監(jiān)測精度高、監(jiān)測范圍大、自動(dòng)化程度高等優(yōu)點(diǎn)。通過實(shí)地測試，經(jīng)本流程得 到的地形復(fù)雜區(qū)域DEM絕對精度可達(dá)到2m左右，其中72 %的點(diǎn)精度在1.5m以內(nèi)（如圖4所 示)。
[0149] 應(yīng)當(dāng)理解的是，對本領(lǐng)域普通技術(shù)人員來說，可以根據(jù)上述說明加以改進(jìn)或變換， 而所有這些改進(jìn)和變換都應(yīng)屬于本發(fā)明所附權(quán)利要求的保護(hù)范圍。
【主權(quán)項(xiàng)】
1. 一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，其特征在于，包括： 步驟1，通過雙基InSAR差分干涉處理單軌影像對，獲取單軌DEM; 步驟2，采取不同軌道的影像對分別獲取DEM，并進(jìn)行逐點(diǎn)融合； 步驟3，將經(jīng)融合的DEM與地面控制點(diǎn)進(jìn)行配準(zhǔn)； 步驟4,校正DEM與坡度、平面位置和坡向相關(guān)的系統(tǒng)誤差。2. 根據(jù)權(quán)利要求1所述的一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，其特 征在于，所述步驟1包括： 步驟1.1，根據(jù)DEM需求區(qū)域坐標(biāo)范圍定制影像對，所述影像對包括主影像和輔影像； 步驟1.2，將所述主影像和所述輔影像共輒相乘生成初始干涉圖； 步驟1.3通過軌道數(shù)據(jù)估計(jì)基線，根據(jù)基線由外部DEM模擬包含平地相位的地形相位 圖； 步驟1.4將所述地形相位圖與所述初始干涉相位圖相減得到差分干涉圖； 步驟1.5對差分干涉圖進(jìn)行相位濾波和相位解纏，將解纏相位轉(zhuǎn)化為高程差； 步驟1.6將高程差與外部DEM疊加得到新DEM; 步驟1.7通過地理編碼將所述新DEM從SAR坐標(biāo)系轉(zhuǎn)換到外部DEM的地理坐標(biāo)系得到所 述單軌DEM。3. 根據(jù)權(quán)利要求2所述的一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，其特 征在于，還包括:將所述步驟1.7所得所述單軌DEM作為所述外部DEM，對步驟1.1至步驟1.7 進(jìn)行迭代，直至獲取原始影像像元大小的DEM作為所述單軌DEM。4. 根據(jù)權(quán)利要求2所述的一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，其特 征在于，所述步驟1.5對差分干涉圖進(jìn)行相位濾波采用的濾波方式包括Nonlocal均值濾波 方式，濾波模型包括：s為濾波窗口內(nèi)的中心像元，t為濾波窗口內(nèi)的其他任意像元，真實(shí)的干涉相位為u，含 噪聲的原始干涉相位為Φ，相位估計(jì)值u( S )為濾波結(jié)果，W( S，t )為由像元t的鄰域像元與像 元S的鄰域像元之間的SAR影像強(qiáng)度和相干性的綜合相似度而決定的權(quán)重。5. 根據(jù)權(quán)利要求1所述的一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，其特 征在于，所述步驟2的具體步驟包括： 步驟21，將所述不同軌道的影像對統(tǒng)一配準(zhǔn)到相同坐標(biāo)系之下，所述不同軌道的影像 對包括升軌和降軌的影像對； 步驟22,采用融合模型進(jìn)行融合，所述融合模型包括其中，HfuslQn為融合后高程，83、13和仏分別為升軌的垂直基線值，疊影或陰影掩膜值及高程 值，Bd、Wd和Hd分別為降軌的垂直基線值，疊影或陰影掩膜值及高程值。6. 根據(jù)權(quán)利要求1所述的一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，其特 征在于，所述步驟3包括： 步驟3.1，將外部控制點(diǎn)的空間直角坐標(biāo)轉(zhuǎn)為大地經(jīng)煒度坐標(biāo)； 步驟3.2,采用雙線性插值提取外部控制點(diǎn)所在DEM高程值； 步驟3.3,將經(jīng)融合的DEM與數(shù)個(gè)外部控制點(diǎn)進(jìn)行統(tǒng)計(jì)配準(zhǔn)，所述配準(zhǔn)模型包括其中，β為控制點(diǎn)位置的局部坡向，α為控制點(diǎn)位置的局部坡度，泛為所有控制點(diǎn)位置坡 度的平均值，a為DEM平面偏移矢量的模，b為DEM偏移矢量的方向角，ΔΖ為DEM垂直偏移，Ah 為DEM任一點(diǎn)點(diǎn)的高程值與控制點(diǎn)高程值的差。7. 根據(jù)權(quán)利要求1或6所述的一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，其 特征在于:對所述配準(zhǔn)過程進(jìn)行迭代，直至高程差標(biāo)準(zhǔn)差改善率低于2 %。8. 根據(jù)權(quán)利要求1所述的一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，其特 征在于，所述步驟4對經(jīng)步驟三改正后DEM與平面位置相關(guān)的誤差進(jìn)行校正，是通過二次曲 面擬合來完成的，二次曲面擬合模型包括 Δ h(x,y) = aix2+a2y2+a3xy+a4X+a5y+a6 其中X和y表示任一個(gè)像元的橫坐標(biāo)和縱坐標(biāo)，△ h(x,y：i為所述像元點(diǎn)的DEM值與控制點(diǎn) 高程差，ai……a6為待估模型參數(shù)，采用解算控制點(diǎn)位的位置與高程差作為觀測值，通過最 小二乘平差方法來求解待估參數(shù)……a 6,然后將模型求取的誤差趨勢值與經(jīng)步驟3改正后 DEM進(jìn)行相減，得到進(jìn)一步改正后的DEM。9. 根據(jù)權(quán)利要求1所述的一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，其特 征在于，所述步驟4對DEM結(jié)果與坡度相關(guān)系統(tǒng)誤差進(jìn)行校正是通過二次或三次多項(xiàng)式模型 來完成的， 二次或三次多項(xiàng)式包括 Δ h(si〇pe) = bia2+b2a+b3或 Δ h(si〇pe) = bia3+b2a2+b3a+b4,其中α表 示某個(gè)像元點(diǎn)的坡度值，Ah(si。[^為該點(diǎn)的DEM值與控制點(diǎn)高程差bi......b4為待估模型參 數(shù)， 采用解算控制點(diǎn)位的坡度值與高程差作為觀測值，通過最小二乘平差方法來求解待估 參數(shù)bi……b4,然后將模型求取的誤差趨勢與經(jīng)過平面位置誤差改正后的DEM進(jìn)行相減，得 到進(jìn)一步改正后的DEM。10. 根據(jù)權(quán)利要求1所述的一種結(jié)合外控點(diǎn)的InSAR高精度高分辨率DEM獲取方法，其特 征在于，所述步驟4對DEM結(jié)果與坡向相關(guān)系統(tǒng)誤差進(jìn)行校正是通過正弦模型來完成的，所 述正弦模型為么]1 ￡1鄧如=418；[11(〇0+(|))+(：1 其中β為某個(gè)像元點(diǎn)的坡向值，Ahasp(3C；t為該點(diǎn)的DEM值與控制點(diǎn)高程差，Αι、ω、φ和Ci 為待估模型參數(shù)。 采用解算控制點(diǎn)位的坡向值與高程差作為觀測值，通過最小二乘平差方法來求解待估 參數(shù)ω、φ和d，然后將模型求取的誤差趨勢與經(jīng)過坡度誤差改正后的DEM進(jìn)行相減，得 到進(jìn)一步改正后的DEM。
【文檔編號(hào)】G01S13/90GK105929398SQ201610246461
【公開日】2016年9月7日
【申請日】2016年4月20日
【發(fā)明人】程正逢, 胡俊, 張健, 李佳, 陳功, 杜亞男, 周毅, 曾渠豐, 劉佳瑩
【申請人】中國電力工程顧問集團(tuán)中南電力設(shè)計(jì)院有限公司