微波遙感集中編程實習(xí)報告二

./微波遙感編程實習(xí)報告實習(xí)目的本次實習(xí)希望通過編程實現(xiàn)基于R-D方程的SAR影像幾何校正,比較校正前后SAR影像特征與區(qū)別,分析距離成像幾何原理,掌握R-D校正步驟。希望通過這次實習(xí)了解SAR成像幾何原理,熟悉距離式成像過程。掌握距離-多普勒方程進(jìn)行校正的原理和難點,認(rèn)識校正所需參數(shù)意義與讀取方法。理解并掌握斜視SAR成像中距離式成像的機(jī)理與方法,了解SAR圖像特殊幾何特征成因與校正必要性；掌握基于R-D構(gòu)像方程的SAR圖像像素坐標(biāo)與坐標(biāo)間的定位關(guān)系；在此基礎(chǔ)上編程實現(xiàn)SAR圖像的幾何校正過程。實習(xí)容數(shù)據(jù)介紹SAR影像的頭文件信息中包括該幅影像上5個成像點的XYZ坐標(biāo),第一點的成像時間和5個點間的成像時間間隔,影像的第一行成像時間與成像間隔,影像距離向的成像時間與時間間隔,SAR影像4個角點的經(jīng)緯度以及投影橢球的長短半軸。這些數(shù)據(jù)用于確定任意時刻下成像點的位置矢量、速度矢量和加速度矢量,確定對應(yīng)的DEM圍以及確定DEM上各點對應(yīng)SAR影像的ij坐標(biāo)和灰度值。這些參數(shù)已經(jīng)存儲在parameter.txt里,可直接讀取使用,但注意單位轉(zhuǎn)換問題。DEM影像是ers2dem.img,可以利用DEM提供的H結(jié)合經(jīng)緯度坐標(biāo)BL,利用投影橢球信息,計算出成像點的坐標(biāo)XYZ,但注意裁剪SAR對應(yīng)到DEM的圍。SAR影像灰度信息則是存儲在DAT_01.001,注意SAR影像的數(shù)據(jù)格式,這里存儲的是復(fù)數(shù)形式,需要計算強(qiáng)度信息作為灰度值。滿足公式：強(qiáng)度2=實部2+虛部2。數(shù)據(jù)的預(yù)處理包括利用影像5個成像點的XYZ坐標(biāo)和第一點成像時間、成像時間間隔,利用多項式擬合的方法,確定軌道系數(shù)ABC,從而可求出任意時刻的衛(wèi)星位置矢量,位置矢量求導(dǎo)可得到速度矢量,再求導(dǎo)可得到加速度矢量,這些都是在之后的R-D方程中有所利用。注意為解求矩陣方便,可將第一點成像時間當(dāng)作0。其余的數(shù)據(jù)預(yù)處理包括單位轉(zhuǎn)換以及由于雙程傳輸引起的時間轉(zhuǎn)換?？紤]到DEM的分辨率是90m,而SAR影像分辨率是30m,所以還要進(jìn)行DEM加密,這里使用的方法是雙線性插方法實現(xiàn)DEM加密。R-D方程已知R-D方程由3條方程組成,包括距離方程、多普勒頻率方程和橢球方程,如下圖。圖1R-D方程利用R-D方程進(jìn)行幾何校正時,使用的是間接校正法。即利用DEM中各點的經(jīng)緯度坐標(biāo)BL結(jié)合DEM影像上的高度信息H,加上投影橢球參數(shù),計算出各點的坐標(biāo)XYZ,再由該點的位置矢量結(jié)合當(dāng)時衛(wèi)星的位置矢量、速度矢量和加速度矢量代入距離-多普勒方程,確定成像時刻,從而找到SAR影像上對應(yīng)的〔i,j,取出此處的灰度值,賦給DEM對應(yīng)的坐標(biāo)處。最后再插得到的灰度影像,即可得到糾正后的SAR影像。實習(xí)步驟下面將根據(jù)上述實習(xí)容一步步地編程實習(xí)SAR影像幾何校正。幾何過程可按概括為以下的流程圖：圖2幾何校正流程圖數(shù)據(jù)讀取和預(yù)處理需要利用的SAR影像頭文件數(shù)據(jù)已經(jīng)保存在parameter.txt中,編程實現(xiàn)數(shù)據(jù)讀取及保存。圖3參數(shù)讀取并存儲結(jié)合5點的位置矢量和第一點的成像時間和成像時間間隔,利用多項式擬合求解擬合參數(shù),由擬合參數(shù)可以求解任意時刻的衛(wèi)星位置矢量、速度矢量和加速度矢量,這里對擬合參數(shù)的求解是利用矩陣最小二乘法求解的。在矩陣計算時,由于時間t較大,高次方運算可能會超過數(shù)據(jù)圍,所以把初始成像時間設(shè)為0。注意一些數(shù)據(jù)需要進(jìn)行轉(zhuǎn)換,如投影橢球的長短半軸長單位是千米,需要轉(zhuǎn)換到米為單位；經(jīng)緯度是以度為單位的,在由經(jīng)緯度坐標(biāo)BLH計算到坐標(biāo)XYZ時需要把經(jīng)緯度轉(zhuǎn)為弧度為單位；考慮到雙程傳輸,時間需要進(jìn)行改正。需要利用到DEM數(shù)據(jù)和SAR影像,利用GDAL庫的RasterIO函數(shù)將影像數(shù)據(jù)讀入。注意SAR影像DAT_01.001是以復(fù)數(shù)形式存儲,分為實部和虛部,數(shù)據(jù)格式是CInt16,而我們需要的是它的灰度信息,灰度信息是影像的強(qiáng)度,符合公式：強(qiáng)度2=實部2+虛部2,利用公式求出其強(qiáng)度信息并讀出來存儲起來。編程實現(xiàn)DEM影像讀取和SAR影像讀取并存儲。圖4DEM的讀取圖5SAR的讀取由R-D方程計算行列號利用R-D方程間接校正,即利用DEM的坐標(biāo),結(jié)合衛(wèi)星位置矢量、速度矢量、加速度矢量代入R-D方程,求解出該點對應(yīng)到SAR影像上的行列號。對多普勒方程求微分,有式中表示衛(wèi)星速度矢量和位置矢量,表示DEM坐標(biāo),對DEM的任一點來說,成像時間未知,多普勒方程不為0,可針對方位向時間進(jìn)行迭代。以中間行成像時間為初始迭代時間,求出當(dāng)前時刻的,結(jié)合DEM上該點的,代入上式可以求出,令,再進(jìn)行下一次迭代,隨著迭代次數(shù)的增加,會減小,成像時間會越趨于真實值?？梢酝ㄟ^設(shè)置一個閾值,當(dāng)時,迭代停止,此時求出的時間為真實方位向時間。為使循環(huán)不陷入死循環(huán),還需設(shè)置一個迭代次數(shù)限制,所以最后迭代次數(shù)到了或是滿足小于閾值的條件時,迭代停止,求出真實方位向成像時間。再利用第一行成像時間和方位向成像時間間隔,求出行號。這里設(shè)置閾值且,公式為求出方位向時間后,可求出正確的衛(wèi)星位置矢量,結(jié)合求出斜距,距離向時間為。再由距離向時間求出列數(shù)最后注意判定行列號在SAR影像包括的圍時,給該點賦上灰度值,若不在區(qū)域時,把灰度值賦為0。行列號計算編程如下：圖6R-D方程間接幾何校正結(jié)果圖像的存儲與輸出利用GDAL庫實現(xiàn)上一步中保存下來的灰度矩陣存為圖像并輸出,可以將校正前后的SAR影像對比。代碼部分如下：圖7結(jié)果的存儲與輸出圖8運行結(jié)果實習(xí)心得通過這次幾何校正編程實習(xí),我弄清楚了利用R-D方程間接幾何校正的具體流程,需要結(jié)合DEM數(shù)據(jù),迭代求出SAR影像上的行列號。并且了解到GDAL庫如何讀取影像已經(jīng)存儲影像,GDAL庫還可用于影像截取指定區(qū)域?qū)崿F(xiàn)影像插,插方式可以選擇最近鄰插和雙線性插等方式。在這次實習(xí)中,我也遇到了許多問題。首先利用GDAL中的RasterIO可以讀取出影像各個像元的灰度值,對DEM來說其灰度值就是它的高程,還可以讀取出影像左上角點的地理坐標(biāo),可以利用這個求出每個像元的地理坐標(biāo),這里是經(jīng)緯度值。我原本希望將計算出的經(jīng)緯度值都以矩陣的方式存儲起來,但程序運行到這里時總是提示錯誤,錯誤理由檢查后發(fā)現(xiàn)是存沒分配到,在詢問老師后,老師說每個程序運行時存都是一定的,我這樣寫存不足沒法實現(xiàn)分配存儲,建議我在每次需要利用BLH求解坐標(biāo)XYZ時再計算經(jīng)緯度BL,也不要把它存儲起來。然后就是計算對應(yīng)SAR影像的行列號處,每次運行完我編寫的程序后得到的影像都是一片黑,設(shè)斷點調(diào)試程序時發(fā)現(xiàn)每次的SAR影像上對應(yīng)的行號i計算出來的都是負(fù)值,不在SAR影像包括的圍。經(jīng)檢查后發(fā)現(xiàn)是自己的經(jīng)緯度BL弄反了,并且由經(jīng)緯度坐標(biāo)BLH計算地理坐標(biāo)XYZ的公式就寫錯了,導(dǎo)致每次運行結(jié)果都是負(fù)值,修改了這兩處錯誤后,再運行編寫的程序即可得到SAR影像幾何校正后的結(jié)果。對比糾正前后的SAR影像,發(fā)現(xiàn)不僅是影像的大小、形狀有變化,主要差別在于幾何校正后的SAR影像包含了地理坐標(biāo)信息,,幾何校正后的影像反映的才是雷達(dá)成像時雷達(dá)波掃過地表的真實區(qū)域圍和回波分布情況?？偟膩碚f,這次編程實習(xí),我收獲很大。首先是真正動手編程去實習(xí)基于R-D方程的幾何糾正,這讓我對課本知識有了更深刻地理解,而