本發(fā)明涉及雷達(dá)干涉處理技術(shù)領(lǐng)域，特別涉及一種基于SAR嚴(yán)密成像模型的子像素級角反射器自動(dòng)定位方法。

背景技術(shù)：

合成孔徑雷達(dá)(SAR)技術(shù)誕生于1951年，它通過收集和記錄目標(biāo)、背景在微波波段內(nèi)的電磁波輻射、散射能量，利用信號處理技術(shù)加工整理，提取識別目標(biāo)物體或現(xiàn)象的有用信息，廣泛應(yīng)用于測繪地形地物、檢測海洋、探測地球資源和偵察監(jiān)視軍事目標(biāo)等。角反射器(CR)在合成孔徑雷達(dá)影像中具有高反射回波、強(qiáng)穩(wěn)定性等特點(diǎn)，已應(yīng)用于SAR影像的幾何糾正和合成孔徑雷達(dá)干涉測量(InSAR)檢校中。

目前對角反射器的設(shè)計(jì)工作主要集中在如角反射器面的選擇、設(shè)計(jì)、加工公差、角度誤差、入射角對準(zhǔn)誤差等，但是涉及對角反射器的自動(dòng)定位的較少，尤其現(xiàn)有技術(shù)中并沒有子像素級的反射器的自動(dòng)定位?，F(xiàn)有的反射器定位技術(shù)多是利用角反射器點(diǎn)具有較強(qiáng)的回波信號與周圍地物具有相對較弱的回波信號特性，利用一定大小的窗口模板去檢測、匹配圖像，計(jì)算相應(yīng)的相關(guān)系數(shù)，相關(guān)系數(shù)高的點(diǎn)即為角反射器所在的位置，這種方法計(jì)算速度較快，然而其探測精度只能到像素級。此外，現(xiàn)有技術(shù)中還有基于提取PS點(diǎn)的方法識別角反射器點(diǎn)，計(jì)算角反射器點(diǎn)周圍一定窗口內(nèi)的相干系數(shù)、振幅離差值，具有高的信噪比點(diǎn)即為角反射器點(diǎn)位置，這種方法優(yōu)于相關(guān)系數(shù)法，精度可控制在1-2像素以內(nèi)，然而影像數(shù)目動(dòng)輒數(shù)十景，成本過高；另外，現(xiàn)有技術(shù)中還有基于信噪比(SCR)的角反射器提取方法，該方法首先估算初始定位的最大偏移，確定角反射器目標(biāo)的大致范圍，然后繪制確定的角反射器目標(biāo)范圍的SCR值圖，SCR值最大的點(diǎn)即為角反射器目標(biāo)。然而，現(xiàn)有技術(shù)的這三種方法由于缺少質(zhì)量評價(jià)與檢核方式，因此無法評估定位結(jié)果的好壞，一旦存在粗差，那么在高精度地形測繪中，會(huì)引入較大的解算誤差。

SAR圖像易受斑點(diǎn)噪聲的影響，如果所采用的窗口模板與實(shí)際角反射器影像特征不一致，則識別、檢測的精度不高；此外，角反射器擺放的位置、入射角、加工大小等都會(huì)對SAR圖像特征產(chǎn)生明顯的影響，比如雷達(dá)傳感器的入射角偏離理想位置，則角反射器的三面角效應(yīng)不明顯，實(shí)際獲取的角反射器特征無法預(yù)測，而現(xiàn)有的檢測角反射器算法則無法識別，無論是使用模板匹配法還是使用PS點(diǎn)的方法去提取所在位置，其檢測精度也不高，因此不具有通用性，且定位精度無法達(dá)到子像素級的精度。此外，使用窗口模板提取角反射器點(diǎn)，并不能有效提供角反射器的提取精度。角反射器點(diǎn)識別的對不對、精度高不高，無法評價(jià)。某些角反射器受到周圍強(qiáng)后向散射回波信號的干擾，會(huì)很容易導(dǎo)致定位錯(cuò)誤，從而影響角反射器在高精度測量檢校過程中的使用效果。

因此，需要提出一種能夠克服上述缺點(diǎn)的新型的反射器自動(dòng)定位方法。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明提出一種基于SAR嚴(yán)密成像模型的子像素級角反射器自動(dòng)定位方法，本發(fā)明采用了嚴(yán)密的SAR成像模型進(jìn)行角反射器自動(dòng)提取，以避免人工參與過程中的低效率和低精度問題。同時(shí)，針對角反射器自動(dòng)提取算法無法提供直觀的質(zhì)量評價(jià)方法的缺陷，本發(fā)明提供逐點(diǎn)誤差，并可剔除粗差點(diǎn)，提高點(diǎn)位提取結(jié)果的穩(wěn)健性，可充分地消除由于現(xiàn)有技術(shù)的限制和缺陷導(dǎo)致的一個(gè)或多個(gè)問題。

本發(fā)明另外的優(yōu)點(diǎn)、目的和特性，一部分將在下面的說明書中得到闡明，而另一部分對于本領(lǐng)域的普通技術(shù)人員通過對下面的說明的考察將是明顯的或從本發(fā)明的實(shí)施中學(xué)到。通過在文字的說明書和權(quán)利要求書及附圖中特別地指出的結(jié)構(gòu)可實(shí)現(xiàn)和獲得本發(fā)明目的和優(yōu)點(diǎn)。

本發(fā)明提供了一種基于SAR嚴(yán)密成像模型的子像素級角反射器自動(dòng)定位方法，其特征在于，所述方法具體包括以下步驟：

S1，采集原始數(shù)據(jù)，所述原始數(shù)據(jù)包括角反射器點(diǎn)的坐標(biāo)、主影像以及所述主影像的成像參數(shù)；

S2，利用采集的原始數(shù)據(jù)構(gòu)建SAR嚴(yán)密成像模型，并基于構(gòu)建的SAR嚴(yán)密成像模型計(jì)算所有角反射器點(diǎn)的初始像素坐標(biāo)(x₀,y₀)，得到所有角反射器點(diǎn)的初始像素坐標(biāo)的集合{P_init}；

S3，基于步驟S2獲取的角反射器點(diǎn)的初始像素坐標(biāo)集合{P_init}，使用最大振幅算法解算所有角反射器點(diǎn)在主影像中的像素級坐標(biāo)，并得到所有角反射器點(diǎn)在主影像中的像素級坐標(biāo)集合{(P_pix)}；

S4，對角反射器點(diǎn)的像素級坐標(biāo)集合{(P_pix)}中所有角反射器點(diǎn)的像素級坐標(biāo)進(jìn)行干涉參數(shù)測量校驗(yàn)，獲得角反射器點(diǎn)的子像素級坐標(biāo)。

本發(fā)明還提供了一種計(jì)算機(jī)可讀介質(zhì)，介質(zhì)中存有計(jì)算機(jī)程序，所述計(jì)算機(jī)程序用于執(zhí)行本發(fā)明所述的方法。

本發(fā)明兼顧了角反射器的影像特征以及幾何特征，可實(shí)現(xiàn)自動(dòng)化子像素級角反射器提取，該方法提高了點(diǎn)位提取效率，同時(shí)保證了點(diǎn)位坐標(biāo)的精度，實(shí)現(xiàn)了子像素級自動(dòng)化的提取角反射器點(diǎn)坐標(biāo)。

附圖說明

圖1為根據(jù)本發(fā)明實(shí)施例的基于SAR嚴(yán)密成像模型的子像素級角反射器自動(dòng)定位方法的流程圖。

圖2為根據(jù)本發(fā)明實(shí)施例的地心空間直角坐標(biāo)系下地形三維重建幾何關(guān)系。

具體實(shí)施方式

下面參照附圖對本發(fā)明進(jìn)行更全面的描述，其中說明本發(fā)明的示例性實(shí)施例。

SAR嚴(yán)密成像模型利用距離方程、多普勒方程和干涉相位方程求解目標(biāo)點(diǎn)的三維坐標(biāo)，描述了干涉SAR基本測量值和地面目標(biāo)點(diǎn)位置之間的關(guān)系，融合了有關(guān)地面目標(biāo)位置的所有信息，是公認(rèn)的嚴(yán)密成像模型。本發(fā)明提出的基于SAR嚴(yán)密成像模型的子像素級角反射器自動(dòng)定位方法，可綜合考慮影像的后向散射特征以及干涉幾何特征，獲取角反射器的精確定位結(jié)果，并給出每個(gè)角反射器的定位誤差，輔助進(jìn)行角反射器的篩選，進(jìn)一步保證干涉測量檢校方案的有效性。

如圖1所示，基于SAR嚴(yán)密成像模型的子像素級角反射器自動(dòng)定位方法，其可實(shí)現(xiàn)角反射器的高精度自動(dòng)化提取。該方法具體包括以下步驟：

步驟S1，采集原始數(shù)據(jù)，所述原始數(shù)據(jù)包括角反射器點(diǎn)的坐標(biāo)、主影像、從影像以及所述主影像和從影像的成像參數(shù)。

根據(jù)本發(fā)明的實(shí)施例，可以適用于任何星載/機(jī)載影像；優(yōu)選的，本發(fā)明的實(shí)施案例以機(jī)載Ka波段干涉SAR影像為例，但并不限于此。數(shù)據(jù)來源于航天科工二院，飛行時(shí)間為2015年3月22日，飛行區(qū)域?yàn)殛兾魇￠惲际?。?yōu)選的，本發(fā)明以航帶中的兩景影像為例。其主影像、從影像成像參數(shù)如表1所示。眾所周知，在星載/機(jī)載SAR影像數(shù)據(jù)處理中，為了獲取高精度的數(shù)據(jù)處理結(jié)果，需要利用控制點(diǎn)數(shù)據(jù)精化SAR處理過程，因此往往通過布設(shè)三面角角反射器點(diǎn)來獲得穩(wěn)定、可靠、易于量測的控制點(diǎn)坐標(biāo)，在后續(xù)表述中均使用“角反射器點(diǎn)”表達(dá)“控制點(diǎn)”，提供的角反射器點(diǎn)坐標(biāo)格式以最為常用的經(jīng)度、緯度和高程形式表示，是測繪遙感領(lǐng)域最為常用的坐標(biāo)表達(dá)形式。在飛機(jī)飛行過程中，同時(shí)布設(shè)了角反射器并采用GPS進(jìn)行了坐標(biāo)測量，GPS測量坐標(biāo)如表2所示，角反射器點(diǎn)坐標(biāo)包括：第一列表示角反射器點(diǎn)的名稱，第二列表示經(jīng)度，第三列表示緯度和第四列表示高程信息。

表1主、從影像成像參數(shù)

表2角反射器點(diǎn)的大地坐標(biāo)

步驟S2，利用采集的原始數(shù)據(jù)構(gòu)建SAR嚴(yán)密成像模型，并基于構(gòu)建的SAR嚴(yán)密成像模型計(jì)算所有角反射器點(diǎn)的初始像素坐標(biāo)(x₀,y₀)，得到所有角反射器點(diǎn)的初始像素坐標(biāo)的集合{P_init}；下面以針對采集的主影像數(shù)據(jù)為例。

步驟S2主要涉及兩部分，第一部分是將角反射器點(diǎn)的地理坐標(biāo)轉(zhuǎn)換到地心空間直角坐標(biāo)，第二部分是獲取角反射器點(diǎn)的初始像素坐標(biāo)，這個(gè)坐標(biāo)表達(dá)是在雷達(dá)坐標(biāo)系下描述的，描述了影像空間坐標(biāo)，角反射器點(diǎn)的初始像素坐標(biāo)是進(jìn)行后續(xù)坐標(biāo)精化的基礎(chǔ)。初始像素坐標(biāo)的解算精度要求并不是很高，一般來說，角反射器定位精度可以達(dá)到數(shù)十個(gè)像素就可以。解算過程中一般采用斜距-多普勒方程。它描述了成像過程中的斜距、多普勒、地球橢球的幾何關(guān)系，是SAR嚴(yán)密成像模型。

具體來說，步驟S2具體包括以下子步驟：

步驟S2.1，將所述角反射器點(diǎn)的大地坐標(biāo)轉(zhuǎn)換到地心空間直角坐標(biāo)；

角反射器點(diǎn)的大地坐標(biāo)是野外測量人員實(shí)測的地面點(diǎn)地理坐標(biāo)，將它轉(zhuǎn)換到地心空間直角坐標(biāo)是因?yàn)樵诤罄m(xù)的處理中，很多處理步驟是在地心空間直角坐標(biāo)系下完成，因此需要進(jìn)行坐標(biāo)轉(zhuǎn)換。

大地坐標(biāo)系是采用經(jīng)緯高描述一個(gè)點(diǎn)的坐標(biāo)信息，其中大地經(jīng)度L是過地面點(diǎn)的橢球子午面與格林尼治天文臺(tái)子午面的夾角，大地緯度B是過地面點(diǎn)的橢球法線與橢球赤道面的夾角，大地高H是地面點(diǎn)沿橢球法線到地面橢球面的距離。而在SAR嚴(yán)密成像模型中，一般采用的是地心空間直角坐標(biāo)系，坐標(biāo)系原點(diǎn)與地球質(zhì)心重合，Z軸指向地球北極，X軸指向格林尼治平均自屋面與地球赤道的角點(diǎn)，Y軸垂直于XOZ平面構(gòu)成右手坐標(biāo)系?？梢妰勺鴺?biāo)系的定義完全不同，一般需要將大地坐標(biāo)轉(zhuǎn)換為地心空間直角坐標(biāo)，以便在同一坐標(biāo)框架下進(jìn)行計(jì)算。

具體的，根據(jù)下述公式(1)將角反射器點(diǎn)的大地坐標(biāo)轉(zhuǎn)換為地心空間直角坐標(biāo)：

其中，(X_P,Y_P,Z_P)為角反射器點(diǎn)的地心空間直角坐標(biāo)，N為卯酉圈曲率半徑，(B,L,H)分別為角反射器點(diǎn)的大地坐標(biāo)的緯度、經(jīng)度和大地高，e為參考橢球的第一偏心率。

以角反射器點(diǎn)D002為例來說明將角反射器點(diǎn)的大地坐標(biāo)轉(zhuǎn)換為地心空間直角坐標(biāo)的過程。通過步驟S1獲得的D002點(diǎn)的大地坐標(biāo)(B,L,H)為(110.0559,35.10238,667.631)，對應(yīng)的橢球?yàn)閃GS-84橢球，其橢球長半軸a為6378137m，第一偏心率e為6.6944×10-3，對應(yīng)的卯酉圈曲率半徑N可表達(dá)為N＝a(1-e²sin²B)^-1/2＝6385208.2m。獲取了D002點(diǎn)的大地坐標(biāo)(B,L,H)以及各項(xiàng)橢球參數(shù)之后，代入上述坐標(biāo)轉(zhuǎn)換公式可得：

步驟S2.2，計(jì)算主影像對應(yīng)的衛(wèi)星成像時(shí)刻t；具體的，根據(jù)下述公式(2)計(jì)算衛(wèi)星成像時(shí)刻t：

t＝t₀+Δt·y (2)

其中，t₀為主影像的起始成像時(shí)間(即，第一行影像的成像時(shí)間)，Δt為成像時(shí)間間隔，y為角反射器點(diǎn)在雷達(dá)坐標(biāo)系中的方位向坐標(biāo)。

需要說明的是，雷達(dá)坐標(biāo)系是以雷達(dá)飛行軌跡方向?yàn)榉轿幌蜃鴺?biāo)y，以雷達(dá)掃描方向?yàn)榫嚯x向坐標(biāo)x進(jìn)行二維成像的。

同樣的，以D002為例，說明衛(wèi)星成像時(shí)刻的計(jì)算過程。將通過步驟S1獲得的起始成像時(shí)間(12079.8304)以及成像時(shí)間間隔(0.00779883)代入上述公式(2)，可得：

t＝12079.8304+0.00779883y。

步驟S2.3，擬合衛(wèi)星位置矢量與所述衛(wèi)星成像時(shí)刻t之間的二階多項(xiàng)式，所述二階多項(xiàng)式如公式(3)所示：

其中，(X_S,Y_S,Z_S)為待擬合的衛(wèi)星位置矢量，為位置擬合二階多項(xiàng)式系數(shù)矩陣。由于衛(wèi)星一般沿規(guī)則軌道飛行，在短時(shí)間內(nèi)用二階多項(xiàng)式能夠完美擬合衛(wèi)星位置。

具體而言，為了得到多項(xiàng)式系數(shù)矩陣，需要獲取至少三組觀測量。以主影像為例，對通過步驟S1獲取的影像中共包含13組觀測量。每一組的觀測時(shí)間可以通過狀態(tài)矢量起始時(shí)間(12079.8304s)以及對應(yīng)的時(shí)間間隔(2.193425s)，計(jì)算得到：t_i＝12079.8304+2.193425*i(i＝0,1,2...,12)

以X方向坐標(biāo)為例，使用步驟S1獲取的13個(gè)狀態(tài)位置矢量可構(gòu)建13個(gè)方程，即

X＝(a₀ a₁ a₂)T

其中，X＝(-1790516.1 -1790657.7 -1790799.2 -1790940.7 -1791082.3 -1791223.8 -1791365.3 -1791506.9 -1791648.4 -1791789.9 -1791931.4 -1792073 -1792214.5)'，T＝(t₀ t₁ … t₁₂)′。解算可獲得多項(xiàng)式系數(shù)

(a_x0 a_x1 a_x2)＝(-992845.37 -67.5383 0.00012459)。

針對Y和Z分別重復(fù)上述步驟，即可構(gòu)建衛(wèi)星位置和成像時(shí)刻之間的關(guān)系矩陣：

步驟S2.4，擬合衛(wèi)星速度矢量與所述衛(wèi)星成像時(shí)刻t之間的二階多項(xiàng)式，所述二階多項(xiàng)式如公式(4)所示：

其中，為待擬合的衛(wèi)星速度矢量，t為衛(wèi)星成像時(shí)刻，為速度擬合二階多項(xiàng)式系數(shù)矩陣。

具體而言，為了得到多項(xiàng)式系數(shù)矩陣，需要獲取至少三組觀測量。以主影像為例，對通過步驟S1獲取的影像中共包含13組觀測量。每一組的觀測時(shí)間可以通過狀態(tài)矢量起始時(shí)間(12079.8304s)以及對應(yīng)的時(shí)間間隔(2.193425s)，計(jì)算得到：t_i＝12079.8304+2.193425*i(i＝0,1,2...,12)。

以V_X方向速度為例，使用步驟S1獲取的13個(gè)狀態(tài)位置矢量可構(gòu)建13個(gè)方程，即

V_X＝(b₀ b₁ b₂)T

其中，V_X＝(-64.552448 -64.551767 -64.579169 -64.578488 -64.577807 -64.549033 -64.548349 -64.547665 -64.546981 -64.684219 -64.655433 -64.654748 -64.544246)'，T＝(t₀ t₁ … t₁₂)′。解算可獲得多項(xiàng)式系數(shù)：

(b_x0 b_x1 b_x2)＝(-5693.1802 0.93345 -0.0000387)。

針對V_Y和V_Z分別重復(fù)上述步驟，即可構(gòu)建衛(wèi)星速度和成像時(shí)刻之間的關(guān)系矩陣：

步驟S2.5，構(gòu)建SAR嚴(yán)密成像模型，所述SAR嚴(yán)密成像模型如公式(5)所示：

其中，Vx_S·(X_S-X_P)+Vy_S·(Y_S-Y_P)+Vz_S·(Z_S-Z_P)＝-f_d·λ·(R₀+ΔR·x)/2為多普勒方程，為距離方程，(X_S,Y_S,Z_S)為衛(wèi)星位置矢量，(X_P,Y_P,Z_P)為角反射器點(diǎn)的地心空間直角坐標(biāo)，為衛(wèi)星速度矢量，x為角反射器點(diǎn)的距離向坐標(biāo)，R₀為近地點(diǎn)斜距，ΔR為距離向分辨率，f_d為多普勒中心頻率，λ為雷達(dá)波長。

對通過步驟S1獲取的表1所示的主影像為例，其近地點(diǎn)斜距R₀(4885.9532m)、距離向分辨率ΔR(0.5353m)、多普勒中心頻率f_d(560.9245)以及雷達(dá)波長(λ＝c/f＝3×10⁸/(3×10⁵)＝0.008565m)均為已知量。

步驟S2.6，使用迭代法解算步驟S2.5構(gòu)建的SAR嚴(yán)密成像模型，求解角反射器點(diǎn)的初始像素坐標(biāo)。

步驟S2.6具體包括以下子步驟：

步驟S2.6.1，基于公式(5)的多普勒方程迭代計(jì)算角反射器點(diǎn)的初始方位向坐標(biāo)y；所述步驟S2.6.1包括以下子步驟：

S2.6.1.1，設(shè)定角反射器點(diǎn)的像素坐標(biāo)初值，優(yōu)選的，設(shè)定的像素坐標(biāo)初值為SAR主影像的中心點(diǎn)(x_center,y_center)；

根據(jù)本發(fā)明的實(shí)施例，在表1提供的主影像參數(shù)列表里，行數(shù)為3376，列數(shù)為5400，因此主影像的中心點(diǎn)(x_center,y_center)＝(2700,1688)。

S2.6.1.2，根據(jù)角反射器點(diǎn)的方位向y_center計(jì)算衛(wèi)星成像時(shí)間t_y，并由擬合多項(xiàng)式計(jì)算相應(yīng)的衛(wèi)星傳感器位置矢量(X_S,Y_S,Z_S)和速度矢量

其中，根據(jù)上述步驟S2.2計(jì)算衛(wèi)星成像時(shí)間t_y，根據(jù)步驟S2.3和S2.4獲取的衛(wèi)星位置和速度矢量多項(xiàng)式，在此不再贅述。

根據(jù)本發(fā)明的實(shí)施例，方位向坐標(biāo)y_center＝1688，計(jì)算得到衛(wèi)星成像時(shí)間為t_y＝12092.9948s，衛(wèi)星傳感器的位置矢量為(X_S,Y_S,Z_S)＝(-1791365.573,4907370.081,3652941.435)，衛(wèi)星傳感器的速度矢量為

S2.6.1.3，利用衛(wèi)星傳感器位置矢量(X_S,Y_S,Z_S)、速度矢量和角反射器點(diǎn)位置矢量(X_P,Y_P,Z_P)，根據(jù)多普勒方程計(jì)算得到多普勒值f_d，同時(shí)計(jì)算得到多普勒變化率f_d'，結(jié)合多普勒測量值f_dop，根據(jù)下式計(jì)算出需要調(diào)整的時(shí)間變量Δt：

Δt＝-(f_d-f_dop)/f_d'

根據(jù)本發(fā)明的實(shí)施例，計(jì)算得到需要調(diào)整的時(shí)間變量Δt＝-0.1403s。

其中，多普勒變化率f_d'是通過數(shù)值微分的方法近似求出，假設(shè)要計(jì)算當(dāng)前方位向坐標(biāo)y對應(yīng)的時(shí)間t_y的f_d'，設(shè)dt＝0.01s，利用多普勒方程分別計(jì)算t_y和t_y+dt處的多普勒值f_d(t_y)和f_d(t_y+dt)，則

f_d'＝(f_d(t_i+dt)-f_d(t_i))/dt

S2.6.1.4，利用時(shí)間變量Δt對衛(wèi)星成像時(shí)間t_y進(jìn)行修改，得到新的成像時(shí)刻t_y+1＝t_y+Δt；

根據(jù)本發(fā)明的實(shí)施例，得到新的成像時(shí)刻

t_y+1＝t_y+Δt＝12092.9948-0.1403＝12092.8544s。

S2.6.1.5，利用t_y+1進(jìn)行步驟S2.6.1.2-S2.6.1.4的迭代計(jì)算，直到Δt絕對值都小于閾值ε，迭代結(jié)束，輸出主影像成像時(shí)刻t_y；根據(jù)本發(fā)明的優(yōu)選實(shí)施例，閾值ε為10^-7。

根據(jù)本發(fā)明的實(shí)施例，計(jì)算得到主影像成像時(shí)刻t_y＝12092.8535s。

S2.6.1.6，由解算的主影像成像時(shí)刻t_y，根據(jù)下式計(jì)算角反射器點(diǎn)的方位向坐標(biāo)y：y＝(t_y-t₀)/Δt_y，其中，t₀為起始時(shí)間，Δt_y為方位向時(shí)間間隔(如表1所示)。

根據(jù)本發(fā)明的實(shí)施例，計(jì)算角反射器點(diǎn)的方位向坐標(biāo)y＝1669.8915。

步驟S2.6.2，根據(jù)公式(5)的距離方程計(jì)算角反射器點(diǎn)的初始距離向坐標(biāo)x。

由角反射器點(diǎn)的方位向坐標(biāo)y，計(jì)算衛(wèi)星傳感器位置，結(jié)合已知角反射器點(diǎn)地理坐標(biāo)計(jì)算衛(wèi)星傳感器到角反射器點(diǎn)的斜距值，進(jìn)而計(jì)算出地面角反射器點(diǎn)的距離向坐標(biāo)x。

所述步驟S2.6.2包括以下子步驟：

S2.6.2.1，由步驟S2.6.1獲得的角反射器點(diǎn)的初始方位向坐標(biāo)y，根據(jù)下述公式計(jì)算衛(wèi)星傳感器位置矢量(X_S,Y_S,Z_S)：

其中，t為當(dāng)前計(jì)算主影像對應(yīng)的衛(wèi)星成像時(shí)刻，t₀為主影像的起始成像時(shí)間(即，第一行影像的成像時(shí)間)，Δt為成像時(shí)間間隔，y為角反射器點(diǎn)在雷達(dá)坐標(biāo)系中的方位向坐標(biāo)，(X_S,Y_S,Z_S)為待擬合的衛(wèi)星位置矢量，為位置擬合二階多項(xiàng)式系數(shù)矩陣。

根據(jù)本發(fā)明的實(shí)施例，已知角反射器點(diǎn)的方位向坐標(biāo)y＝1669.89，計(jì)算得到衛(wèi)星傳感器位置矢量(X_S,Y_S,Z_S)＝(-1791356.456,4907373.417,3652941.425)。

S2.6.2.2，由距離方程結(jié)合角反射器點(diǎn)的位置矢量(X_P,Y_P,Z_P)、近地點(diǎn)斜距R₀、距離向分辨率ΔR，計(jì)算角反射點(diǎn)的初始距離向坐標(biāo)x。

優(yōu)選的，由角反射器點(diǎn)的位置矢量(X_P,Y_P,Z_P)＝(-1791657.3,4907631.1,3647548.8)、近地點(diǎn)斜距R₀＝4885.9532m、距離向分辨率ΔR＝0.5353m，計(jì)算得到角反射點(diǎn)的初始距離向坐標(biāo)x＝977.72736。

如上所述，根據(jù)計(jì)算獲取的角反射器點(diǎn)的方位向坐標(biāo)y和距離向坐標(biāo)x，最終求得角反射器點(diǎn)的初始像素坐標(biāo)(x,y)。以角反射器點(diǎn)D002為例，解算出來的初始像素坐標(biāo)結(jié)果為(977.72736,1669.8915)。

由步驟2.2–步驟2.6可知，在構(gòu)建的SAR嚴(yán)密成像模型中，衛(wèi)星位置矢量(X_S,Y_S,Z_S)(13個(gè)離散點(diǎn)衛(wèi)星位置信息)、衛(wèi)星速度矢量(13個(gè)離散點(diǎn)衛(wèi)星速度信息)、衛(wèi)星成像時(shí)間段(起始時(shí)間t₀，中間時(shí)間t₁，終止時(shí)間t₂)t₀-t₂、角反射器點(diǎn)的大地坐標(biāo)(B,L,H)、橢球參數(shù)(如第一偏心率e，卯酉圈曲率半徑N)、雷達(dá)波長λ、近地點(diǎn)斜距R₀、距離向分辨率ΔR、多普勒中心頻率f_d等參數(shù)或中間輔助變量均通過步驟S1獲得，具體參數(shù)如表1、2所示。在公式(5)中涉及到衛(wèi)星位置矢量和衛(wèi)星速度矢量，它們是可以根據(jù)公式(2)和(3)計(jì)算得到，在公式(2)中y為角反射器在雷達(dá)坐標(biāo)系中的方位向坐標(biāo)。所以在公式(5)中，待求解的未知參數(shù)只有兩個(gè)(x,y)，由于是第一次獲取的初始像點(diǎn)坐標(biāo)，所以計(jì)算得到的(x,y)是角反射器的初始坐標(biāo)，如表3所示，以(x₀,y₀)表示根據(jù)公式(5)計(jì)算得到的角反射器點(diǎn)的初始像素坐標(biāo)。

此外，步驟S1采集的角反射器點(diǎn)的大地坐標(biāo)(B,L,H)，為野外實(shí)測地理坐標(biāo)，并非角反射器點(diǎn)的初始像素坐標(biāo)，所以經(jīng)過坐標(biāo)轉(zhuǎn)換為空間直角坐標(biāo)(X_P,Y_P,Z_P)，參與步驟S2.5的解算過程。使用步驟S2.5的兩個(gè)方程，可解算這兩個(gè)未知參數(shù)(x,y)。

步驟S2.7，重復(fù)步驟S2.1-S2.6，計(jì)算所有角反射器點(diǎn)的初始像素坐標(biāo)，得到所有角反射器點(diǎn)的初始像素坐標(biāo)的集合{P_init}。計(jì)算出的各個(gè)角反射器點(diǎn)的初始像素坐標(biāo)見表3。

表3角反射器點(diǎn)的初始像素坐標(biāo)

步驟S3，基于步驟S2獲取的角反射器點(diǎn)的初始像素坐標(biāo)集合{P_init}，使用最大振幅算法解算所有角反射器點(diǎn)在主影像中的像素級坐標(biāo)，并得到所有角反射器點(diǎn)在主影像中的像素級坐標(biāo)集合{(P_pix)}；

步驟S2的精度一般在數(shù)十個(gè)像素，這對于精確的角反射器定位來說是遠(yuǎn)遠(yuǎn)不夠的。此時(shí)一般需要采用SAR影像本身的特性，將坐標(biāo)精度控制在數(shù)個(gè)像素以內(nèi)?？紤]到角反射器在影像上屬于高亮的點(diǎn)目標(biāo)，因此可采用影像振幅進(jìn)行像素級定位。定位過程中，需要采用快速傅立葉變換進(jìn)行插值，因此所有用到的窗口尺寸均為2的指數(shù)。解算過程中，依然以D002為例，其他角反射器像素級坐標(biāo)解算方法與D002一致。

步驟S3具體包括以下子步驟：

步驟S3.1，以通過步驟S2獲取的角反射器點(diǎn)的初始像素坐標(biāo)(x₀,y₀)為中心，獲取窗口大小為W₁×W₁的復(fù)數(shù)矩陣；

以D002為例，以通過步驟S2獲得的初始像素坐標(biāo)(977.72736,1669.8915)為中心，獲取窗口大小為W₁×W₁的復(fù)數(shù)矩陣。優(yōu)選的，W₁＝128。

步驟S3.2，計(jì)算所述W₁×W₁的復(fù)數(shù)矩陣的振幅值，得到第一振幅值矩陣。

具體的，根據(jù)下述公式計(jì)算所述W₁×W₁的復(fù)數(shù)矩陣的振幅值：

其中，Re(i,j)表示W(wǎng)₁×W₁的復(fù)數(shù)矩陣的像元值(i,j)的實(shí)部，Im(i,j)表示W(wǎng)₁×W₁的復(fù)數(shù)矩陣的像元值(i,j)的虛部。

優(yōu)選的，以W₁＝128為例，獲取128×128的復(fù)數(shù)影像，對復(fù)數(shù)矩陣中每一個(gè)像元值(i,j)進(jìn)行振幅值計(jì)算，計(jì)算振幅值公式為最后得到128×128的振幅值矩陣。

步驟S3.3，找出所述第一振幅值矩陣中的最大振幅值，并計(jì)算所述最大振幅值所在點(diǎn)在主影像中的像素坐標(biāo)

例如，在128×128的振幅值矩陣中，以第一個(gè)像素點(diǎn)g(1,1)為最大振幅值點(diǎn)，與剩余所有振幅值點(diǎn)進(jìn)行比較g(i,j)，其中i＝2…128,j＝2…128，如果g(i,j)大于g(1,1)，則將g(i,j)振幅值賦值給g(1,1)，繼續(xù)參與后續(xù)計(jì)算，直至所有點(diǎn)遍歷結(jié)束，輸出最大振幅值點(diǎn)。

對通過步驟S3.2獲得的第一振幅矩陣，計(jì)算最大值點(diǎn)在原始影像中的像素坐標(biāo)振幅最大值點(diǎn)提供的是角反射器的進(jìn)一步精化坐標(biāo)。如果不考慮粗差的影響，振幅最大值點(diǎn)提供的角反射器定位精度可達(dá)到一個(gè)像素。為了達(dá)到子像素級精度，一般需要對坐標(biāo)進(jìn)行進(jìn)一步精化，即采用插值的方式獲取子像素級坐標(biāo)。

步驟S3.4，以所述像素坐標(biāo)為中心，獲取窗口大小為W2×W₂的復(fù)數(shù)矩陣。

舉例來說，以步驟S3.3獲得的D002坐標(biāo)(1001,1638)為中心，獲取窗口大小為W₂×W₂的復(fù)數(shù)矩陣。因?yàn)檫^小的窗口尺寸將導(dǎo)致插值算法的結(jié)果出現(xiàn)偏差，所以，此時(shí)選取的窗口尺寸應(yīng)小于步驟S3.3選取的窗口尺寸，即，W₂<W₁。優(yōu)選的，W₂＝32。

步驟S3.5，計(jì)算所述W₂×W₂的復(fù)數(shù)矩陣的振幅值，得到第二振幅值矩陣；

具體的，根據(jù)下述公式計(jì)算所述W₂×W₂的復(fù)數(shù)矩陣的振幅值：

其中，Re(i,j)表示W(wǎng)₂×W₂的復(fù)數(shù)矩陣的像元值(i,j)的實(shí)部，Im(i,j)表示W(wǎng)₂×W₂的復(fù)數(shù)矩陣的像元值(i,j)的虛部。

優(yōu)選的，以W₂＝32為例，獲取32×32的復(fù)數(shù)影像，對復(fù)數(shù)矩陣中每一個(gè)像元值(i,j)進(jìn)行振幅值計(jì)算，計(jì)算振幅值公式為最后得到32×32的振幅值矩陣。

步驟S3.6，對所述第二振幅矩陣中的振幅值進(jìn)行N₁倍過采樣，得到逆變換矩陣，并計(jì)算所述逆變換矩陣中每一像元的振幅值，得到第三振幅矩陣；

為了保證解算速度以及解算精度，一般采用快速傅立葉變換進(jìn)行過采樣。其基本原理在于，影像中的突變性信息(例如紋理、噪聲、斑點(diǎn)等)多為高頻信息，傅立葉變換之后，這些高頻信息在影像邊緣，低頻信息在影像中心。在傅氏空間無論如何變換，只要保證高頻以及低頻信息不損失，信號特征就可以保留下來。步驟S3.6具體包括以下子步驟：

步驟S3.6.1，對所述第二振幅矩陣中的振幅值進(jìn)行快速傅里葉變換，轉(zhuǎn)到復(fù)數(shù)空間；

步驟S3.6.2，對轉(zhuǎn)到復(fù)數(shù)空間的第二振幅矩陣中的振幅值進(jìn)行N₁倍過采樣；

將所述復(fù)數(shù)空間均分為四份，四份的左上角與右下角坐標(biāo)分別為：A[(1,1),(W₂/2,W₂/2)],B[(W₂/2+1,1),(W₂,W₂/2)],C[(1,W₂/2+1),(W₂/2,W₂)],D[(W₂/2+1,W₂/2+1),(W₂,W₂)]，同時(shí)創(chuàng)建一個(gè)空矩陣M，用于存儲(chǔ)經(jīng)快速傅里葉變換的頻域信號，最后對通過賦值之后的矩陣M進(jìn)行快速傅里葉逆變換，得到逆變換之后的復(fù)數(shù)矩陣，這個(gè)過程即為對第二振幅矩陣中的振幅值進(jìn)行N₁倍過采樣；

步驟S3.6.2具體包括以下子步驟：

S3.6.2.1，將所述復(fù)數(shù)空間均分為四份，四份的左上角與右下角坐標(biāo)分別為：A[(1,1),(W₂/2,W₂/2)],B[(W₂/2+1,1),(W₂,W₂/2)],C[(1,W₂/2+1),(W₂/2,W₂)],D[(W₂/2+1,W₂/2+1),(W₂,W₂)]；

優(yōu)選的，以D002角反射器點(diǎn)為例，其四組坐標(biāo)分別為A[(1,1),(16,16)],B[(17,1),(32,16)],C[(1,17),(16,32)],D[(17,17),(32,32)]。

S3.6.2.,2，創(chuàng)建空矩陣M，其用于存儲(chǔ)經(jīng)快速傅里葉變換的頻域信號，其中，所述空矩陣M的大小為(W₂×N₁)×(W₂×N₁)，其中N₁為過采樣倍數(shù)，將所述空矩陣的左上、右上、左下、右下四角處的(W₂/2)×(W₂/2)像素分別賦值為A、B、C、D，賦值之后的矩陣M參與下一步的傅里葉逆變換。

空矩陣是用來存放經(jīng)過快速傅里葉變換的頻域信號的。其根本原理是，頻率域與空間域相互轉(zhuǎn)換的過程中，矩陣的維數(shù)并不會(huì)發(fā)生變化。因此，為了保證過采樣之后的矩陣大小滿足要求，空矩陣M的大小必須與過采樣之后的矩陣大小一致。

優(yōu)選的，所述空矩陣M的左上、右上、左下、右下四角中的每個(gè)角占據(jù)的像素塊為16×16，此外，過采樣倍數(shù)N₁＝4，則(W₂×N₁)×(W₂×N₁)＝128×128。A、B、C、D塊的大小均為16×16，均含有256個(gè)數(shù)據(jù)值，以A塊為例，由于A矩陣塊是復(fù)數(shù)數(shù)據(jù)，所以有[(3.14,4.48)_ii＝1,(5.01,1.70)_ii＝2,(0.43,-1.94)_ii＝3,…,(-4.51,1.70)_ii＝255,(0.02,0.09)_ii＝256]，其中ii＝1到256。

S3.6.2.3，對通過賦值之后的空矩陣M進(jìn)行快速傅里葉逆變換，得到逆變換矩陣，其中，所述逆變換矩陣仍為復(fù)數(shù)矩陣形式；

上述三個(gè)子步驟的處理過程，即為對第二振幅矩陣中的振幅值進(jìn)行N₁倍過采樣。

步驟S3.6.3，計(jì)算所述逆變換矩陣中每一像元的振幅值，得到第三振幅矩陣。

具體的，根據(jù)下述公式計(jì)算所述逆變換矩陣的振幅值：

其中，Re(i,j)表示所述逆變換矩陣的像元值(i,j)的實(shí)部，Im(i,j)表示所述逆變換矩陣的像元值(i,j)的虛部。

步驟S3.7，找出所述第三振幅值矩陣中的最大振幅值，并計(jì)算所述最大振幅值所在點(diǎn)在主影像中的像素級坐標(biāo)(x₁,y₁)。

仍然以角反射器點(diǎn)D002為例，通過上述步驟計(jì)算得到的像素坐標(biāo)(x₁,y₁)＝(1000.5,1638.0625)。

步驟S3.8，對所有角反射器點(diǎn)重復(fù)步驟S3.1–S3.7，獲得所有角反射器點(diǎn)在主影像中的像素級坐標(biāo)集合{(P_pix)}，具體如表4所示。

表4角反射器點(diǎn)像素級影像坐標(biāo)

步驟S4，對角反射器點(diǎn)的像素級坐標(biāo)集合{(P_pix)}中所有角反射器點(diǎn)的像素級坐標(biāo)進(jìn)行干涉參數(shù)測量校驗(yàn)，獲得角反射器點(diǎn)的子像素級坐標(biāo)；

本步驟根據(jù)步驟S2構(gòu)建的SAR嚴(yán)密成像模型并結(jié)合干涉相位方程來進(jìn)行處理，三個(gè)方程組(SAR嚴(yán)密成像模型的兩個(gè)方程組和干涉相位方程)能夠在地心空間直角坐標(biāo)系下描述地面目標(biāo)點(diǎn)的三維坐標(biāo)，稱為地形三維重建模型。而在步驟S2.1中得到的角反射器點(diǎn)坐標(biāo)同樣在地心空間直角坐標(biāo)系下，據(jù)此可以用地形三維重建模型評價(jià)步驟S3獲取的角反射器點(diǎn)的像素級影像坐標(biāo)定位精度，獲取最終的角反射器點(diǎn)的子像素級坐標(biāo)(x₂,y₂)。

通常情況下，地形三維重建模型易受干涉系統(tǒng)參數(shù)誤差的影響，需要利用高精度的角反射器點(diǎn)(即控制點(diǎn))進(jìn)行干涉測量檢校，獲取精度較高的干涉系統(tǒng)參數(shù)。然而，角反射器點(diǎn)易受SAR斑點(diǎn)噪聲的影響，獲取的影像像素坐標(biāo)存在一定的誤差，可以利用地形三維重建模型進(jìn)行角反射器點(diǎn)像素精度的評價(jià)并剔除誤差較大的點(diǎn)。基于上述分析，步驟S4首先構(gòu)建地形三維重建模型的基本關(guān)系式，并基于地形三維重建模型推導(dǎo)三維重建模型對干涉系統(tǒng)參數(shù)的偏導(dǎo)數(shù)，據(jù)此構(gòu)建干涉測量檢校模型；隨后將通過步驟S3得到的像素級影像坐標(biāo)集合{P_pix}對干涉測量檢校模型進(jìn)行最小二乘迭代計(jì)算，優(yōu)化并輸出干涉系統(tǒng)參數(shù)，隨后統(tǒng)計(jì)每個(gè)角反射器點(diǎn)平差前后的精度，剔除誤差較大的點(diǎn)，再次進(jìn)行迭代計(jì)算，直至所有點(diǎn)滿足設(shè)定要求，輸出角反射器點(diǎn)的子像素級定位結(jié)果。

步驟S4具體包括以下子步驟：

步驟S4.1，基于SAR嚴(yán)密成像模式和干涉相位方程，確定基于視向量正交分解的三維重建模型；

在該步驟中，首先確定三維重建模型的基本關(guān)系式，所述三維重建模型的基本關(guān)系式為如公式(1)所示的包括距離方程、多普勒方程和干涉相位方程的方程組，然后，根據(jù)距離方程、多普勒方程和干涉相位方程，推導(dǎo)出基于視向量正交分解的三維重建模型。

其中，Vx_S·(X_S-X_P)+Vy_S·(Y_S-Y_P)+Vz_S·(Z_S-Z_P)＝-f_d·λ·(R₀+ΔR·x)/2為多普勒方程，為距離方程，為干涉相位方程，(X_S,Y_S,Z_S)為衛(wèi)星位置矢量，(X_P,Y_P,Z_P)為角反射器點(diǎn)的地心空間直角坐標(biāo)，為衛(wèi)星速度矢量，x為角反射器點(diǎn)的距離向坐標(biāo)，R₀為近地點(diǎn)斜距，ΔR為距離向分辨率，f_d為多普勒中心頻率，λ為雷達(dá)波長，為絕對干涉相位，φ為纏繞干涉相位，k為整周期數(shù)，Q為干涉模式，為衛(wèi)星主天線到目標(biāo)點(diǎn)的距離，為衛(wèi)星副天線到目標(biāo)點(diǎn)的距離。

距離方程、多普勒方程和干涉相位方程，它們構(gòu)成的方程組是高度非線性的，直接通過數(shù)值解法對它們求解是非常困難的，因此就需要通過分析InSAR幾何模型(InSAR幾何模型里包含三個(gè)方程：距離方程、多普勒方程和干涉相位方程)來將它們(三個(gè)方程)通過三維重建模型表達(dá)形式進(jìn)行求解。而視向量正交分解方法避免了它們之間的非線性數(shù)值求解關(guān)系，具有概念明確，求解簡便等特點(diǎn)?；诖?，本發(fā)明采用視向量正交分解法則，通過視向量的正交分解和坐標(biāo)轉(zhuǎn)換的形式獲取目標(biāo)點(diǎn)的三維坐標(biāo)解析解。

如圖2所示，本發(fā)明基于SAR嚴(yán)密成像模式和干涉相位方程，所確定的基于視向量正交分解的三維重建模型為：

其中，P為目標(biāo)點(diǎn)坐標(biāo)，S₁為衛(wèi)星位置矢量，r₁為視向量，表示主影像對應(yīng)的衛(wèi)星位置到地面目標(biāo)的距離(其中，即，為視向量r₁的模長)，為視向量r₁經(jīng)歸一化后得到的單位視向量。

其中，目標(biāo)點(diǎn)坐標(biāo)P由主影像對應(yīng)的衛(wèi)星位置矢量S₁和視向量r₁所構(gòu)成的三角形之和表示，需計(jì)算主影像對應(yīng)的衛(wèi)星位置矢量S₁，主影像對應(yīng)的衛(wèi)星到地面目標(biāo)的距離r_s1，以及單位視向量

眾所周知，三維重建模型能夠描述角反射器點(diǎn)影像像點(diǎn)坐標(biāo)與地面點(diǎn)地理坐標(biāo)之間的數(shù)學(xué)關(guān)系，然而由于三維重建模型所需要的求解參數(shù)存在誤差，需要通過三維重建模型確定誤差項(xiàng)，并通過三維重建模型構(gòu)建誤差模型(即，干涉測量檢校模型)，以此來消除誤差項(xiàng)，提供準(zhǔn)確的干涉系統(tǒng)參數(shù)信息參與后續(xù)計(jì)算。在下面的步驟中，首先對步驟S4.1確定的三維重建模型的各個(gè)變量進(jìn)行推導(dǎo)，得到其解析解表達(dá)形式(具體包括步驟S4.2-S4.4)；其次，構(gòu)成三維重建模型的解析解表達(dá)形式的參數(shù)必然存在誤差，就需要通過三維重建模型確定主要誤差源，在X、Y、Z三個(gè)方向上推導(dǎo)三維重建模型函數(shù)對主要誤差源的偏導(dǎo)數(shù)，構(gòu)建干涉測量檢校模型，優(yōu)化干涉系統(tǒng)參數(shù)(具體包括步驟S4.5-S4.8)；最后利用三維重建模型和干涉測量檢校模型，對干涉系統(tǒng)參數(shù)進(jìn)行優(yōu)化，剔除粗差，重復(fù)步驟S4.5-S4.8，直至所有角反射器點(diǎn)均滿足閾值要求，輸出角反射器點(diǎn)的子像素級坐標(biāo)(具體包括步驟S4.9-S4.12)。下面對各個(gè)步驟進(jìn)行詳細(xì)說明。

步驟S4.2，計(jì)算主影像對應(yīng)的衛(wèi)星位置矢量S₁，計(jì)算公式為其中，(X_S,Y_S,Z_S)為待擬合的衛(wèi)星位置矢量，為位置擬合二階多項(xiàng)式系數(shù)矩陣，t為主影像對應(yīng)的衛(wèi)星成像時(shí)刻。

優(yōu)選的，根據(jù)步驟S2.3計(jì)算得到擬合后的衛(wèi)星位置矢量

步驟S4.3，計(jì)算主影像對應(yīng)的衛(wèi)星位置到地面目標(biāo)的距離和單位視向量

按照視向量正交分解法，在VNW移動(dòng)坐標(biāo)系中，將單位視向量分解為在V軸、N軸和W軸上的分量，具體表示為：其中，r_v為單位視向量在V軸上的分量，r_n為單位視向量在N軸上的分量，r_w為單位視向量在W軸上的分量。

即，步驟S4.3計(jì)算的是和單位視向量

視向量正交分解法的基本思想是在以天線相位中心為原點(diǎn)的VNW移動(dòng)坐標(biāo)系下，對單位視向量進(jìn)行正交分解。

VNW移動(dòng)坐標(biāo)系的定義：以主天線所在相位中心為坐標(biāo)原點(diǎn)，V軸定義為載機(jī)速度矢量方向，W軸定義為載機(jī)速度矢量與干涉基線矢量的叉積，N軸由V軸和W軸的叉積共同確定，構(gòu)成右手坐標(biāo)系。在干涉SAR數(shù)據(jù)處理中，移動(dòng)坐標(biāo)系一般用于求解單位視向量。通常，將三個(gè)正交基命名為v、n和w，依次表示V軸、N軸和W軸，簡稱為VNW系。

VNW移動(dòng)坐標(biāo)系下的三個(gè)正交基其表達(dá)式為：

上式中，b為干涉基線矢量，v為載機(jī)飛行速度矢量。

步驟S4.3具體包括以下子步驟：

步驟S4.3.1，基于單位視向量計(jì)算單位視向量在V軸上的分量r_v，其表達(dá)式為：

其中，v為衛(wèi)星飛行速度矢量，|v|為速度模長，r為視向量，為視向量r的模長。

步驟S4.3.2，基于單位視向量計(jì)算單位視向量在N軸上的分量r_n，其表達(dá)式為：

其中，表示主天線到目標(biāo)點(diǎn)的距離，|b_pv|為垂直于衛(wèi)星速度方向的基線矢量模長，b_v為基線矢量在衛(wèi)星速度方向的投影值，r_v為單位視向量在V軸上的分量，B為基線模長，λ為雷達(dá)波長，為絕對干涉相位，Q為天線工作模式，其中，當(dāng)天線的工作模式為重復(fù)軌道模式時(shí)，Q＝2，當(dāng)天線的工作模式為一發(fā)雙收模式時(shí)，Q＝1。優(yōu)選的，在本實(shí)施例中，Q＝1。

步驟S4.3.3，根據(jù)角反射器點(diǎn)的像素坐標(biāo){P_pix}獲取角反射器點(diǎn)的絕對干涉相位其中φ₁為纏繞相位，φ₂為平地相位。纏繞相位和平地相位均由相應(yīng)的數(shù)據(jù)文件提供，根據(jù)步驟S3獲取的角反射器的像素坐標(biāo)信息讀取即可。獲取角反射器點(diǎn)的絕對干涉相位用于計(jì)算主、副天線相位中心到地面目標(biāo)點(diǎn)的距離差在步驟S4.3.2的公式(4)中參與計(jì)算。

步驟S4.3.4，在步驟S4.3中描述了單位視向量由于為單位視向量，其表達(dá)式符合右手螺旋法則，單位視向量在W軸上的分量r_w表示為：

以D002角反射器點(diǎn)為例，對通過步驟S1獲得的成像參數(shù)(包括主影像起始時(shí)間、多普勒多項(xiàng)式、近斜距、雷達(dá)波長、13個(gè)衛(wèi)星位置矢量和速度矢量、方位向分辨率、距離向分辨率、基線矢量、地球長半軸、地球短半軸)、步驟S2.3擬合衛(wèi)星位置矢量與所述衛(wèi)星成像時(shí)刻t之間的二階多項(xiàng)式、步驟S2.4擬合衛(wèi)星速度矢量與所述衛(wèi)星成像時(shí)刻t之間的二階多項(xiàng)式，經(jīng)計(jì)算

根據(jù)VNW坐標(biāo)系定義可知，W軸指向下方時(shí)，取正號；W軸指向上方時(shí)，取負(fù)號，這里指向上方或下方，并非嚴(yán)格的指向天頂或地心方向。優(yōu)選的，在本實(shí)施例中取正號。

步驟S4.3.5，在步驟S2.2已經(jīng)說明了SAR嚴(yán)密成像模型是在地心空間直角坐標(biāo)系下完成的，提供的目標(biāo)點(diǎn)P的地理坐標(biāo)也是在地心空間直角坐標(biāo)系下表示，為了得到地心空間直角坐標(biāo)系下的目標(biāo)點(diǎn)P的單位視向量，需要將VNW移動(dòng)坐標(biāo)系下的角反射器點(diǎn)坐標(biāo)轉(zhuǎn)換到地心空間直角坐標(biāo)系，其轉(zhuǎn)換矩陣為：

其中，b為干涉基線矢量，v為衛(wèi)星飛行速度矢量，為衛(wèi)星飛行速度的單位矢量。

優(yōu)選的，對通過步驟S1獲得的成像參數(shù)(包括基線矢量、衛(wèi)星速度矢量)、步驟S2獲得的衛(wèi)星速度矢量，經(jīng)計(jì)算有：

步驟S4.4，基于步驟S4.2得到的主影像對應(yīng)的衛(wèi)星位置矢量S₁，步驟S4.3得到的和推導(dǎo)出所述三維重建模型為：

其中，P為目標(biāo)點(diǎn)坐標(biāo)，S₁為衛(wèi)星位置矢量，r₁為視向量，r_s1表示主影像對應(yīng)的衛(wèi)星到地面目標(biāo)點(diǎn)的距離，為轉(zhuǎn)換矩陣，為單位視向量，r_v為單位視向量在V軸上的分量，r_n為單位視向量在N軸上的分量，r_w為單位視向量在W軸上的分量，b為干涉基線矢量，v為載機(jī)飛行速度矢量，為飛行速度的單位矢量，λ為雷達(dá)波長，為絕對干涉相位，Q為天線工作模式，|b_pv|為垂直于衛(wèi)星速度方向的基線矢量模長，b_v為基線矢量在衛(wèi)星速度方向的投影值，f_dop為多普勒中心頻率，B為干涉基線模長。

以角反射器點(diǎn)D002為例說明計(jì)算的過程，主影像對應(yīng)的衛(wèi)星位置矢量斜距r_s1＝5421.521，旋轉(zhuǎn)矩陣以及單位視向量有：

步驟S4.2-S4.4對步驟S4.1提供的基本表達(dá)式進(jìn)行了公式推導(dǎo)和描述，確定了三維重建模型的解析表達(dá)形式，為后續(xù)步驟S4.5-S4.8提供了理論參考依據(jù)，將根據(jù)步驟S4.2-S4.4確定的三維重建模型，依據(jù)該三維重建模型對干涉系統(tǒng)參數(shù)進(jìn)行偏導(dǎo)數(shù)的計(jì)算，干涉測量檢校模型(誤差模型)的構(gòu)建以及最小二乘求解等步驟。下面具體描述步驟S4.5-S4.8的處理過程。

步驟S4.5，對步驟S4.4的三維重建模型進(jìn)行改寫，得到三維重建模型的觀測方程；

在S4.1-S4.4建立的三維重建模型過程中，各項(xiàng)系統(tǒng)參數(shù)不可避免的會(huì)受誤差的影響，為了獲得精確的各項(xiàng)系統(tǒng)參數(shù)，需要利用角反射器點(diǎn)進(jìn)行檢校處理，優(yōu)化干涉系統(tǒng)參數(shù)，輸出高精度的干涉系統(tǒng)參數(shù)，為最終獲取高精度的角反射器坐標(biāo)提供可靠干涉系統(tǒng)參數(shù)。干涉測量檢校模型以三維重建模型為基礎(chǔ)，推導(dǎo)各系統(tǒng)參數(shù)的偏導(dǎo)數(shù)，進(jìn)而利用最小二乘原理求解系統(tǒng)參數(shù)，優(yōu)化和提高角反射器點(diǎn)的像素坐標(biāo)，并剔除粗差。具體地，對公式(7)進(jìn)行改寫，得到其函數(shù)表達(dá)形式，即三維重建模型的觀測方程：

式中F_x(X)、F_y(X)、F_z(X)分別為X、Y、Z方向的函數(shù)，為待優(yōu)化的干涉系統(tǒng)參數(shù)，b_l為基線長度、b_α為基線傾角、r₀為近斜距、為干涉相位，f_d為多普勒頻率，這些參數(shù)均參與三維重建模型的計(jì)算，而且是存在一定誤差的，需要對這些參數(shù)進(jìn)行優(yōu)化處理；S_1x,S_1y,S_1z為主影像對應(yīng)的衛(wèi)星位置坐標(biāo)，為主影像到地面角反射器(即地面目標(biāo)點(diǎn))的斜距，u₁₁,u₁₂,u₁₃,u₂₁,u₂₂,u₂₃,u₃₁,u₃₂,u₃₃為旋轉(zhuǎn)矩陣的計(jì)算變量；r_v為單位視向量在V軸上的分量，r_n為單位視向量在N軸上的分量，r_w為單位視向量在W軸上的分量；P_x,P_y,P_z為角反射器的地心空間直角坐標(biāo)，均通過步驟S1提供的主、副影像參數(shù)信息獲得。

步驟S4.6，計(jì)算所述三維重建模型的觀測方程中的干涉系統(tǒng)參數(shù)的偏導(dǎo)數(shù)；

根據(jù)步驟S4.5的公式(8)，公式(8)對干涉系統(tǒng)參數(shù)(包括b_l為基線長度、b_α為基線傾角、r₀為近斜距、為干涉相位，f_d為多普勒頻率)求導(dǎo)進(jìn)行線性化得到誤差方程(即，干涉測量檢校模型)，因此需要計(jì)算干涉系統(tǒng)參數(shù)的偏導(dǎo)數(shù)，其具體包括以下子步驟：

步驟S4.6.1，根據(jù)式(8)，分別計(jì)算角反射器點(diǎn)在X、Y、Z三個(gè)方向上F_x、F_y和F_z函數(shù)對基線長度b_l的偏導(dǎo)數(shù)。我們以F_x對基線長度b_l求偏導(dǎo)數(shù)為例說明求偏導(dǎo)數(shù)的過程，經(jīng)計(jì)算有其中為主影像到地面角反射器(即地面目標(biāo)點(diǎn))的斜距(距離)，旋轉(zhuǎn)矩陣單位視向量r_v為單位視向量在V軸上的分量，r_n為單位視向量在N軸上的分量，r_w為單位視向量在W軸上的分量，經(jīng)分析發(fā)現(xiàn)：和u₁₁,u₁₂,u₁₃的計(jì)算過程不涉及到基線長度b_l，然而的計(jì)算需要用到基線長度b_l，所以和u₁₁,u₁₂,u₁₃對基線長度b_l的偏導(dǎo)數(shù)為0，單位視向量對基線長度b_l的偏導(dǎo)數(shù)不為0，是需要計(jì)算的，最后F_x對基線長度b_l求偏導(dǎo)數(shù)就轉(zhuǎn)換為單位視向量對基線長度b_l求偏導(dǎo)數(shù)。F_y和F_z對基線長度b_l求偏導(dǎo)數(shù)與F_x對基線長度b_l求偏導(dǎo)數(shù)類似。因此，在X、Y、Z三個(gè)方向上F_x、F_y和F_z函數(shù)對基線長度b_l求偏導(dǎo)數(shù)轉(zhuǎn)化為單位視向量對基線長度b_l求導(dǎo)數(shù)，得到：

其中，b_l為基線長度；單位視向量r_v為單位視向量在V軸上的分量，r_n為單位視向量在N軸上的分量，r_w為單位視向量在W軸上的分量；為主影像到地面角反射器點(diǎn)的斜距，Δr為主影像對應(yīng)的衛(wèi)星到角反射器點(diǎn)的距離與主影像對應(yīng)的衛(wèi)星到角反射器點(diǎn)的距離之差。

步驟S4.6.2，根據(jù)式(8)，分別計(jì)算X、Y、Z三個(gè)方向上F_x、F_y和F_z函數(shù)對基線傾角b_α的偏導(dǎo)數(shù)，在移動(dòng)坐標(biāo)系中單位視向量對基線傾角b_α的導(dǎo)數(shù)都是0，基線傾角b_α僅在建立移動(dòng)坐標(biāo)系時(shí)對其產(chǎn)生影響，對b_α的導(dǎo)數(shù)為：

其中，b_α為基線傾角，單位視向量r_v為單位視向量在V軸上的分量，r_n為單位視向量在N軸上的分量，r_w為單位視向量在W軸上的分量；，為主影像到地面角反射器的斜距，為旋轉(zhuǎn)矩陣。

式(11)中對b_α求偏導(dǎo)數(shù)，有：

其中，和為VNW移動(dòng)坐標(biāo)系下的三個(gè)正交基，b_α為基線傾角，為旋轉(zhuǎn)矩陣。

步驟S4.6.3，根據(jù)式(8)，分別計(jì)算X、Y、Z三個(gè)方向上F_x、F_y和F_z函數(shù)對初始斜距的偏導(dǎo)數(shù)。主天線到地面角反射器點(diǎn)的斜距與初始斜距r₀之間的關(guān)系：沿距離向上地面點(diǎn)的影像坐標(biāo)為n_rg，斜距上的像元分辨率為r_ps。式(9)對r₀求導(dǎo)數(shù)，得：

其中，F(xiàn)_x、F_y和F_z分別為X、Y、Z方向的函數(shù)，單位視向量r_v為單位視向量在V軸上的分量，r_n為單位視向量在N軸上的分量，r_w為單位視向量在W軸上的分量，u₁₁,u₁₂,u₁₃,u₂₁,u₂₂,u₂₃,u₃₁,u₃₂,u₃₃為旋轉(zhuǎn)矩陣的計(jì)算變量，r₀為近斜距，為主天線到地面角反射器點(diǎn)的斜距。

再求單位視向量對r₀的導(dǎo)數(shù)：

上式中，基線矢量與衛(wèi)星速度方向垂直的分量將式(14)代入式(13)中求r₀的偏導(dǎo)數(shù)。其中，為主天線到地面角反射器點(diǎn)的斜距，單位視向量r_v為單位視向量在V軸上的分量，r_n為單位視向量在N軸上的分量，r_w為單位視向量在W軸上的分量，b_pv為基線矢量與衛(wèi)星速度方向垂直的分量，Δr為主副影像衛(wèi)星位置到角反射器點(diǎn)的距離之差，b_l為基線長度。

步驟S4.6.4，根據(jù)式(8)，分別計(jì)算X、Y、Z三個(gè)方向上F_x、F_y和F_z函數(shù)對干涉相位的偏導(dǎo)數(shù)，得：

其中，單位視向量對的導(dǎo)數(shù)為：

其中，為干涉相位，為主天線到地面角反射器點(diǎn)的斜距，單位視向量r_v為單位視向量在V軸上的分量，r_n為單位視向量在N軸上的分量，r_w為單位視向量在W軸上的分量，b_pv為基線矢量與衛(wèi)星速度方向垂直的分量，為絕對干涉相位，b_l為基線長度；Q為天線工作模式，其中，當(dāng)天線的工作模式為重復(fù)軌道模式時(shí)，Q＝2，當(dāng)天線的工作模式為一發(fā)雙收模式時(shí)，Q＝1；b_pv為基線矢量與衛(wèi)星速度方向垂直的分量，λ為雷達(dá)波長。

步驟S4.6.5，根據(jù)式(8)，分別計(jì)算X、Y、Z三個(gè)方向上F_x、F_y和F_z函數(shù)對多普勒中心頻率f_d的偏導(dǎo)數(shù)，得：

其中，單位視向量對多普勒中心頻率f_d的導(dǎo)數(shù)為：

式中，f_d為多普勒中心頻率，b_v為干涉基線與衛(wèi)星速度方向平行的投影分量；單位視向量r_v為單位視向量在V軸上的分量，r_n為單位視向量在N軸上的分量，r_w為單位視向量在W軸上的分量，b_pv為基線矢量與衛(wèi)星速度方向垂直的分量，λ為雷達(dá)波長，|v|為衛(wèi)星速度矢量的模長。

步驟S4.7，基于步驟S4.5得到的觀測方程和S4.6推導(dǎo)的干涉系統(tǒng)參數(shù)的偏導(dǎo)數(shù)，利用若干個(gè)角反射器點(diǎn)對干涉系統(tǒng)參數(shù)進(jìn)行非線性最小二乘估計(jì)；

利用角反射器點(diǎn)構(gòu)建干涉測量檢校模型。假設(shè)在SAR影像中存在m_G個(gè)角反射器點(diǎn)，在各角反射器點(diǎn)處有誤差方程(即，干涉測量檢校模型)：

式中，v_xG(i),v_yG(i),v_zG(i)分別為角反射器點(diǎn)G_x(i),G_y(i),G_z(i)在X、Y、Z方向的殘差，分別為F_xG(i),F_yG(i),F_zG(i)函數(shù)對基線長度b_l在X、Y、Z方向的偏導(dǎo)數(shù)；分別為F_xG(i),F_yG(i),F_zG(i)函數(shù)對基線傾角b_α在X、Y、Z方向的偏導(dǎo)數(shù)；分別為F_xG(i),F_yG(i),F_zG(i)函數(shù)對近斜距r₀在X、Y、Z方向的偏導(dǎo)數(shù)；分別為F_xG(i),F_yG(i),F_zG(i)函數(shù)對干涉相位在X、Y、Z方向的偏導(dǎo)數(shù)；分別為F_xG(i),F_yG(i),F_zG(i)函數(shù)對多普勒頻率f_dop在X、Y、Z方向的偏導(dǎo)數(shù)；l_xG(i),l_yG(i),l_zG(i)為角反射點(diǎn)已知地心坐標(biāo)與三維重建模型解算出的地心坐標(biāo)之差構(gòu)成的常數(shù)項(xiàng)，角反射器點(diǎn)G(i)表示第i個(gè)角反射器點(diǎn)，i＝1,2,…,n_G，總計(jì)有n_G個(gè)角反射器點(diǎn)。

步驟S4.8，基于最小二乘原理對干涉系統(tǒng)參數(shù)進(jìn)行解算，其具體包括以下子步驟：

步驟S4.8.1，根據(jù)式(18)所列出的誤差模型，列出誤差方程式

V＝AΔx-l (19)

其中，

矩陣V為殘差項(xiàng)矩陣，其中v_xG(1)，v_yG(1)，v_zG(1)表示第G(1)個(gè)角反射器點(diǎn)在X、Y、Z方向的殘差，表示第G(m_G)個(gè)角反射器點(diǎn)在X、Y、Z方向的殘差；矩陣A為偏導(dǎo)數(shù)系數(shù)矩陣，其中分別為F_xG(i),F_yG(i),F_zG(i)函數(shù)對基線長度b_l在X、Y、Z方向的偏導(dǎo)數(shù)；分別為F_xG(i),F_yG(i),F_zG(i)函數(shù)對基線傾角b_α在X、Y、Z方向的偏導(dǎo)數(shù)；分別為F_xG(i),F_yG(i),F_zG(i)函數(shù)對近斜距r₀在X、Y、Z方向的偏導(dǎo)數(shù)；分別為F_xG(i),F_yG(i),F_zG(i)函數(shù)對干涉相位在X、Y、Z方向的偏導(dǎo)數(shù)；分別為F_xG(i),F_yG(i),F_zG(i)函數(shù)對多普勒頻率f_dop在X、Y、Z方向的偏導(dǎo)數(shù)，i＝1,2,…,n_G；矩陣l為常數(shù)項(xiàng)矩陣l_xG(i),l_yG(i),l_zG(i)分別為L_xG(i)-L_0,xG(i)，L_yG(i)-L_0,yG(i)，L_zG(i)-L_0,zG(i)計(jì)算得到，i＝1,2,…,n_G為計(jì)算的第i個(gè)角反射器點(diǎn)，n_G為角反射器點(diǎn)總個(gè)數(shù)；Δx為待求解的干涉系統(tǒng)參數(shù)改正量矩陣，基線長度改正量Δb_l、基線傾角改正量Δb_α、近斜距改正量Δr₀、干涉相位改正量和多普勒頻率改正量Δf_d。

步驟S4.8.2，通過步驟S1獲取的影像成像參數(shù)(表1所示)，計(jì)算系數(shù)矩陣A，在步驟S4.7和S4.8中已描述其計(jì)算公式，帶入計(jì)算即可；通過步驟S2獲取的角反射點(diǎn)地心空間直角坐標(biāo)，計(jì)算常數(shù)項(xiàng)l＝L₀-L，L由三維重建模型公式(7)計(jì)算得到，L₀由野外測量人員實(shí)測角反射器點(diǎn)位置得到。

計(jì)算系數(shù)矩陣A，其計(jì)算公式如下式所示：

我們以角反射器點(diǎn)D002為例，分別計(jì)算D002點(diǎn)在X、Y、Z三個(gè)方向上F_x、F_y和F_z函數(shù)對基線長度b_l的偏導(dǎo)數(shù)，說明具體的計(jì)算過程。在步驟S4.5的公式(8)中描述了三維重建模型的函數(shù)表達(dá)形式，如下式所示：

計(jì)算D002點(diǎn)在X、Y、Z三個(gè)方向上F_x、F_y和F_z函數(shù)對基線長度b_l的偏導(dǎo)數(shù)，有涉及到主影像到地面角反射器的斜距旋轉(zhuǎn)矩陣單位視向量對基線長度b_l的偏倒數(shù)

計(jì)算根據(jù)步驟S2.5可知，根據(jù)表4中提供的D002點(diǎn)像素坐標(biāo)(x,y)＝(1000.5,1638.0625)；步驟S提供的表1中近斜距R₀＝4885.9532m，距離向采樣間隔ΔR＝0.536m；最后得到

計(jì)算旋轉(zhuǎn)矩陣根據(jù)步驟S4.3.5可知，結(jié)合步驟S1提供的速度矢量和基線矢量，得

計(jì)算單位視向量對基線長度b_l的偏倒數(shù)的計(jì)算，可根據(jù)步驟S4.3.1，S4.3.2和S4.3.3得到的數(shù)值，以及b_l＝0.312816m，可得

針對基線傾角b_α、近斜距r₀、干涉相位多普勒頻率f_dop分別重復(fù)上述步驟，即可得到系數(shù)矩陣A。下面列出D002點(diǎn)在X、Y、Z三個(gè)方向上F_x、F_y和F_z函數(shù)對基線長度b_l、基線傾角b_α、近斜距r₀、干涉相位多普勒頻率f_dop的偏導(dǎo)數(shù)，如下所示：

計(jì)算常數(shù)項(xiàng)矩陣l，其計(jì)算公式如下式所示：

我們以角反射器點(diǎn)D002為例，有其中L_x,L_y,L_z由步驟S4.4的三維重建模型公式(7)計(jì)算得到；L_0,x,L_0,y,L_0,z由表2提供的經(jīng)度、緯度、高程即(110.0559,35.10238,667.631)，經(jīng)步驟S2.1將角反射器的控制點(diǎn)的大地坐標(biāo)轉(zhuǎn)換為地心空間直角坐標(biāo)；最終得到

l_x,D002＝L_0,x-L_x＝-1791657.022-(-1791633.082)＝-23.940

l_y,D002＝L_0,y-L_y＝4907631.293-4907598.887＝32.406。

l_z,D002＝L_0,z-L_z＝3647548.734-3647534.327＝14.406

步驟S4.8.3，平差準(zhǔn)則選擇常用的最小二乘準(zhǔn)則：V^TV＝min；求解干涉系統(tǒng)參數(shù)的改正數(shù)，有：

Δx＝(A^TA)^-1A^Tl (20)

其中，Δx為待求解的干涉系統(tǒng)參數(shù)改正量矩陣，A為步驟S4.8.2計(jì)算得到的系數(shù)矩陣，A^T為矩陣A的轉(zhuǎn)置矩陣，l為常數(shù)項(xiàng)矩陣。

步驟S4.8.4，計(jì)算每次迭代前后角反射器點(diǎn)的坐標(biāo)差的均方根，據(jù)此判斷是否進(jìn)行下一步的迭代計(jì)算，公式為：

式中，eps為每次迭代前后角反射器點(diǎn)的坐標(biāo)差的均方根，表示第i個(gè)角反射器點(diǎn)根據(jù)步驟S4.8構(gòu)建的誤差模型優(yōu)化后的第k次迭代計(jì)算出的坐標(biāo)，表示第i個(gè)角反射器點(diǎn)根據(jù)步驟S4.8構(gòu)建的誤差模型優(yōu)化后的第k-1次迭代計(jì)算出的坐標(biāo)，n_G為角反射器點(diǎn)總個(gè)數(shù)，當(dāng)eps小于設(shè)定的閾值時(shí)平差結(jié)束，統(tǒng)計(jì)各角反射器點(diǎn)平差前后的殘差。若平差過程中，異常終止，則需檢查相關(guān)參數(shù)，重新計(jì)算，直至滿足要求。

優(yōu)選的，eps取10^-4。

在步驟S4.5-S4.8中，依據(jù)該三維重建模型對干涉系統(tǒng)參數(shù)進(jìn)行偏導(dǎo)數(shù)的計(jì)算，干涉測量檢校模型(誤差模型)的構(gòu)建以及最小二乘求解等步驟，完整描述了誤差模型的推導(dǎo)、構(gòu)建、求解、輸出結(jié)果等處理過程，為后續(xù)步驟S4.9-S4.12提供了角反射器點(diǎn)坐標(biāo)。而在步驟S4.9-S4.12中，將對角反射器點(diǎn)的像素精度進(jìn)行評價(jià)并剔除粗差，輸出最終的角反射器點(diǎn)子像素級坐標(biāo)。

步驟S4.9，利用步驟S4.8獲得的干涉系統(tǒng)參數(shù)，對步驟S3中的角反射器像點(diǎn)坐標(biāo)進(jìn)行干涉參數(shù)測量校驗(yàn)，剔除誤差較大的角反射器點(diǎn)；

由步驟S4.8可知，步驟S3獲取的像素級影像坐標(biāo){P_pix}，由于SAR影像易受斑點(diǎn)噪聲的影響，肯定存在誤差較大的點(diǎn)，這些點(diǎn)會(huì)影響干涉測量檢校的精度，據(jù)此可以根據(jù)干涉測量檢校模型進(jìn)行粗差探測和剔除，自動(dòng)探測角反射器的粗差點(diǎn)，迭代計(jì)算，剔除粗差點(diǎn)，實(shí)現(xiàn)子像素級自動(dòng)定位，其具體包括以下子步驟：

步驟S4.9.1，步驟S3獲取的像素級影像坐標(biāo)集合{P_pix}中所有角反射器像點(diǎn)坐標(biāo)均參與干涉測量檢校平差計(jì)算，直至滿足最小二乘約束條件，迭代結(jié)束，此時(shí)將通過模型反演的干涉系統(tǒng)參數(shù)，根據(jù)步驟S4.4的公式(7)計(jì)算角反射器點(diǎn)對應(yīng)的地心空間三維坐標(biāo)，結(jié)果記為{P_{XYZ_m}}。

計(jì)算角反射器點(diǎn)對應(yīng)的地心空間三維坐標(biāo)基本公式如下式所示：

那么根據(jù)優(yōu)化后的干涉系統(tǒng)參數(shù)(包括基線長度b_l、基線傾角b_α、近斜距r₀、干涉相位和多普勒頻率f_d)和步驟S1提供的主、副影像參數(shù)列表(包括基線矢量、13個(gè)衛(wèi)星位置矢量和速度矢量、雷達(dá)波長、近斜距)以及步驟S3.8提供的表4角反射器點(diǎn)像素坐標(biāo)。逐一將角反射器點(diǎn)帶入公式計(jì)算，結(jié)果記為{P_{XYZ_m}}。

步驟S4.9.2，將通過步驟S2.1獲得的角反射點(diǎn)已知地心空間直角坐標(biāo)集合為{P_XYZ}，計(jì)算兩集合{P_XYZ}與{P_{XYZ_m}}在X方向與Y方向的平面點(diǎn)位中誤差{P_{XYZ_err}}，并獲取平面點(diǎn)位坐標(biāo)誤差的均值和標(biāo)準(zhǔn)差μ、σ，有如下計(jì)算公式：

P_{XY_err}＝SQRT((P_{XYZ_m}(X)-P_XYZ(X))²+(P_{XYZ_m}(Y)-P_XYZ(Y))²)

其中，P_{XYZ_err}為X方向與Y方向的平面點(diǎn)位中誤差，SQRT()為開根號運(yùn)算，P_{XYZ_m}(X)為步驟S4.9.1得到的地心空間三維坐標(biāo){P_{XYZ_m}}的X方向，P_XYZ(X)為外業(yè)實(shí)測的角反射器點(diǎn)的地心空間三維坐標(biāo){P_XYZ}的X方向，P_{XYZ_m}(Y)為步驟S4.9.1得到的地心空間三維坐標(biāo){P_{XYZ_m}}的Y方向，P_XYZ(Y)為外業(yè)實(shí)測的角反射器點(diǎn)的地心空間三維坐標(biāo){P_XYZ}的Y方向。μ為平面點(diǎn)位坐標(biāo)誤差的均值，P_{XY_err}(i)為{P_{XYZ_err}}集合里的元素，n_G為角反射器點(diǎn)總個(gè)數(shù)。σ為平面點(diǎn)位坐標(biāo)誤差的標(biāo)準(zhǔn)差。

步驟S4.9.3，根據(jù)步驟S1獲得的影像的距離向、方位向分辨率為(x_Res,y_Res)，可以得到對應(yīng)的像點(diǎn)地面分辨率若步驟S4.9.2的平面點(diǎn)位的標(biāo)準(zhǔn)差σ大于像點(diǎn)地面分辨率(xy)_Res，則說明存在誤差較大的角反射器點(diǎn)。

優(yōu)選的，(xy)_Res＝0.75m。

步驟S4.9.4，計(jì)算平面點(diǎn)位的標(biāo)準(zhǔn)差σ與像點(diǎn)地面分辨率k*(xy)_Res之差，考慮到角反射器點(diǎn)的定位精度要滿足子像素級的要求，必須優(yōu)于一個(gè)像素或者一個(gè)像素對應(yīng)的地面分辨單位大小，k表示優(yōu)于一個(gè)地面分辨單元的比例系數(shù)k∈(0.1,1)。按照1倍標(biāo)準(zhǔn)差σ閾值剔除粗差，如下公式所示，若σ＞k*(xy)_Res，則保留{e}中數(shù)值范圍滿足下列條件的點(diǎn)，不滿足條件的點(diǎn)自動(dòng)設(shè)置為無效點(diǎn)，不再參與平差計(jì)算，直至參與平差的所有角反射器點(diǎn)滿足設(shè)定的閾值條件；

e∈(μ-σ,μ+σ)

優(yōu)選的，k＝0.5，k*(xy)_Res＝0.375m。

步驟S4.9.5，經(jīng)過步驟S4.7-S4.9.4進(jìn)行干涉測量檢校平差和粗差剔除運(yùn)算，經(jīng)過三次終止運(yùn)算，滿足子像素級定位精度要求，獲得模型優(yōu)化參數(shù)(即，干涉系統(tǒng)參數(shù))，其統(tǒng)計(jì)結(jié)果如下表5所示：

表5統(tǒng)計(jì)模型運(yùn)算結(jié)果

步驟S4.10，利用步驟S4.9.5獲取的模型優(yōu)化參數(shù)(包括基線長度b_l、基線傾角b_α、近斜距r₀、干涉相位和多普勒頻率)，根據(jù)前述步驟S2構(gòu)建的SAR嚴(yán)密成像模型重新計(jì)算角反射器點(diǎn)的影像像點(diǎn)坐標(biāo)，最終獲得的子像素級角反射器點(diǎn)坐標(biāo)為(x₂,y₂)，并將所有坐標(biāo)集合記為{P_{pix_f}}見表6。

表6幾何關(guān)系優(yōu)化后的角反射器點(diǎn)影像坐標(biāo)

以上內(nèi)容僅為本發(fā)明的較佳實(shí)施例，對于本領(lǐng)域的普通技術(shù)人員，依據(jù)本發(fā)明的思想，在具體實(shí)施方式及應(yīng)用范圍上均會(huì)有改變之處，本說明書內(nèi)容不應(yīng)理解為對本發(fā)明的限制。