雷達(dá)衛(wèi)星自動成圖的精密干涉測量關(guān)鍵技術(shù)

唐新明，李 濤，高小明，陳乾福，張 祥

國家測繪地理信息局衛(wèi)星測繪應(yīng)用中心，北京 100048

數(shù)字?jǐn)z影測量技術(shù)經(jīng)歷了從半數(shù)字時代到全數(shù)字時代的變革，正在由傳統(tǒng)光學(xué)攝影測量手段進(jìn)一步轉(zhuǎn)向?yàn)榫哂懈呔取⒏咧悄芑?、寬?nèi)涵化的多傳感器聯(lián)合測量技術(shù)?，F(xiàn)階段，以高自動化、高精度為主的智能精密干涉測量技術(shù)，已經(jīng)成為了數(shù)字?jǐn)z影測量技術(shù)的必要內(nèi)涵。精密干涉測量技術(shù)使用合成孔徑雷達(dá)干涉(interferometric SAR，InSAR)獲取地表的干涉相位，從而反演高精度高程信息，現(xiàn)已成為了最有效的全球測圖手段之一[1]。InSAR的高自動化體現(xiàn)在其干涉鏈路處理過程中。InSAR數(shù)據(jù)處理過程中采用具備高相干性的兩景單視復(fù)(single look complex，SLC)數(shù)據(jù)，完成自動化干涉、去平、解纏、編碼等操作，生成具有較高精度的數(shù)字高程模型，整個鏈路完全無需人工參與，在多節(jié)點(diǎn)并行條件下，可在一分鐘之內(nèi)完成整個鏈路的處理。InSAR的高精度體現(xiàn)在其干涉鏈路之前的參數(shù)設(shè)計(jì)與檢校過程，以及之后的DEM后處理過程。其中最重要的是參數(shù)的設(shè)計(jì)過程，即衛(wèi)星平臺與載荷穩(wěn)定性設(shè)計(jì)，測量檢校精度設(shè)計(jì)等。隨后的DEM后處理需要采用外部數(shù)據(jù)修正基線、斜距、軌道、整周未知數(shù)等殘差帶來的微小的高程梯度誤差及高程常數(shù)誤差，修正之后的數(shù)據(jù)能夠達(dá)到極高的精度要求。例如TanDEM-X的DEM數(shù)據(jù)修正完之后，其精度能夠由4 m提升到1.3 m(不包含南極地區(qū))。

由于精密干涉測量技術(shù)具有高智能化特性，因此多次被用于全球測圖任務(wù)中。早在2000年，美國國家航空航天局就啟動了航天飛機(jī)雷達(dá)地形測繪任務(wù)(shuttle radar topography mission，SRTM)，采用C波段SAR數(shù)據(jù)獲取了56°S至60°N的全球DEM數(shù)據(jù)，并于2003年發(fā)布了第一個版本，其標(biāo)稱格網(wǎng)大小為3″，約90 m，絕對高程精度優(yōu)于16 m[2]。雖然已經(jīng)過去十?dāng)?shù)年，但是SRTM依然是迄今為止在覆蓋范圍、高程精度、公開程度等各方面綜合評價最好的全球DEM數(shù)據(jù)。表1中給出了現(xiàn)有全球DEM的獲取技術(shù)及主要參數(shù)。其中最早的是美國地質(zhì)調(diào)查局(U.S.Geological Survey，USGS)于1996年發(fā)布的GTOPO30數(shù)據(jù)，其格網(wǎng)大小為30″，約900 m，高程精度優(yōu)于30 m。其后，美國國家海洋和大氣管理局在2008年發(fā)布了全球1′的地形模型ETOPO1，其中包含了冰蓋及基巖數(shù)據(jù)。然而ETOPO1的數(shù)據(jù)來源更加復(fù)雜多樣，無法給出具體的精度[3]。USGS隨后發(fā)布了GMTED2010(global multi-resolution terrain elevation data 2010)，這是由11種數(shù)據(jù)融合得到的DEM，其產(chǎn)品包含3種分辨率，即30″、15″及7.5″，精度優(yōu)于10 m[4]。2009年ASTER GDEM發(fā)布了第一個版本，然而由于質(zhì)量不佳，只能支撐部分科研應(yīng)用，隨后在2011年，ASTER GDEM發(fā)布了第二個版本，其標(biāo)稱格網(wǎng)大小為15 m，絕對高程精度接近像素尺寸大小[5]。2015年，JAXA發(fā)布了ALOS DSM，其標(biāo)稱精度為15 m。2016年，DLR獲取了全球包括南極地區(qū)的DEM數(shù)據(jù)，其高程精度達(dá)到了驚人的4 m以內(nèi)[6]，為全球測圖拉開了新的序幕，然而DLR獲取的DEM數(shù)據(jù)由Airbus進(jìn)行代理售賣，已經(jīng)不再是完全公開的數(shù)據(jù)。

表1 全球DEM及其主要參數(shù)

從表1可以看出，除了多數(shù)據(jù)源融合之外，獲取全球DEM的技術(shù)主要包括兩類，即光學(xué)立體攝影測量和InSAR。其中光學(xué)立體攝影測量技術(shù)使用具有一定視差的像對完成高程解算，其視差大小一般不小于30°，并確?；弑仍?左右，從而平衡同名點(diǎn)匹配精度和立體幾何的穩(wěn)健性。現(xiàn)階段，我國第一顆民用立體測繪衛(wèi)星資源三號，已經(jīng)完成了國內(nèi)近400萬平方千米的高精度DEM生產(chǎn)[7]。DOM更是覆蓋了924萬平方千米，有效覆蓋率達(dá)到96%。然而在西南地區(qū)，由于全年都有云霧覆蓋，始終無法獲得有效數(shù)據(jù)。因此我國對干涉SAR衛(wèi)星的需求日益迫切。InSAR技術(shù)使用視差極小的像對完成高程解算，其視差大小不會超過0.1°，并確保模糊高在35～55 m之間，從而平衡相位解纏精度和干涉幾何的穩(wěn)健性。相對于光學(xué)來說，InSAR無需進(jìn)行逐像素匹配，粗差較少，不受云霧干擾，且能夠直接提供逐點(diǎn)的高程信息，無需借助三角網(wǎng)構(gòu)建DEM，因此有極大的測繪優(yōu)勢。然而其陰影、疊掩、穿透性等問題，則需要在任務(wù)設(shè)計(jì)過程中予以高度重視?，F(xiàn)階段，我國雖然已經(jīng)先后發(fā)射了環(huán)境一號和高分三號兩顆民用SAR衛(wèi)星，然而受限于軌道回歸過程中垂直基線的控制精度，衛(wèi)星并不具備較好的業(yè)務(wù)化干涉能力。在“十三五”期間，我國規(guī)劃了第一對以干涉為主要任務(wù)的民用SAR衛(wèi)星，即L波段差分干涉SAR衛(wèi)星。衛(wèi)星擬采用雙星繞飛模式獲取DEM數(shù)據(jù)，為我國全球測圖開啟新的篇章。

InSAR獲取全球DEM的技術(shù)在國內(nèi)外諸多文獻(xiàn)中已經(jīng)有詳細(xì)介紹。國外以NASA和DLR為主的科研團(tuán)隊(duì)已經(jīng)公開了SRTM及TanDEM-X的精密干涉測量指標(biāo)設(shè)計(jì)方法[8-10]、幾何參數(shù)標(biāo)定方法[8，11-13]、數(shù)據(jù)處理方法[14-15]、精度的詳細(xì)評估[2]等。國內(nèi)也有學(xué)者針對精密干涉測量開展了衛(wèi)星性能分析[16]、數(shù)據(jù)處理方法研究[17]等。然而SRTM的原始SLC數(shù)據(jù)并不公開發(fā)布，TanDEM-X的SLC數(shù)據(jù)也必須通過科研途徑申請獲取，因此數(shù)據(jù)來源的限制使得國內(nèi)針對雙天線及分布式SAR衛(wèi)星的研究多處于模擬仿真階段[18]，或者使用機(jī)載數(shù)據(jù)替代星載數(shù)據(jù)進(jìn)行先期研究[19-22]。本文從精密干涉測量的全球測圖應(yīng)用出發(fā)，給出了基于分布式SAR衛(wèi)星的原始數(shù)據(jù)到最終DEM產(chǎn)品的整個技術(shù)鏈路，并闡述了SAR衛(wèi)星的測量檢校技術(shù)、數(shù)據(jù)處理技術(shù)及DEM數(shù)據(jù)后處理技術(shù)，為數(shù)字?jǐn)z影測量的自動化和智能化奠定了基礎(chǔ)，并為國產(chǎn)SAR衛(wèi)星的全球測圖任務(wù)提供技術(shù)參考。

1 測量檢校技術(shù)

SAR衛(wèi)星參數(shù)是數(shù)據(jù)處理的必要輸入，參數(shù)的精度決定了數(shù)據(jù)處理的直接精度。然而，高精度測量參數(shù)無法通過衛(wèi)星直接給定。例如基線參數(shù)是干涉過程中的重要參數(shù)，毫米級的基線誤差會帶來米級的高程誤差。衛(wèi)星的基線測量設(shè)備直接給出的是設(shè)備的幾何中心的位置，而并非天線相位中心的位置，天線相位中心基線參數(shù)則需要聯(lián)合人工角反射器進(jìn)行檢校。因此，檢校技術(shù)是確定星上參數(shù)誤差并進(jìn)行誤差修正的必要途徑。SAR衛(wèi)星的檢校鏈路包括內(nèi)定標(biāo)、天線指向標(biāo)定、輻射定標(biāo)、幾何檢校、干涉測量檢校等。其中幾何檢校和干涉測量檢校技術(shù)確定了應(yīng)用過程中所能獲取的干涉幾何參數(shù)的精度，二者共同決定了全球DEM的直接精度水平。

1.1 幾何檢校技術(shù)

幾何檢校的目的是檢驗(yàn)并校正與幾何定位相關(guān)的SAR系統(tǒng)參數(shù)。圖1中給出了SAR干涉過程中的基本幾何關(guān)系，其中P為地面點(diǎn)，S1為主影像成像時刻的相位中心位置，S2為從影像成像時刻的相位中心位置，B為雙星之間的基線長度，α為基線傾角，B⊥為垂直基線長度，θ為側(cè)視角，RH為主影像相位中心到地心的距離，H為衛(wèi)星高度，r為P點(diǎn)對應(yīng)的斜距，Re為地球曲率半徑，h為P點(diǎn)的高程，β為本地入射角，Rh為P點(diǎn)到地心的距離。為了便于描述，所有的參數(shù)統(tǒng)一在WGS-84坐標(biāo)系中進(jìn)行探討。SAR一般采用斜距-多普勒(Range-Doppler，R-D)模型完成像素坐標(biāo)與地理坐標(biāo)的轉(zhuǎn)換，即

(1)

圖1 干涉幾何示意圖Fig.1 Diagrammatic drawing of interferometric geometry

斜距的誤差來源主要包括兩部分，一部分為大氣誤差，另一部分是距離向計(jì)時系統(tǒng)誤差。其中大氣誤差包含流體靜力學(xué)分量(氣壓)、濕分量(水汽)及流體力學(xué)分量(水滴)。對于TanDEM-X來說，三類分量在天頂方向的延遲可采用經(jīng)驗(yàn)?zāi)Ｐ徒票磉_(dá)[23]

(2)

當(dāng)高程小于9000 m時，經(jīng)驗(yàn)?zāi)Ｐ团c精確的大氣模型之間的差異一般在0.7 m以內(nèi)，對于1∶50 000比例尺要求的25 m平面定位精度來說，其誤差在允許的范圍之內(nèi)。星上距離向計(jì)時系統(tǒng)誤差一般來說是數(shù)百ns，這會帶來數(shù)十米的斜距測量誤差，不過這部分誤差是系統(tǒng)誤差，可通過與式(2)類似的過程加以檢校，即

(3)

檢校后的時間誤差基本在1 ns左右，能夠完成無控條件下的1∶50 000比例尺定位要求。

SAR衛(wèi)星一般需要在星上采用偏航導(dǎo)引，并在地面采用零多普勒成像，以確保SLC中多普勒中心頻率會在0 Hz附近。實(shí)際過程中多普勒中心頻率無法嚴(yán)格歸零，其誤差一般在數(shù)個赫茲，對于ScanSAR成像模式來說，中心頻率誤差最大，不過其限差也會控制在10 Hz以內(nèi)。在實(shí)際處理過程中，多普勒中心頻率誤差為時間域高頻分量，無法精確標(biāo)定，因此一般不做過多處理。除了多普勒中心頻率之外，多普勒方程中的未知參數(shù)還包括主星位置及速度矢量。它依賴于星上GPS精度，現(xiàn)階段雙頻GPS可提供優(yōu)于5 cm的事后定軌精度，以及mm/s級的速度測量精度。在這種情況下，軌道的誤差大部分來源于方位向時間誤差，如果軌道運(yùn)行在500 km高度，那么1 ms的方位向時間誤差會引入約7 m的定軌偏差，這對于精密定位來說，是不可忽略的誤差源。方位向時間誤差Δta一般通過地面角反射器進(jìn)行修正，即

(4)

幾何檢校過程中需要使用人工角反射器確保地面點(diǎn)的坐標(biāo)與影像坐標(biāo)嚴(yán)格對應(yīng)，從而減少刺點(diǎn)誤差，提高檢校精度。理論上來說，一個角反射器就可以完成檢校工作，但是為了減少粗差的影響，一般至少采用6個角反射器，采用兩行三列的方式均勻分布在影像刈幅范圍內(nèi)，同時為了減小其他誤差的影響，檢校區(qū)域的高程變化應(yīng)盡量小。幾何檢校之后，需要保證SAR影像在全球范圍內(nèi)無控定位精度達(dá)到25 m以內(nèi)，且衛(wèi)星參數(shù)應(yīng)至少在3個月內(nèi)保持極好的穩(wěn)定性，從而支撐全球測圖任務(wù)。

1.2 干涉測量檢校技術(shù)

干涉測量檢校的目的是修正基線誤差，確保無控條件下的高程精度滿足指定標(biāo)準(zhǔn)。從資源三號的成功應(yīng)用來看，衛(wèi)星測繪的精度一般不低于1∶50 000，即丘陵地區(qū)的高程精度優(yōu)于5 m。對于InSAR全球測圖來說，一般也遵循類似的標(biāo)準(zhǔn)，即SAR衛(wèi)星需要在完成干涉測量檢校及后處理之后，保證全球高程精度達(dá)到5 m。這對于我國InSAR全球測圖來說，是極為嚴(yán)苛的精度水平。首先，我國民用干涉SAR衛(wèi)星剛剛起步，雖然已經(jīng)有了高分三號衛(wèi)星，衛(wèi)星成像質(zhì)量也達(dá)到了極高的水平，但是高分三號尚不具備業(yè)務(wù)化干涉能力，無法用于全球測圖。其次，我國L-SAR雖然是以干涉為主的衛(wèi)星，但是依然是科研星，部分研制指標(biāo)相比TanDEM-X依然有較大的差距，在這種情況下，進(jìn)行可靠的干涉測量檢校，測量和修正干涉基線誤差，就顯得尤為重要。

SAR衛(wèi)星基線測量方式與干涉模式相關(guān)?，F(xiàn)階段InSAR地形測繪主要通過兩種干涉模式實(shí)現(xiàn)，一種是以SRTM為代表的雙天線干涉模式，這種情況下，基線需采用姿軌測量儀(Attitude and Orbit Determination Avoinics，AODA)進(jìn)行測量，AODA提供的基線數(shù)據(jù)精度約2 mm，姿態(tài)精度約9″，能夠滿足DETD-2的高程精度要求(16 m)。另一種是以TanDEM-X為代表的雙星編隊(duì)干涉模式，這種情況下，基線需采用星間雙差GPS(double differential GPS，DDGPS)進(jìn)行測量，測量精度要求達(dá)到1 mm，從而滿足HRTI-3的高程精度要求(4 m)。

然而值得說明的是，無論何種測量方式，獲取的始終是天線的幾何中心對應(yīng)的基線，而并非SAR天線相位中心對應(yīng)的基線。在InSAR地形測繪過程中使用的則是后者，因此在實(shí)際測繪過程中，依然需要對基線參數(shù)進(jìn)行標(biāo)定，消除基線參數(shù)測量和轉(zhuǎn)換過程中的誤差。L-SAR采用的雙星編隊(duì)繞飛模式，其基線測量方式與TanDEM-X相似。DDGPS直接獲取的參數(shù)是基線的三維分量，即[BxByBz]，在InSAR地形測繪過程中，一般采用更為直觀的基線分量[BTB⊥B‖]，即沿軌基線、垂直基線及平行基線。其中沿軌基線并不影響InSAR交軌干涉的地形測繪精度，一般不予考慮。垂直基線誤差對高程誤差的傳播可描述為

(5)

即垂直基線誤差傳播過程與地表高程相關(guān)。因此檢校過程中應(yīng)盡量選擇平坦地區(qū)，以避免高程差異引入額外的誤差。平行基線誤差對高程誤差的傳播分為兩個方面，一方面，平行基線會通過誤差傳播直觀傳遞到高程誤差中，即

(6)

(7)

式中，Δs為地距向兩點(diǎn)間的距離；Δht為兩點(diǎn)間的相對高程誤差，在基線傾角小于側(cè)視角的情況下，側(cè)視角越大，垂直基線越小，平行基線誤差帶來的高程誤差越大?；谄叫谢€的這種誤差特性，一般需要采用相隔上千千米的高程較為一致的區(qū)域進(jìn)行平行基線誤差標(biāo)定。

2 數(shù)據(jù)處理技術(shù)

衛(wèi)星下傳的數(shù)據(jù)是回波數(shù)據(jù)，數(shù)據(jù)處理系統(tǒng)需完成零多普勒成像及距離徙動校正等，生成經(jīng)嚴(yán)格配準(zhǔn)的SLC數(shù)據(jù)。同時使用SLC完成干涉鏈路的處理，生成高精度的經(jīng)編碼的DEM初始產(chǎn)品。干涉鏈路的處理過程中，嚴(yán)格來講不應(yīng)采用任何外部數(shù)據(jù)，即可將高程誤差控制在一定精度標(biāo)準(zhǔn)以內(nèi)(一般為最終產(chǎn)品精度的2～3倍)。這對于數(shù)據(jù)處理的精度和效率提出了極大的挑戰(zhàn)，特別是高精度相位解纏技術(shù)及快速地理編碼技術(shù)。

2.1 高精度相位解纏技術(shù)

高精度相位解纏技術(shù)需要同時確保解纏結(jié)果的穩(wěn)健性和可靠性。首先，算法應(yīng)能夠滿足不同地形條件下的相位解纏精度需求，不會由于地形的變化或者解纏參數(shù)的變化出現(xiàn)較大的解纏精度差異。其次，算法需能夠確保相位的可靠性，即解纏相位應(yīng)與真實(shí)相位保持最小的L-P范數(shù)。一般來說，在工程中常采用1-P范數(shù)，即最小費(fèi)用流進(jìn)行相位解纏，以平衡穩(wěn)健性、可靠性與解纏效率。解纏之后的相位常數(shù)需要采用雷達(dá)攝影測量的手段解算，即獲取主輔影像同名點(diǎn)之間的雷達(dá)攝影測量時間延遲Δts和Δtm，從而用下式解算獲得[15]

(8)

式中，f0為雷達(dá)中心頻率。由于時間延遲差異為未知量，一般在配準(zhǔn)過程中采用配準(zhǔn)偏移量確定。推薦在整景影像中采用32×32個點(diǎn)進(jìn)行解算，以便同時兼顧精度和效率。在獲取了每個點(diǎn)的配準(zhǔn)偏移量，從而解算出絕對相位常數(shù)之后，應(yīng)使用穩(wěn)健估計(jì)算法，獲取最終的絕對相位對應(yīng)的整周未知數(shù)，從而解算參考點(diǎn)的絕對相位。

2.2 快速地理編碼技術(shù)

快速地理編碼技術(shù)是確保整個數(shù)據(jù)處理鏈路效率的重要環(huán)節(jié)。地理編碼算法一般分為前向編碼和后向編碼。前向編碼通過R-D模型將像素坐標(biāo)解算到地理坐標(biāo)，解算之后的DEM數(shù)據(jù)并不具備規(guī)則的格網(wǎng)，因此一般采用后向編碼。后向編碼從規(guī)則的DEM格網(wǎng)中獲取地理坐標(biāo)，通過R-D模型解算到對應(yīng)的像素坐標(biāo)，從而獲取對應(yīng)高程值，迭代優(yōu)化模型參數(shù)，獲取精確的點(diǎn)位，最后以16點(diǎn)sinc截斷函數(shù)進(jìn)行點(diǎn)位高程插值，獲取可靠的點(diǎn)位高程信息。如果不進(jìn)行算法優(yōu)化，單線程條件下30×50 km2的范圍內(nèi)完成10 m格網(wǎng)插值會消耗長達(dá)二三十分鐘的時間，這會大大降低全球測繪的效率，因此需要進(jìn)行工程優(yōu)化，將編碼時間控制在1 min以內(nèi)。編碼之后的DEM為初始DEM，這些殘差的大小一般不超過10 m，需要在后續(xù)的區(qū)域網(wǎng)平差及多數(shù)據(jù)后處理過程中進(jìn)一步消除。

3 數(shù)據(jù)后處理技術(shù)

初始DEM的高程精度雖然已經(jīng)能夠控制在較好的水平，但是對于1∶50 000比例尺地形圖測繪來說，需達(dá)到平地、丘陵、山地和高山地分別為3 m、5 m、8 m及14 m的高程精度水平，還需要進(jìn)一步的數(shù)據(jù)處理。首先，需借助區(qū)域網(wǎng)平差技術(shù)消除部分衛(wèi)星參數(shù)殘差；其次，需使用多軌數(shù)據(jù)進(jìn)行聯(lián)合處理，減小低相干區(qū)域面積，提高DEM的有效覆蓋范圍。

3.1 DEM區(qū)域網(wǎng)平差技術(shù)

幾何和干涉測量檢校并不能消除所有誤差，初始DEM中將包含基線殘差、斜距殘差、定軌殘差及絕對相位解算殘差等。這些殘差量呈現(xiàn)一定的系統(tǒng)性，可在鑲嵌拼接過程中，通過連接點(diǎn)進(jìn)行平差消除。為了提高DEM精度，需要采用適量的地面控制點(diǎn)完成區(qū)域網(wǎng)平差過程。一般來說，各類殘差在影像中可采用一階多項(xiàng)式進(jìn)行擬合，其多項(xiàng)式表達(dá)為[24]

he=a+bx+cy

(9)

式中，a、b、c為多項(xiàng)式系數(shù)；x、y為距離向及方位向坐標(biāo)。在平差過程中，可用的控制數(shù)據(jù)主要包括三類，即ICESat-GLAS數(shù)據(jù)、SRTM數(shù)據(jù)及資源三號控制點(diǎn)庫數(shù)據(jù)。由于L-SAR的設(shè)計(jì)精度為1∶50 000，在高山地的精度需要達(dá)到14 m，因此SRTM(16 m)并不能滿足平差的參考數(shù)據(jù)精度要求。而資源三號控制點(diǎn)庫只能覆蓋國內(nèi)區(qū)域，且點(diǎn)位密度不均，不能確保參數(shù)估計(jì)的穩(wěn)健性。因此建議采用ICESat-GLAS數(shù)據(jù)完成平差處理，ICESat-GLAS的標(biāo)稱精度為2 m，可用于高程控制。然而ICESat-GLAS的光斑直徑為70 m，L-SAR的分辨率為米級，平差過程中需要將InSAR生成的DEM對應(yīng)的直徑為70 m圓覆蓋范圍內(nèi)的點(diǎn)的高程進(jìn)行平均，以建立合理的對應(yīng)關(guān)系。

3.2 長短基線組合技術(shù)

長短基線組合算法的目的是提高相位解纏的可靠性[25]。相位解纏精度與相位梯度相關(guān)，按照奈奎斯特采樣定律，當(dāng)相位梯度超過π時無法恢復(fù)真實(shí)相位。在實(shí)際相位解纏過程中，最大可探測相位梯度(maximum detectable phase gradient，MDPG)與多視數(shù)和相干系數(shù)之間的關(guān)系可表達(dá)為[26]

(10)

3.3 升降軌融合技術(shù)

SAR衛(wèi)星的側(cè)視特性會引入其特有的陰影、疊掩等畸變，因此使用InSAR技術(shù)進(jìn)行全球測圖的過程中，至少需要對地面進(jìn)行兩次觀測，以彌補(bǔ)測試過程中陰影和疊掩帶來的信息損失。例如SRTM完成了兩次對地全覆蓋，TanDEM-X完成了四次對地覆蓋(兩次全覆蓋，兩次重點(diǎn)地區(qū)重復(fù)覆蓋)。在一軌之內(nèi)，獲取升軌數(shù)據(jù)和降軌數(shù)據(jù)的時間各占一半。升軌數(shù)據(jù)提供指定區(qū)域由西向東的側(cè)視信息，降軌數(shù)據(jù)提供指定區(qū)域由東向西的側(cè)視信息。升降軌數(shù)據(jù)的陰影和疊掩區(qū)域基本互補(bǔ)，數(shù)據(jù)融合后可提高DEM數(shù)據(jù)的完整性。此外，對其他地區(qū)的數(shù)據(jù)來說，采用相干系數(shù)進(jìn)行加權(quán)計(jì)算，能夠在一定程度上剔除粗差，進(jìn)一步提高DEM精度。

4 試驗(yàn)分析

4.1 試驗(yàn)數(shù)據(jù)選擇

試驗(yàn)過程中采用了德國宇航中心(German Aerospace Center，DLR)提供的陜西省的TanDEM-X數(shù)據(jù)，數(shù)據(jù)是經(jīng)測量檢校處理的配準(zhǔn)的斜距復(fù)圖像(coregistered single look slant range complex，CoSSC)。研究區(qū)在西安市的東北側(cè)，涵蓋了閻良區(qū)、富平縣、蒲城縣及白水縣等。整個區(qū)域南側(cè)地形起伏極小，屬于典型的平地地形，東北側(cè)有一定的地形起伏，符合丘陵地的特征，而西北側(cè)地形起伏較大，符合山地的特征。試驗(yàn)過程中用到的數(shù)據(jù)共6景。其中成像日期為20130903(YYYYMMDD)和20130925的四景影像為降軌數(shù)據(jù)，20121021和20140402為升軌數(shù)據(jù)。影像的方位向采樣間隔為2.2 m，距離向采樣間隔為1.4 m，入射角為42.6°，對應(yīng)地距采樣間隔為2.1 m，刈幅大小為35.5×56.9 km2，其他概要信息如表2所示。由于影像的西北側(cè)地形起伏極其復(fù)雜，因此需聯(lián)合采用升降軌融合以及長短基線組合算法進(jìn)行聯(lián)合解算。詳細(xì)的解算過程如4.2節(jié)所示。

圖2 研究區(qū)及所用影像分布圖Fig.2 Images distribution map of the research regions

4.2 試驗(yàn)結(jié)果及分析

由于DLR已經(jīng)對CoSSC影像進(jìn)行了處理，分發(fā)的是檢校之后的數(shù)據(jù)，因此研究過程中無法對結(jié)果進(jìn)行進(jìn)一步的檢校處理。然而通過數(shù)據(jù)處理過程中可對DLR的檢校結(jié)果進(jìn)行側(cè)面的驗(yàn)證。首先，4景影像的無控地面定位精度可通過與SRTM的配準(zhǔn)確定。以20130903_3為例，其距離向的常數(shù)項(xiàng)偏移為2.65像素，方位向常數(shù)項(xiàng)偏移為0.12像素，處理過程中為了最大程度的抑制相干斑，采用了8×8的多視數(shù)，多視之后的像素采樣間隔為16 m，這意味著影像距離向和方位向直接定位精度分別為42.4 m及1.92 m，解算得到大氣的影響為5.52 m，對應(yīng)的距離向定位誤差為8.16 m，那么由斜距延遲帶來的誤差至少為34.2 m，并不能滿足高精度直接定位要求。處理過程中為了獲取更精確幾何位置，將影像與10倍過采樣之后的SRTM數(shù)據(jù)做了精配準(zhǔn)，其配準(zhǔn)精度為0.1像素以內(nèi)，以確保平面定位精度不超過0.9 m。同時，計(jì)算過程中采用精配準(zhǔn)之后的SRTM均勻選取了10 000個高程點(diǎn)進(jìn)行基線的精估計(jì)，精估計(jì)后的交軌方向和垂軌方向基線變化量分別為11.9 cm和5.7 cm，這意味著采用TanDEM-X提供的CoSSC數(shù)據(jù)只能獲取與定軌精度(5 cm)同數(shù)量級的基線精度，在實(shí)際的數(shù)據(jù)處理過程中，依然需要采用地面控制點(diǎn)完成高精度基線的估計(jì)。

表2 試驗(yàn)所用影像概要信息

CoSSC數(shù)據(jù)處理過程中采用的常規(guī)干涉處理流程，即多視、干涉、去平、濾波、相位解纏、相高轉(zhuǎn)換、地理編碼等。其中的關(guān)鍵技術(shù)之一就是無控條件下整周未知數(shù)的解算技術(shù)。論文從CoSSC中均勻選取了1024個同名點(diǎn)，通過配準(zhǔn)偏移量解算時間延遲和整周未知數(shù)，并求取整周未知數(shù)的眾數(shù)，以獲取穩(wěn)健結(jié)果。結(jié)果表明，這6景影像的整周未知數(shù)分別為π、12π、8π、11π、0π、3π弧度。在解纏結(jié)果中添加整周未知數(shù)，并進(jìn)行相高轉(zhuǎn)換以及地理編碼，即可得到初始DEM數(shù)據(jù)。

本文采用了SRTM、ICESat-GLAS數(shù)據(jù)及資源三號控制點(diǎn)庫數(shù)據(jù)進(jìn)行了初始DEM的精度評估，不同影像的高程精度如表2所示。其中使用SRTM及ICESat-GLAS評估的結(jié)果較為一致，精度會同步提升或下降，而控制點(diǎn)庫數(shù)據(jù)量較少，并不能全面反映高程的精度水平，綜合考慮下，推薦使用SRTM進(jìn)行InSAR處理鏈路的精度評估(從回波數(shù)據(jù)到初始DEM數(shù)據(jù))，使用GLAS對區(qū)域網(wǎng)平差、升降軌融合及長短基線組合等中間結(jié)果進(jìn)行精度的二次評估，而GCP適合對最終的DEM產(chǎn)品進(jìn)行精度評估，給出最終的產(chǎn)品精度。圖3給出了ICESat-GLAS點(diǎn)及控制點(diǎn)庫數(shù)據(jù)的分布情況。其中包含激光點(diǎn)17 439個，控制點(diǎn)58個。激光點(diǎn)在20130925_6中分布較少，只有411個，且集中在影像的西南方向，從分布上來看，其精度的評估并不穩(wěn)健。

其中深色為ICESat-GLAS點(diǎn)，淺色為DEM平差連接點(diǎn)，黑色為資三控制點(diǎn)庫數(shù)據(jù),多邊形選定區(qū)域?yàn)樯弟壓烷L短基線共同覆蓋的區(qū)域。圖3 試驗(yàn)用到的矢量點(diǎn)數(shù)據(jù)分布情況Fig.3 Vector data distribution map of the experimental regions

DEM區(qū)域網(wǎng)平差針對的是具有相同側(cè)視角度、相同飛行方向的數(shù)據(jù)，因此試驗(yàn)過程中使用表3 中的后4景影像開展了平差試驗(yàn)。平差過程中所用連接點(diǎn)為西側(cè)影像的最后一列(升軌情況下)或者北側(cè)影像的最后一行(降軌情況下)，首先在影像邊緣均勻選取100個連接點(diǎn)，并使用R-D模型解算連接點(diǎn)的地理坐標(biāo)信息，隨后將地理坐標(biāo)轉(zhuǎn)換為鄰接影像的像素坐標(biāo)，最終完成像素坐標(biāo)空間的平差。4景影像的平差結(jié)果如表4所示，由于20130925_6的點(diǎn)位分布較為集中，出現(xiàn)了秩虧問題，并未解算成功。由影像的平差結(jié)果可知，DEM結(jié)果在距離向及方位向均會有2～6 m/萬像素，即10～30 cm/km的高程斜坡。這種高程斜坡除了基線誤差的貢獻(xiàn)之外，還包括斜距、定軌誤差、相位整周未知數(shù)估算誤差等。在進(jìn)行全球測繪的過程中，如果能夠進(jìn)行精確的斜距和軌道標(biāo)定，并解決無控條件下整周未知數(shù)精估計(jì)問題，將有望進(jìn)一步縮小高程斜坡對應(yīng)的誤差大小。4景影像平差之后，使用GLAS點(diǎn)進(jìn)行了精度的重新評估，其高程精度分別提升到1.79 m、0.46 m、1.06 m及0.67 m，這說明平差算法有效提高了DEM的高程精度水平。

表3影像DEM精度評估結(jié)果

Tab.3DEMaccuracygenereatedfromsingleinterferogram

影像編號SRTMICESat-GLAS控制點(diǎn)庫20140402_13.125.961.2820121021_22.424.753.4420130903_31.394.723.7520130925_40.983.140.7420130903_51.031.671.8020130925_60.651.951.36

表4 DEM區(qū)域網(wǎng)平差解算結(jié)果

由于影像西北區(qū)域地形復(fù)雜，精度較差，多為3～5 m(表2)，因此論文針對這部分?jǐn)?shù)據(jù)進(jìn)行了長短基線組合及升降軌融合試驗(yàn)。試驗(yàn)結(jié)果如圖4 所示。由于第一景及第二景影像的側(cè)視幾何較為一致，因此本文使用這兩景影像進(jìn)行長短基線組合試驗(yàn)。試驗(yàn)過程中首先完成短基線的相位解纏及相高轉(zhuǎn)換，然后將短基線高程數(shù)據(jù)配準(zhǔn)到主影像空間中，進(jìn)而完成主影像干涉圖的去平去地形操作，減小相位梯度，提高相位解纏精度。論文使用GLAS數(shù)據(jù)對長短基線的公共區(qū)域進(jìn)行了精度評估，長基線的精度為5.90 m，短基線的精度為3.77 m，經(jīng)長短基線組合之后，精度達(dá)到5.04 m，長基線精度有較為明顯的提升。進(jìn)一步融合降軌數(shù)據(jù)之后，精度變?yōu)?.07 m，此時精度雖然已經(jīng)沒有提升，但是低相干區(qū)域面積已經(jīng)極大降低，統(tǒng)計(jì)結(jié)果表明，升軌數(shù)據(jù)經(jīng)長短基線組合之后，其低相干區(qū)域(經(jīng)多視及濾波之后相干性<0.9)的面積為253.6 km2，而升降軌融合之后的低相干面積僅為0.8 km2，這對于實(shí)現(xiàn)全球測圖中DEM的全覆蓋來說具有重要意義。

其中(a)—(e)的DEM分別對應(yīng)長基線、短基線、長短基線組合、降軌數(shù)據(jù)、長短基線和降軌數(shù)據(jù)融合結(jié)果。圖4 升降軌以及長短基線組合試驗(yàn)Fig.4 Results of the descending-and-ascending，long-and-short baseline combination experiments

5 結(jié)論與展望

本文提出了面向全球測圖的智能化精密干涉測量系統(tǒng)技術(shù)，并采用6景覆蓋陜西區(qū)域的TanDEM-X數(shù)據(jù)進(jìn)行了技術(shù)的處理與驗(yàn)證。首先，論文對數(shù)據(jù)進(jìn)行了定標(biāo)精度的評估。由于缺乏角反射器數(shù)據(jù)，本文采用SRTM數(shù)據(jù)對定位精度進(jìn)行了側(cè)面驗(yàn)證。結(jié)果表明，TanDEM-X在距離向的誤差較大，為40 m以上，而方位向能夠控制在2 m以內(nèi)。因此，在后續(xù)處理過程中，為了確保高程精度，采用SRTM數(shù)據(jù)為參考。其次，論文完成了數(shù)據(jù)的干涉處理，獲取了初始DEM，并開展了區(qū)域網(wǎng)平差試驗(yàn)，將DEM高程精度提升了0.61～2.93 m不等。隨后，本文進(jìn)行了長短基線和升降軌融合的試驗(yàn)分析。試驗(yàn)表明，長短基線組合算法能夠?qū)㈤L基線精度從5.9 m提升到5.04 m，進(jìn)一步融合降軌數(shù)據(jù)之后，精度并沒有明顯的提升，然而低相干區(qū)域的面積已經(jīng)從253.6 km2降低到0.8 km2，這極大提升了DEM數(shù)據(jù)的有效性。在數(shù)據(jù)處理及區(qū)域網(wǎng)平差過程中，不包含人工交互過程，這意味著使用InSAR進(jìn)行全球測圖時，除了DEM編輯過程之外，可極大提高整個鏈路的智能化水平，達(dá)到高精度、高效率、高智能化的攝影測量新要求。本文處理結(jié)果表明，TanDEM-X可提供5 m左右的高精度DEM數(shù)據(jù)，這能夠滿足我國1∶50 000比例尺山地地形8 m的測繪精度，能夠基本滿足我國1∶25 000比例尺山地地形5 m的測繪精度。

在使用InSAR進(jìn)行全球測圖的過程中，除了本文重點(diǎn)闡述的測繪精度之外，還需解決測繪效率問題、測繪標(biāo)準(zhǔn)問題及測繪產(chǎn)品生產(chǎn)問題。一方面，需要確保在多節(jié)點(diǎn)并行運(yùn)算情況下，從回波數(shù)據(jù)開始，完成整個干涉鏈路處理，生成經(jīng)編碼的初始DEM數(shù)據(jù)，能夠?qū)r間控制在1 min左右。另一方面，亟需針對L波段差分干涉SAR衛(wèi)星開展權(quán)威的國家標(biāo)準(zhǔn)的制定。最后，需積累測繪人才，去除房屋、植被、水體等的影響，并完成DEM空洞區(qū)域的插值，完成DEM的后期編輯，生成可靠的全球DEM產(chǎn)品。

本文提出的面向全球測圖的精密干涉測量技術(shù)，是面向未來數(shù)字?jǐn)z影測量智能化下的自動成圖技術(shù)，精密干涉測量技術(shù)的高自動化、高精度化等特點(diǎn)，使得它能夠在無人工參與的情況下，1 min內(nèi)完成數(shù)千平方千米的DEM數(shù)據(jù)粗生產(chǎn)，并在自動化智能化篩選ICESat-GLAS點(diǎn)之后，完成大范圍數(shù)據(jù)的后處理，生成具備極高精度的全球DEM數(shù)據(jù)。這種技術(shù)將有望打破我國全球高精度DEM長期依賴SRTM數(shù)據(jù)的現(xiàn)狀，為我國全球測圖提供大范圍、高精度、高一致性的基礎(chǔ)數(shù)據(jù)。

參考文獻(xiàn)：

[1] MASSONNET D, ELACHI C. High-resolution Land Topography[J]. Comptes Rendus Geoscience, 2006, 338(14-15): 1029-1041.

[2] RODRIGUEZ E, MORRIS C S, BELZ J E, et al. An Assessment of the Srtm Topographic Products[R]. JPL D-31639, Pasadena, California: JPL, 2005.

[3] AMANTE C, EAKINS B W. Etopo1 1 Arc-minute Global Relief Model: Procedures, Data Sources and Analysis[R]. NESDIS NGDC-24, NOAA, 2009.

[4] DANIELSON J J, GESCH D B. Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010)[R]. Open-File Report 2011-1073, Reston: U.S. Department of the Interior, U.S. Geological Survey, 2011.

[5] NIKOLAKOPOULOS K G, KAMARATAKIS E K, CHRY-SOULAKIS N. SRTM Vs ASTER Elevation Products. Comparison for Two Regions in Crete, Greece[J]. International Journal of Remote Sensing, 2006, 27(21): 4819-4838.

[6] BUCKREU S. Terrasar-X/Tandem-X Mission Overview[C]∥Proceedings of TerraSAR-X/TanDEM-X Science Meeting. München: [s. n.], 2016.

[7] 國家測繪地理信息局衛(wèi)星測繪應(yīng)用中心. 國家測繪地理信息局衛(wèi)星測繪應(yīng)用中心2010-2015科技報告[R].國家測繪地理信息局衛(wèi)星測繪應(yīng)用中心，北京：2016.

Satellite Surveying and Mapping Application Center. National Administration of Surveying Mapping and Geoinformation, 2010-2015 Scientific and Technical Report[R].Beijing：Satellite Surveying and Mapping Application Center.2016.

[8] FARR T G, ROSEN P A, CARO E, et al. The Shuttle Radar Topography Mission[J]. Reviews of Geophysics, 2007, 45(2): RG2004.

[9] KRIEGER G, MOREIRA A, FIEDLER H, et al. TanDEM-X: A Satellite Formation for High-resolution SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(11): 3317-3341.

[10] BAMLER R. The SRTM Mission: A World-Wide 30 M Resolution DEM from Sar Interferometry in 11 Days[C]∥FRITSCH D, SPILLER R. Photogrammetric Week’99. Heidelberg: Wichmann, 1999.

[11] FREY O, MEIER E, NüESCH D, et al. Geometric Error Budget Analysis for Terrasar-X[C]∥Proceedings of the 5th European Conference on Synthetic Aperture Radar. Ulm, Germany: EUSAR, 2004.

[12] HUESO GONZLEZ J, WALTER ANTONY J M, BACH-MANN M, et al. Bistatic System and Baseline Calibration in TanDEM-X to Ensure the Global Digital Elevation Model Quality[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2012, 73: 3-11.

[13] SCHWERDT M, GONZALEZ J H, BACHMANN M, et al. In-orbit Calibration of the TanDEM-X System[C]∥Proceedings of 2011 IEEE International Geoscience and Remote Sensing Symposium. Vancouver, BC, Canada: IEEE, 2011: 2420-2423.

[14] DEO R, ROSSI C, EINEDER M, et al. Framework for Fusion of Ascending and Descending Pass Tandem-X Raw DEMs[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015, 8(7): 3347-3355.

[15] ROSSI C, RODRIGUEZ GONZALEZ F, FRITZ T, et al. TanDEM-X Calibrated Raw DEM Generation[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2012, 73: 12-20.

[16] 汪丙南, 向茂生. TanDEM-X系統(tǒng)相對測高性能分析[J]. 遙感學(xué)報, 2009, 13(1): 49-53.

WANG Bingnan, XIANG Maosheng. Relative Height Accuracy Analysis of TanDEM-X System[J]. Journal of Remote Sensing, 2009, 13(1): 49-53.

[17] 杜亞男, 馮光財, 李志偉, 等. Terrasar-X/Tandem-X獲取高精度數(shù)字高程模型技術(shù)研究[J]. 地球物理學(xué)報, 2015, 58(9): 3089-3102.

DU Yanan, FENG Guangcai, LI Zhiwei, et al. Generation of High Precision DEM from TerraSAR-X/TanDEM-X[J]. Chinese Journal of Geophysics, 2015, 58(9): 3089-3102.

[18] 張永勝, 黃海風(fēng), 梁甸農(nóng), 等. 星載分布式insar測高性能的理論及系統(tǒng)仿真評價方法[J]. 電子學(xué)報, 2008, 36(7): 1273-1278, 1255.

ZHANG Yongsheng, HUANG Haifeng, LIANG Diannong, et al. Theoretic and Simulation Experimental Performance Evaluation Methods of Spaceborne Distributed InSAR System[J]. Acta Electronica Sinica, 2008, 36(7): 1273-1278, 1255.

[19] 靳國旺, 張薇, 向茂生, 等. 一種機(jī)載雙天線inSAR干涉參數(shù)定標(biāo)新方法[J]. 測繪學(xué)報, 2010, 39(1): 76-81.

JIN Guowang, ZHANG Wei, XIANG Maosheng, et al. A New Calibration Algorithm of Interferometric Parameters for Dual-antenna Airborne InSAR[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(1): 76-81.

[20] 李芳芳, 胡東輝, 丁赤飚, 等. 機(jī)載雙天線inSAR對飛數(shù)據(jù)處理與分析[J]. 雷達(dá)學(xué)報, 2015, 4(1): 38-48.

LI Fangfang, HU Donghui, DING Chibiao, et al. Antiparallel Aspects of Airborne Dual-antenna InSAR Data Processing and Analysis[J]. Journal of Radars, 2015, 4(1): 38-48.

[21] 王萌萌, 黃國滿, 花奮奮, 等. 機(jī)載雙天線inSAR聯(lián)合定標(biāo)算法[J]. 測繪學(xué)報, 2014, 43(12): 1259-1265. DOI: 10.13485/j.cnki.11-2089.2014.0139.

WANG Mengmeng, HUANG Guoman, HUA Fenfen, et al. Joint Calibration Method of Airborne Dual-antenna Interferometric SAR[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(12): 1259-1265. DOI: 10.13485/j.cnki.11-2089.2014.0139.

[22] 張延冰, 郭華東, 韓春明. 利用機(jī)載雙天線InSAR數(shù)據(jù)生成高精度DEM的試驗(yàn)研究——以大面積丘陵地區(qū)為例[J]. 國土資源遙感, 2014, 26(1): 97-102.

ZHANG Yanbing, GUO Huadong, HAN Chunming. High Precision DEM Generation Using Airborne Dual-antenna InSAR Data: A Case Study of Large Hilly Areas[J]. Remote Sensing for Land & Resources, 2014, 26(1): 97-102.

[23] JEHLE M, PERLER D, SMALL D, et al. Estimation of Atmospheric Path Delays in Terrasar-X Data Using Models Vs. Measurements[J]. Sensors, 2008, 8(12): 8479-8491.

[24] GRUBER A, WESSEL B, HUBER M, et al. Operational Tandem-X DEM Calibration and First Validation Results[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2012, 73: 39-49.

[25] FRITZ T, BALSS U, BAMLER R, et al. Phase Unwrapping Correction with Dual-baseline Data for the TanDEM-X Mission[C]∥Proceedings of 2012 IEEE International Geoscience and Remote Sensing Symposium. Munich, Germany: IEEE, 2012: 5566-5569.

[26] GAO Xiaoming, LIU Yaolin, LI Tao, et al. High Precision DEM Generation Algorithm Based on Insar Multi-look Iteration[J]. Remote Sensing, 2017, 9(7): 741.