一種預先確定星載激光測高儀的足印位置的方法與流程
本發明涉及星載激光測高儀在軌幾何檢校領域，特別涉及一種預先確定星載激光測高儀的足印位置的方法，應用于星載激光測高儀在軌幾何檢校過程中的激光地面探測器布設位置的預報。
背景技術：
在軌幾何檢校試驗的成功開展是星載激光測高儀獲得高精度測高數據的必要條件之一。星載激光測高儀足印位置的預報是在軌幾何檢校的第一步，也是極為關鍵的一步。預報的精度不僅關系著整個檢校試驗外業布設工作的難易程度，甚至決定著試驗的成敗。
目前，僅有美國ICESat-GLAS(Ice Cloud Elevation Satellite Geoscience Laser Altimeter System)開展過相關的研究和工程實踐，其激光脈沖發射重頻為40Hz，激光沿軌相鄰兩個足印中心間距為170米左右，地面足印平均大小為65米左右，沿軌布設大于170米的探測器陣列，必然可以捕獲到光斑。由于其硬件的優勢，使其在足印預報方法上并不需要開展過多的研究。但是其高重頻是根據衛星的使用目的而設計的，即對冰川的厚度進行高精度測量；而對于光學測繪衛星而言，星上本身攜帶了高分辨率相機、姿態測控系統和軌道確定系統等多個載荷，激光測高儀的脈沖出射頻率必定要控制在一定的范圍內，以保證整星的能源分配問題，且用于符合光學測繪的激光測高數據在實際應用中也不需要過高的激光出射頻率。低頻率的脈沖出射為星載激光測高儀的在軌幾何檢校增加了試驗難度，最突出地體現在激光足印的預報上。由于國際上暫時沒有涉及該方面技術的測繪衛星在軌運行，激光足印預報方面的研究幾乎屬于空白。
因此，本領域需要一種能夠準確預報星載激光測高儀足印位置的方法。
技術實現要素：
為此，本發明提出一種預先確定星載激光測高儀的足印位置的方法，彌補了目前激光測高儀高精度足印位置預報方法的空白，并解決了現有簡單預報技術的限制和缺陷導致的一個或多個問題。
本發明另外的優點、目的和特性，一部分將在下面的說明書中得到闡明，而另一部分對于本領域的普通技術人員通過對下面的說明的考察將是明顯的或從本發明的實施中學到。通過在文字的說明書和權利要求書及附圖中特別地指出的結構可實現和獲得本發明目的和優點。
本發明提供了一種預先確定星載激光測高儀的足印位置的方法，其特征在于，所述方法包括以下步驟：
步驟1，計算激光測高儀在衛星本體坐標系中的指向Mb_l，步驟1具體包括以下子步驟：
步驟1.1，預估激光測高儀與平臺之間的安裝角(αr,βr)；
步驟1.2，基于步驟1.1預估的安裝角，計算激光測高儀在衛星本體坐標系中的指向Mb_l；
步驟2，將激光測高儀的指向從衛星本體坐標系旋轉至軌道坐標系Mo_l；
步驟3，將激光測高儀的指向從軌道坐標系旋轉至地球固定坐標系，獲得激光測高儀在地球固定坐標系下的指向MT_l；
步驟4，通過計算地球固定坐標系下的激光測高儀的指向與地球橢球的交點，確定星載激光測高儀的足印位置，其中，所述步驟4具體包括以下子步驟：
步驟4.1，根據下述公式計算激光指向與地球橢球相交的地球固定坐標系下的地面點坐標(X1,Y1,Z1)：
其中，MT＝[X1 Y1 Z1]T，X1、Y1和Z1為待求的激光指向與地球橢球相交的地球固定坐標系下的地面點坐標，μ為待求的比例系數，MT_l為激光測高儀在地球固定坐標系下的指向，P0為衛星激光質心在地球固定坐標系的坐標值，a和b分別為地球橢球長半軸和短半軸長度，h為地面大地高程；
步驟4.2，根據下述公式將地球固定坐標系下的地面點坐標(X1,Y1,Z1)轉換成大地經緯度坐標系下的地面點坐標(L,B,H)：
其中，e1和e2分別為第一偏心率和第二偏心率，a和b分別為地球橢球長半軸和短半軸長度，X1、Y1、Z1為通過步驟4.1計算出的激光指向與地球橢球相交的地球固定坐標系下的地面點坐標。
本發明還提供了一種計算機程序，其中，當執行所述計算機程序時，執行本發明任一實施例所述的方法。
本發明針對星載激光測高儀足印位置預報領域的空白，提供了一種預先確定星載激光足印位置的方法，該方法結合了傳統的光學相機成像的原理，并充分顧及了激光測高儀自身的測高機理，預報結果能夠較好地符合衛星真實的飛行過境情況，具有較高的預報精度，并且極大限度地減少不必要的人力物力，節省了大量的試驗開銷。排除激光指向角、衛星平臺姿態預報數據和衛星質心軌道預報數據精度損失的情況下，該方法的星載激光測高儀足印位置預報精度能夠控制在5米。
附圖說明
圖1為根據本發明實施例的確定星載激光測高儀的足印位置的方法的流程圖。
具體實施方式
下面參照附圖對本發明進行更全面的描述，其中說明本發明的示例性
實施例。
根據本發明的實施例，可以隨機抽取星載激光測高儀進行足印位置預報模型的構建；優選的，本發明的實施案例以資源三號02星的激光測高儀的足印位置預報模型構建為例，但并不限于此。
如圖1所示，激光測高儀預報模型的構建和光學嚴密幾何成像模型的構建存在同樣的過程，涉及到各個坐標系之間的旋轉，通過各坐標系的旋轉最終將星上指向傳輸至地面點，本發明所提出的星載激光測高儀足印位置預報模型構建方法具體包括以下步驟：
步驟1，計算激光測高儀在衛星本體坐標系中的指向Mb_l，步驟1具體包括以下子步驟：
步驟1.1，預估激光測高儀與平臺之間的安裝角(αr,βr)；
由于衛星發射升空過程中的劇烈震蕩和在軌后太空環境與地面實驗室的差異較大，會引起激光測高儀與平臺之間的安裝關系發生變化，導致激光測高儀的安裝角發生變化，因此需要計算真實的激光測高儀與平臺之間安裝角。
參考資源三號02星激光測高儀設計指標，獲取激光測高儀的實驗室標定的安裝角α＝90.86229°，β＝0.861985°。以實驗室標定安裝角為基礎，采用基于地形匹配的安裝角預估方法計算真實的指向角。這里采用金字塔計算策略，逐步縮小計算區域和遍歷間隔，將時間限定在一定的區間內，保證計算的精度。將金字塔分為3層來進行計算，第一層金字塔的區域間隔為2°，遍歷間隔為0.1°，第二層金字塔的區域間隔為0.4°，遍歷間隔為1′，第三層金字塔的區域間隔為0.02°，遍歷間隔為1″，具體的計算公式如下式所示：
其中，(αr,βr)為預估的安裝角(其中，αr為激光指向與本體坐標系X軸之間的夾角，βr為激光指向與本體坐標系Z軸之間的夾角)，表示從1層至3層逐層計算，αi、βi表示第i層遍歷的安裝角，表示括號內第i層取最小值時對應的αi、βi，并將其賦值給α0、β0，當i＝1首次計算時，α0、β0的值取實驗室測量值90.86229°、0.861985°，q(i)為第i層遍歷間隔，p(i)為第i層的區域間隔，表示對括號內進行遍歷，行遍歷序號m，列遍歷序號n，E(αi,βi)為αi、βi指向下計算某軌的高程值與已知高程值之間的標準差，具體的高程值計算公式可參照激光測高儀的嚴格幾何測高模型，已知高程采用公開全球范圍的AW3D30數據。
以資源三號02星的第382軌激光數據為例，計算過程如下：首先進行第一層金字塔遍歷，區域間隔為2°，遍歷間隔為0.1°，
經過第一層的遍歷，可以獲得第一層的預估結果α0、β0分別為89.9°、0.0°，下面進行第二層金字塔的遍歷，區域間隔為0.4°，遍歷間隔為1′，
經過第二層的遍歷，可以獲得第二層的預估結果α0、β0分別為89.945°、0.046°，下面進行第三層金字塔的遍歷，區域間隔為0.02°，遍歷間隔為1″，
通過上式可以求得預估的安裝角αr＝89.94445°，βr＝0.04692°
步驟1.2，基于步驟1.1預估的安裝角，根據下述公式計算激光測高儀在衛星本體坐標系中的指向Mb_l：
式中，Mb_l為步驟1待求的激光測高儀在衛星本體坐標系中的指向，αr、βr為步驟1.1求得的激光測高儀的安裝角。
步驟1.1求得的激光測高儀安裝角預估值代入到上式可以求得預估指向
步驟2，將激光測高儀的指向從衛星本體坐標系旋轉至軌道坐標系Mo_l。
根據光束的物理傳播路徑，傳輸的第一步是將激光測高儀的指向從衛星本體坐標系旋轉至軌道坐標系，步驟2具體包括以下子步驟：
步驟2.1，估計衛星平臺的本體坐標系相對于軌道坐標系的姿態值數值單位為弧度，其中，roll、pitch、yaw分別為衛星本體坐標系相對于軌道坐標系的三軸歐拉角：橫滾角、俯仰角和偏航角，即，roll為衛星本體坐標系相對于軌道坐標系的橫滾角，pitch為衛星本體坐標系相對于軌道坐標系的俯仰角，yaw為衛星本體坐標系相對于軌道坐標系的偏航角。
一般情況下，由于衛星姿控系統的高精度實時控制，平臺的本體坐標系相對于軌道坐標系的姿態變化較小，因此這里采用了取零值的方法，即三軸姿態均為零
步驟2.2，根據下述公式計算本體坐標系相對于軌道坐標系的旋轉矩陣Mo_b；
a1＝cos(pitch)*cos(yaw)
a2＝-cos(pitch)*sin(yaw)
a3＝sin(pitch)
b1＝-sin(roll)*sin(pitch)*cos(yaw)+cos(roll)*sin(yaw)
b2＝sin(roll)*sin(pitch)*sin(yaw)+cos(roll)*cos(yaw)
b3＝sin(roll)*cos(pitch)
c1＝-cos(roll)*sin(pitch)*cos(yaw)-sin(roll)*sin(yaw)
c2＝cos(roll)*sin(pitch)*sin(yaw)-sin(roll)*cos(yaw)
c3＝cos(roll)*cos(pitch)
將步驟2.1得到的衛星平臺的本體坐標系相對于軌道坐標系的姿態值(即，)代入到Mo_b，可得本體坐標系相對于軌道坐標系的旋轉矩陣
本實施例的步驟2.1優選的采用了的方法，即代入到步驟2.2計算。
步驟2.3，根據下述公式將激光測高儀的指向由本體坐標系轉換到軌道坐標系；
Mo_l＝Mo_b*Mb_l
其中，Mo_l為待求的激光測高儀在軌道坐標系下的指向，Mo_b為步驟2.2計算的本體坐標系相對于軌道坐標系的旋轉矩陣，Mb_l為步驟1.2計算的激光測高儀的指向。
將步驟1.2和步驟2.2計算得到的激光測高儀載荷與平臺之間的預估指向與本體坐標系相對于軌道坐標系的旋轉矩陣相乘，可以得到激光測高儀的指向由本體坐標系轉換到軌道坐標系的向量
步驟3，將激光測高儀指向從軌道坐標系旋轉至地球固定坐標系。
根據光束的物理傳播路徑，傳輸的第二步是將激光測高儀指向從軌道坐標系旋轉至地球固定坐標系，步驟3具體包括以下子步驟：
步驟3.1，計算衛星激光質心在地球固定坐標系的坐標值P0＝[X0 Y0 Z0]T和位移速度值V0＝[VX0 VY0 VZ0]T；
根據本發明的優選實施例，步驟3.1的具體計算過程如下：
式中，為指定時刻的慣性系衛星位置關于時間的二階導數即加速度，為待求量t為某一指定時刻，r為指定時刻的慣性系衛星位置，v為指定時刻的慣性系衛星速度，CD為大氣阻力系數，CR為太陽光壓系數，aiR為作用在衛星軌道徑向上的第i段經驗加速度，aiT為作用在衛星軌道切向上的第i段經驗加速度，aiN為作用在衛星軌道法向上的第i段經驗加速度，eR(t)為指定時刻衛星軌道徑向方向單位向量，eT(t)為指定時刻衛星軌道切向方向單位向量，eN(t)為指定時刻衛星軌道法向方向單位向量，ti為第i段經驗加速度參數開始時刻，ti+1為第i+1段經驗加速度開始時刻，也即第i段經驗加速度結束時刻。
參考時刻衛星位置速度、大氣阻力系數CD、太陽光壓系數CR、經驗加速度參數aiR、aiT、aiN等來自于預報區間之前的精密軌道擬合。
通過上式可求得新的加速度因此可以利用下式求待估參數(P0,V0)，已知初值位置和速度P0'，V0'，CD，CR及加速度，即
式中，CD和CR作為常值估計。
根據待預報地區，粗略估計衛星過境時間t，利用軌道預報計算公式，結合衛星過境前24小時的長時段軌道數據作為初始值，對過境的軌道參數進行預報。由于該預報一個漸進遞推的長時段處理，過程涉及的數據較多，這里僅給出一例作為參考，24小時前的軌道數據初值為：P0'＝[-728959.1 6827631.7 473238.5]T和V0'＝[1543.3 -367.0 7528.4]T。結合其他可獲取的參數，可計算得到過境預報區域的衛星的激光質心在地球固定坐標系的坐標值P0＝[-1855244.6 4669501.6 4693461.4]T和位移速度值V0＝[-287.4 5397.1 -5468.8]T。
步驟3.2，根據下述公式計算軌道坐標系相對于地球固定坐標系的旋轉矩陣MT_o；
其中，P0＝[X0 Y0 Z0]T為衛星激光質心在地球固定坐標系的坐標值，V0＝[VX0 VY0VZ0]T為衛星激光質心在地球固定坐標系的位移速度值，A1,A2,A3,B1,B2,B3,C1,C2,C3為旋轉矩陣MT_o的各個元素。
將步驟3.1得到的過境預報區域的衛星的激光質心在地球固定坐標系的坐標值P和位移速度值V代入到步驟3.2的公式，可以求得軌道坐標系相對于地球固定坐標系的旋轉矩陣
步驟3.3，根據下述公式將激光測高儀的指向由軌道坐標系轉換到地球固定坐標系；
MT_l＝MT_o*Mo_l
式中，Mo_l表示激光測高儀在軌道坐標系下的指向，MT_o表示軌道坐標系相對于地球固定坐標系的旋轉矩陣，MT_l表示激光測高儀在地球固定坐標系下的指向。
將步驟2.3和步驟3.2計算得到的激光測高儀的指向由本體坐標系轉換到軌道坐標系的旋轉矩陣與軌道坐標系相對于地球固定坐標系的旋轉矩陣相乘，可以得到激光測高儀的指向由本體坐標系轉換到地球固定坐標系的向量
步驟4，通過計算地球固定坐標系下的激光測高儀的指向與地球橢球的交點，確定星載激光測高儀的足印位置。
根據光束的物理傳播路徑，在激光到達地面后，需要求地球固定坐標系下的激光測高儀的指向與地球橢球的交點，所述交點即本發明要確定的星載激光測高儀的足印位置，即激光探測器的布設位置的坐標。步驟4具體包括以下子步驟：
步驟4.1，根據下述公式計算激光指向與地球橢球相交的地球固定坐標系下的地面點坐標(X1,Y1,Z1)：
其中，MT＝[X1 Y1 Z1]T，X1、Y1和Z1為待求的激光指向與地球橢球相交的地球固定坐標系下的地面點坐標，μ為待求的比例系數，MT_l為激光測高儀在地球固定坐標系下的指向，P0為衛星激光質心在地球固定坐標系的坐標值，a和b分別為地球橢球長半軸和短半軸長度(a＝6378137.0米，b＝6356752.3米)，h為地面大地高程，可以通過其他輔助高程數據獲取，初始值給0。
其中，MT＝P0+μ*MT_l可以表示為：或表示為將X1、Y1、Z1代入到公式中，可以得到下述方程
式中，X0、Y0和Z0為步驟3.1計算得到衛星激光質心在地球固定坐標系的坐標值，MT_l為步驟3.3計算得到的值。
解上述一元二次方程，得到比例系數μ，因為方程會有兩個解，最終選擇絕對值小的u＝506437.3為最終解，根據MT＝P0+μ*MT_l從而可以計算出MT＝[-1718742.3 4325848.3 4347414.8]T。
步驟4.2，根據下述公式將地球固定坐標系下的地面點坐標(X1,Y1,Z1)轉換成大地經緯度坐標系下的地面點坐標(L,B,H)；
(X1,Y1,Z1)對于我們來說，沒法直接使用，需要轉換為我們經常用的經緯度格式(L,B,H)，這個也是從天上到地面的必要轉換。
其中，e1和e2分別為第一偏心率和第二偏心率，a和b分別為地球橢球長半軸和短半軸長度，X1、Y1、Z1為通過步驟4.1計算出的激光指向與地球橢球相交的地球固定坐標系下的地面點坐標。
根據上述公式，可以求得地面坐標(L,B,H)＝(111.66887,43.23643,1079.99)。
根據本發明的優選實施例，為了得到更精確的星載激光測高儀的足印位置，步驟4還可以包括：
步驟4.3，通過迭代計算，獲得精確的激光探測器的布設位置；步驟4.3具體包括：將步驟4.2計算得到的地面點坐標(L,B)代入到高程數據中，重新內插得到更精確的高程值H'，將計算得到的H'代入到步驟4.1和步驟4.2進行迭代運算，以得到精確的星載激光測高儀的足印位置。
具體來說，將步驟4.2得到的地面點坐標(L,B)代入到高程數據中，重新內插得到更精確的高程值H'，內插公式采用雙線性內插：
H'＝(1-Δx)(1-Δy)H1+(1-Δx)ΔyH2+Δx(1-Δy)H3+ΔxΔyH4
式中，H1、H2、H3和H4為H'附近x、y方向上相鄰點的高程值，Δx為H'距H3和H4的x方向距離差，Δy為H'距H1和H2的y方向距離差。
將計算得到的H'代入到步驟4.1和步驟4.2，重復計算，此過程一般迭代三次，最終得到精確的星載激光測高儀的足印位置(L,B,H)，即，足印中心點的坐標。
將步驟4.2計算得到的地面點坐標(111.66887,43.23643)代入到我們使用的AW3D30數據(公開數據)中，采用上述內插公式，可以內插得到新的高程值H＝951.21；將H＝951.21代入到步驟4.1和步驟4.2得到新的地面點坐標(111.66879,43.23638，951.21),將(111.66879,43.23638)代入到AW3D30數據中，可以內插得到新的高程值H＝950.14；將H＝950.14代入到步驟4.1和步驟4.2得到新的地面點坐標(111.66879,43.23638，950.08)，該位置即為待鋪設探測器陣的中心點。
本發明結合了傳統的光學相機成像的原理，并充分顧及了激光測高儀自身的測高機理，預報結果能夠較好地符合衛星真實的飛行過境情況，具有較高的預報精度。
以上內容僅為本發明的較佳實施例，對于本領域的普通技術人員，依據本發明的思想，在具體實施方式及應用范圍上均會有改變之處，本說明書內容不應理解為對本發明的限制。
- 還沒有人留言評論。精彩留言會獲得點贊！