專利名稱：基于線陣推掃式異步采樣衛星影像幾何模型的正射校正方法
技術領域：
本發明涉及一種正射校正方法，特別是一種基于雙中心線投影的線陣推掃式異步采樣衛星影像正射校正方法，屬于衛星影像處理領域。
背景技術：
從上世紀九十年代末至今，由于高分辨率商業遙感衛星影像的大量獲取和分發，人們開始將這些衛星影像用于地理空間數據的采集和建庫。衛星遙感影像不可避免地存在各種變形，地理空間數據庫要求消除變形，得到精確的二維或者三維地理空間信息。高分辨率光學遙感衛星大都采用線陣推掃式傳感器，獲取的每一行影像有獨立的外方位元素，不能采用傳統的中心投影共線方程進行整體幾何處理。因此，對線陣推掃式衛星影像的幾何建模一直是學術界研究的一個重要課題。
線陣推掃式遙感衛星的成像模式主要有兩種同步采樣成像模式和異步采樣成像模式。在同步采樣模式模式下，傳感器的地面采樣速度與衛星的地速相等，在地圖投影坐標系中傳感器的外方位元素盡量保持不變。在異步采樣模式下，傳感器的地面采樣速度與衛星的地速不相等。
針對線陣推掃式遙感衛星不同的成像模式，學術界主要提出了嚴密幾何模型、簡化幾何模型和主觀模型。
嚴密幾何模型根據成像光束在傳感器線陣方向符合中心投影原理，對每一行影像使用中心投影共線方程 其中
在嚴密模型中不同行之間的外方位元素的變化采用一般多項式進行描述。嚴密模型需要傳感器的物理和幾何參數以及衛星星歷數據，并且存在參數之間的高度相關的現象，導致參數解不穩定。
在嚴密模型的基礎上假設傳感器在成像過程中姿態不變，衛星沿直線軌道運行，或者舍棄一些在具體條件下不重要的參數，得到簡化的幾何模型，其中包括平行推掃式成像幾何模型和平行投影模型。簡化模型只能描述衛星影像的主要變形因素，在地面起伏較大以及衛星姿態不穩定的條件下，簡化么模型得到的幾何精度較低。
主觀模型不考慮成像的幾何原理和傳感器的物理參數，直接采用一般多項式或者有理多項式來建立影像坐標與地面坐標的關系，其參數與實際的成像幾何沒有直接的對應關系。采用一般多項式作為模型時，與確定的地面控制點擬合很好，但在其它點的內插值則可能有明顯的偏離，即可能在某些點處產生振蕩，導致用多項式近似計算的中誤差明顯超出平均誤差。有理多項式模型也是一種主觀模型，從理論上講有理多項式模型所描述的數學關系只在控制點處成立，而在其它地方都是近似的。因此，有理多項式模型的精度與控制點的精度、分布、數量以及糾正范圍密切相關。
由于線陣推掃式異步采樣遙感衛星出現得比較晚，其成像幾何尚未引起足夠的重視，1B級異步采樣影像的幾何模型還未建立。
發明內容
本發明目的就是針對線陣推掃式異步采樣衛星的1B級影像的幾何建模。
實現本發明目的采用的技術方案是利用雙中心線投影模型描述線陣推掃式異步采樣衛星1B級影像的影像坐標與對應點地面坐標之間的數學關系。線陣推掃式異步采樣衛星影像正射糾正方法包括以下步驟 (1)雙中心投影模型參數初始值的確定； (2)誤差方程系數的計算與模型參數改正數的解算將誤差方程系數矩陣標準化，將常數項中心化，用標準化的誤差方程解算標準化的模型參數改正值，再采用標準化的逆過程求出模型參數改正值； (3)輸出正射影像范圍的計算； (4)正射影像格網點灰度值的內插計算。
本發明通過下式建立1B級影像的影像坐標與對應的地面坐標之間的關系 上式1B級影像雙中心線投影共線方程，該方程為隱函數方程組，除了像點觀測值x和y以外，還有xo、yo、Xo、Yo、Zo、
ω、κ、
ω′、κ′、f、f′、α共14個參數，只要有7個以上的像點坐標和對應點的地面坐標，就可以用最小二乘法唯一確定上述參數的解。
本發明具有以下優點 雙中心線投影模型不但適用于處理經過幾何改正的線陣推掃式異步采樣衛星影像，而且可以通過用來對中心投影影像、線陣推掃式同步采樣的原始衛星影像進行幾何處理。
圖1為本發明方法的流程圖。
圖2為異步采樣對稱位置成像光束交會示意圖。
圖3為線陣方向與運行軌道垂直時成像光束的交會示意圖。
具體實施例方式 中心投影是廣為熟悉的一種投影模式，實現中心投影必須具備2個要素投影中心、投影面。從物點發出的光線通過投影中心，在投影面上形成影像，從而實現中心投影。在中心投影中，投影中心是幾何上的一個點，光學上以一個小孔替代小孔。設想該幾何點為兩條直線的交點，光學上可以用兩個相交的狹縫來實現，投影面不變，中心投影仍然成立。將上述假設中的兩個相交的狹縫分離，使它們既不平行，也不相交，同時保持投影面不變，從物點發出的光線依次通過兩個狹縫，在投影面上仍然可以成像。這里把這種由物點發出的光線依次通過兩個既不平行、也不相交的狹縫，然后在投影面上形成影像的過程稱作雙中心線投影。幾何上可以把這兩個既不平行、也不相交的狹縫看作兩條既不平行、也不相交的直線，這里把它們稱做投影中心線。由此可見，雙中心線投影包含3個要素2條投影中心線和1個投影面。
線陣推掃式傳感器異步采樣的雙中心線投影原理 探討異步采樣衛星遙感影像成像的幾何原理，設想如圖3所示，衛星首先從S1點開始，對某目標區域進行掃描，此時衛星傳感器投影中心是S1，主點為o1，主光軸S1o1向前傾斜一個角度φ，傳感器線陣列為a1b1，傳感器的焦距為f。衛星在向前運行的過程中，對目標區域進行掃描，同時為了提高傳感器的感應時間，衛星以一個接近均勻的角速度向后旋轉，到達位置S的時候，衛星傳感器線陣列為ab，此時主光軸So垂直于S1S。衛星繼續向前運行，到達位置S2時，傳感器線陣列為a2b2，主光軸S2o2向后傾斜φ。假設衛星平臺向后轉的速度是是對稱的，衛星從S1運行到S所用的時間和從S運行到S2所用的時間相等，又因為衛星運行的速度也是一個常數，所以有 S1S＝S2S(3) 這里繼續采用一個基本假設，即在一景影像范圍內，衛星平臺質心的運行軌道近似為一條直線。由于衛星在進行異步采樣時，要不停的向后轉動，因此，在采樣過程中，傳感器的投影中心相對于衛星質心來說，會有一定的變化。衛星本身的尺寸較小，與掃描一景影像的過程中衛星的運行距離相比，傳感器投影中心相對于衛星質心位置變化是一個很小的量，所以，在一景影像范圍內，傳感器投影中心的運行軌跡S1S2也可以看作近似為一條直線。
衛星平臺在采集一景影像的過程中，由于運行時間較短，可以認為衛星平臺的側偏角和自旋角不發生變化，也就是說，衛星傳感器在對目標區域的掃描過程中，傳感器線陣方向不發生變化，而是保持相互平行的狀態。
在上述兩個基本前提條件下，以圖3中的S為原點，以oS方向為w軸正方向，以垂直于面Sab且與SS2成銳角的方向為u軸正方向，建立起右手直角坐標系S-uvw作為傳感器參考坐標系。設S1在S-uvw坐標系中的坐標為S1(-Δu，-Δv，0)，因為S1、S2相對于S是對稱的，則有S2在S-uvw坐標系中的坐標為S2(Δu，Δv，0)，因此，傳感器在S1、S、S2三點的外方位元素分別為
將上述參數代入式(1)的第一個等式和(2)分別得到 將上述第三式與第一式聯立求解，得到 將(7)的第三式與第二式聯立求解，也得到(8)。
從式(8)可以看出，它實際上是一個直線方程，該直線在u軸方向和w軸方向的分量為零。由此證明了Sab、S1a1b1、S2a2b2三個面交于同一條直線，該直線為平行于ab，與S的距離取決于SS1和SS2在垂直于Sab的直線方向的分量以及φ的大小。
由以上的分析可以看出，衛星在掃描成像的過程中，投影中心處于點S的時候，以及處于與S成對稱位置的S1點和S2的時候，成像光束交會于一條平行于傳感器線陣的直線。進一步可以證明，當衛星運行軌道與傳感器線陣相互垂直時，經過主點的成像光束交會于一個點，與主光束對稱的成像光束也分別交會于各自的點，這些點構成一條平行于傳感器線陣的直線，如圖3所示，So、So1、So2交會于點S′，Sa1、Sa2交會于點S1′，Sb1、Sb2交會于點S2′。
從圖3可以看出，在掃描一景影像的過程中，所有的成像光束必須經過直線S1S2，S1S2的作用類似于中心投影中的投影中心，這里把它稱作第一投影中心線。為了探討異步采樣衛星遙感影像的整體構像方程，這里再進一步假設衛星從S1開始，掃描到主光軸與衛星運行軌跡相垂直的位置S，然后繼續掃描，直到與S1對稱的位置S2，完成一景影像的掃描，整個過程中的成像光束全都交會于一條位于Sab平面內，且平行于ab的直線，也就是說，整景影像的成像光束除了經過S1S2直線外，還同時經過另外一條直線。由于該條線的作用也類似于中心投影中的投影中心，這里把它稱作第二投影中心線。從幾何原理可知，過一個已知空間點和兩條既不平行也不相交的空間直線，可以唯一地確定一條空間直線。因為過地面點且同時經過兩條投影中心線的直線與過像點且同時經過兩條投影中心線的直線是同一成像光線，它們應該共線，可形成共線方程。由于上述結論是從成像光束的空間關系得出的，因此不會因為影像空間參照系的改變而改變。設想已知一定數量的點的像點坐標和像點對應的地面點的地面坐標，就可以用空間后方交會法，求出這兩條直線的方程，從而確立從像點坐標到地面坐標的對應關系，實現異步采樣衛星遙感影像的幾何處理。
本發明利用雙中心線投影模型進行正射糾正的步驟 (1)雙中心線投影模型參數初始值的確定 基于雙中心線投影模型的EROS A1衛星的1B級影像的模型參數的解算步驟如下 根據雙中心線模型的定義，參數中f實際上是以影像比例尺表示的傳感器焦距，由于EROS A1衛星的1B級影像是根據衛星星歷和傳感器姿態參數，將原始影像投影到WGS84 UTM坐標系后得到的，在處理EROS A1衛星的1B級影像時，可以用影像中心對應的地面點與衛星之間的距離作為f的初始值。而f′是傳感器線陣與傳感器投影中心組成的平面在不同位置的交線與虛擬像平面之間的距離，可以根據衛星地面采樣弧長、衛星軌道弧長在地面的投影長度以及f等3個參數按下式計算 式中，D為衛星的飛行弧長，d為地面采樣距離，f為第一焦距。上述參數可以從影像的元數據中得到或者根據影像元數據推算出來。
由于EROS A1的1B級影像是經過幾何糾正的，并且投影到WGS84 UTM坐標系，而正射糾正也在WGS84 UTM坐標系中進行，所以主點o的坐標可以直接在EROS A1的1B級影像上量取。具體做法是找出影像掃描行方向長度最短的一條線的中點，量出該點的X、Y坐標，作為該景影像的主點坐標Xo和Yo的初始值。而影像掃描行方向可以由EROS A1的1B級影像的上邊緣方向確定。
轉角
是衛星處于軌道平面內軌道的垂線方向與WGS84 UTM坐標系的Z軸的夾角，實際上反映了衛星軌道與水平面的夾角，一般情況下近似為零。轉角ω是衛星成像時主光軸與本地大地法線之間的夾角，可以用元數據中傳感器的側偏角作為初始值。轉角κ是衛星軌道方向在WGS84 UTM坐標系中的方位角，可以由EROS A1的1B級影像的上邊緣的向下的垂線方向與WGS84 UTM坐標系的X軸之間的夾角來確定。而轉角α是衛星實際飛行軌跡與衛星軌道方向的夾角，是地球自轉和傳感器線陣方向與衛星運行方向不垂直兩個因素共同造成的，其初始值可以根據衛星的軌道參數、地球自轉速度以及影像所處的地理位置來計算；也可以在影像上直接量取，其做法是在EROS A1衛星的1B級影像上畫出第一行影像和最后一行影像的中點之間的連線，量取該連線與第一行影像的垂線之間的夾角β，然后利用下面的近似公式計算α (2)誤差方程系數的計算與模型參數改正數的解算 誤差方程的系數是根據上一個步驟計算出的模型參數初始值、量測的像點坐標、對應于像點的地面點坐標等三組數據計算的。由于EROS A1的1B級影像經過了幾何改正，并且投影到WGS84 UTM坐標系中，所以像點坐標與地面坐標同屬于一個坐標系，像點坐標也用WGS84 UTM坐標表示。將模型參數初始值、像點坐標、對應于像點的地面點坐標代入雙中心線投影共線方程，就可以得到誤差方程的系數矩陣和常數項。設一景影像中量測的控制點有n個，就可以得到2n個誤差方程，當n＞4時，需要用最小二乘法求解模型參數的改正值。
由于模型中各個參數值在數量級上差別相當大，直接將誤差方程法化后進行平差計算會遇到法方程病態化問題。為了使法方程各個元素在同一個數量級，需要將誤差方程系數矩陣標準化，將常數項中心化，用標準化的誤差方程解算標準化的模型參數改正值，再采用標準化的逆過程求出模型參數改正值，具體步驟如下。
將上面求出誤差方程系數矩陣記為 常數項記為 W＝[w1 w2 … w2n]′ (12) 求出誤差方程系數矩陣每一列的平均值 引進每一列的縮放系數 然后將矩陣A的每一列元素減去該列元素的平均值，再除以該列的縮放系數Rj得到 上式即為標準化的誤差方程系數矩陣，用標準化的誤差方程計算出了法方程系數矩陣實際上就是模型之間的相關系數矩陣。同樣將誤差方程常數項中心化，記 W＝[w1-w w2-w…w2n-w]′ (17) 標準化的誤差方程系數矩陣和中心化的誤差方程常數項構成新的誤差方程 V＝A·δX+W (18) 用最小二乘法解上式，得到經過縮放后的參數改正量。
δX＝[A′·P·A]-1A′·P·W (19) 參數的實際改正量為 δxj＝δxj×Rj(20) 在式(19)中，P為觀測值權陣。如果像控點的地面坐標測量誤差對空間后方交會平差精度有比較顯著的影響，可以對地面點坐標列立誤差方程，給定其權矩陣，采用上面一樣的方法解算。
解算完畢后，將各項參數改正量與對應的參數初始值相加，得到改正后的參數值。將改正后的參數值、像點坐標以及像點對應的地面點坐標代入雙中心線投影共線方程，重新計算誤差方程系數項和常數項，按照前面的過程重新平差計算，這樣反復迭代，直到各項參數的改正量小于允許值為止。
平差結束后，用各項參數最后的值和像點坐標、地面坐標重新計算出誤差方程的常數項，得到每個像點坐標的殘差。利用像點坐標的殘差統計像點坐標測量精度，作為精度評定的依據。在控制點較多的情況下，可以采用一部分點作為控制點，另外一部分點作為獨立檢查點，利用檢查點的地面坐標和模型參數平差值，求出這些檢查點的像點坐標，將實際量測的檢查點的像點坐標與計算出的像點坐標進行比較，以評定模型定向的精度。
(3)輸出正射影像范圍的計算 糾正后圖像邊界范圍的確定過程如下 3-1求解原始圖像的四個角點1、2、3、4對應的地面點在地圖投影坐標系中的坐標值，得到8個坐標值(X1，Y1)、(X2，Y2)、(X3，Y3)、(X4，Y4)。
由于雙中心線投影模型從影像的二維坐標到地面的三維坐標變換需要用到地面點的高程，可以用EROS A1 1B級影像四個角點在WGS84 UTM坐標系中的坐標值，直接從DEM中內插出對應的地面的高程值。將角點的影像坐標、對應的高程值和解算出的模型參數值代入共線方程，求解出每一個角點對應的地面點的坐標，此時算出的地面點坐標與從影像上直接量取的坐標值不會相同。用解算出的地面點坐標再次從DEM中內插出對應的高程值，然后將角點在影像上的坐標值、模型參數和內插出的地面點高程再次代入共線方程，再次求出角點對應的地面點坐標。采用上述方法反復迭代，直到兩次計算出的地面坐標相等，即得出影像四個角點對應的地面點坐標。
3-2利用8個坐標值計算正射影像的范圍 式中Xmin，Ymin為數字正射影像的起始點(左下角)坐標，Xmax，Ymax為數字正射影像的終止點(右上角)坐標，d為數字正射影像的地面分辨率。
正射影像格網點對應于原影像中像點坐標的計算 求出數字正射影像的范圍后，即可計算出正射影像上每一個像元中心對應的地面坐標值。以正射影像的左下角開始，行坐標向上取正值，列坐標向左取正值，對應第i行第j列像元中心的坐標值為 利用上式計算出Xij和Yij坐標后，就可根據DEM內插出對應的高程值Zij。將數字正射影像每個格網點的三維坐標(Xij，Yij，Zij)代入共線方程，即可求出對應的像點坐標。
(4)正射影像格網點灰度值的內插計算 求出數字正射影像格網點中心在原影像中對應的像點坐標后，即可內插出該格網點的影像灰度值。影像灰度值的內插方法可以采用SINC函數法、雙三次卷積重采樣法和最鄰近像元法等。
如果計算出的影像坐標超出原影像的范圍，可采用空白值填充。采用上述方法將正射影像上每一個格網點的灰度值計算出來，并賦給該像元，即得到完整的數字正射影像，整個正射糾正過程結束。
權利要求
1.一種基于線陣推掃式異步采樣衛星影像幾何模型的正射校正方法，其特征在于采用雙中心線投影的投影模式來描述線陣推掃式異步采樣衛星影像的幾何變形規律，通過恢復雙中心線投影的各項參數，建立衛星影像像素坐標與對應的地面坐標之間的關系，從而實現對線陣推掃式異步采樣衛星影像正射校正，包括以下步驟
(1)確定雙中心線投影模型參數初始值；
(2)誤差方程系數的計算與模型參數改正數的解算；
(3)計算輸出正射影像的范圍；
(4)正射影像格網點灰度值的內插計算。
2.根據權利要求1所述基于線陣推掃式異步采樣衛星影像幾何模型的正射校正方法，其特征在于通過下式建立1B級影像的影像坐標與對應的地面坐標之間的關系
上式1B級影像雙中心線投影共線方程，該方程為隱函數方程組，除了像點觀測值x和y以外，還有xo、yo、Xo、Yo、Zo、
ω、κ
ω′、κ′、f、f′、α共14個參數，只要有7個以上的像點坐標和對應點的地面坐標，就可以用最小二乘法唯一確定上述參數的解。
全文摘要
本發明公開了一種基于線陣推掃式異步采樣衛星影像幾何模型的正射校正方法，采用雙中心線投影的投影模式來描述線陣推掃式異步采樣衛星影像的幾何變形規律，通過恢復雙中心線投影的各項參數，建立衛星影像像素坐標與對應的地面坐標之間的關系，從而實現對線陣推掃式異步采樣衛星影像正射校正，包括以下步驟(1)確定雙中心線投影模型參數初始值；(2)誤差方程系數的計算與模型參數改正數的解算；(3)計算輸出正射影像的范圍；(4)正射影像格網點灰度值的內插計算。
文檔編號G06T7/00GK101609551SQ20091006332
公開日2009年12月23日 申請日期2009年7月24日 優先權日2009年7月24日
發明者龔健雅, 胡興樹, 眭海剛, 馬國銳 申請人:武漢大學