CN116203562A - 一种光学影像辅助的星载立体sar影像控制点自动生成方法 - Google Patents

Abstract

本发明涉及一种光学影像辅助的星载立体SAR影像控制点自动生成方法，与现有技术相比实现了从不同升降轨立体SAR影像中自动化的完成灯柱这一类控制点的匹配、提取、误差剔除，并形成可用的控制点。本发明包括以下步骤：基于高分光学影像上提取立体SAR影像上同名灯柱点的近似坐标；SAR影像灯柱PS点像方坐标粗略提取；SAR影像灯柱PS点坐标精确提取；SAR影像灯柱PS点异常值剔除；SAR影像灯柱PS点坐标误差补偿；立体平差处理生成灯柱控制点三维坐标。本发明实现立体SAR影像同名PS点自动化匹配、提取、误差剔除，生成高精度可用的控制点，为基础测绘以及其他卫星几何校正提供控制数据。

Description

一种光学影像辅助的星载立体SAR影像控制点自动生成方法

技术领域

本发明属于星载合成孔径雷达影像地面应用领域，具体来说是一种光学影像辅助的星载立体SAR影像控制点自动生成方法。

背景技术

星载合成孔径雷达SAR(Synthetic Aperture Radar)由于不受云雨影响，可全天时、全天候工作，已成为一种不可或缺的空间对地观测手段。星载SAR相对光学卫星几何定位不受姿态的影响，全球无控几何定位精度较好。国外TerraSAR卫星经几何定标后，绝对定位精度优于0.3m,国产SAR卫星经几何定标后，绝对定位精度可优于1.5m。因此，可应用星载立体SAR生成高精度控制点，服务于基础测绘。

针对利用星载立体SAR生成地面控制点研究，目前研究主要集中在地面控制点选择和SAR影像匹配两方面。在控制点选取方面，有研究将立体SAR技术应用于城市区域中少量天然永久散射体PS(Permanent Scatter)点，这些PS点主要来自建筑物的外墙立面上的窗角或路灯的底部，所有参与立体平差的同名点影像坐标均通过人工在立体SAR影像中提取。另外也有研究尝试自动化提取和匹配有限数量的建筑物外墙立面上的窗角点，但方法只适用于相同升降轨且入射角差异不大的立体SAR影像。在SAR影像匹配方面，有研究提出SAR-SIFT方法，该方法通过修改SIFT(Scale-invariant Feature Transform)算法以考虑斑点的统计特性，从而在SAR影像中有效提取局部描述符，但是该方法仅适用于从相同升降轨的SAR影像并且要求入射角差异较小。另外，有研究提出基于Harris的局部特征来检测相同的散射体点，通过使用外部数字高程模型DEM(Digital Elevation Model)和轨道信息对影像进行地理编码来限制搜索空间，并最终使用SIFT进行匹配，尽管该方法在不同升降轨的立体SAR影像上有一定适用性，但仅适用于孤立的PS点。

现有研究中，从不同观测几何的立体SAR影像中识别和提取同名PS点，目前多为手动提取，应发展一种自动化的提取过程。另外采用不同升降轨的SAR立体影像配置可获得大的基高比，从而提升立体SAR精度，但要求能够稳健的从不同升降轨配置的立体SAR影像中检测和识别同名PS点。国内暂无对立体SAR影像自动生成高精度控制点的论文与专利。

发明内容

本发明的目的是为了解决现有技术中无法对具有不同升降轨的立体SAR影像自动化处理生成可用的控制点，提供一种光学影像辅助的星载立体SAR影像控制点自动生成方法来解决上述问题。

为了实现上述目的，本发明技术方案如下：

一种光学影像辅助的星载立体SAR影像控制点自动生成方法，包括以下步骤：

11)基于高分光学影像上提取立体SAR影像上同名灯柱点的近似坐标：获取立体SAR影像覆盖区域所对应的高分光学遥感影像I₀，基于I₀检测到灯柱PS点后并得到其三维坐标，利用SAR影像的距离多普勒几何定位模型将灯柱的三维坐标反算至SAR影像像方，完成立体SAR影像上同名灯柱PS点的近似像点坐标(τ₀,t₀)提取。

12)SAR影像灯柱PS点像方坐标粗略提取：对SAR影像进行阈值化处理，得到二值图，将光学影像中的路灯点坐标利用SAR的几何定位模型反算至其像方，应用迭代最近点法实现灯柱点集与SAR阈值处理后的亮点点集匹配，从而实现SAR影像灯柱PS点的匹配，并得到灯柱点对应SAR影像上的粗略坐标。

13)SAR影像灯柱PS点坐标精确提取：以近似坐标(x_0,pi,SAR,y_0,pi,SAR)为中心截取SAR影像块，进行点目标分析，精确提取灯柱p_i点的方位和距离坐标(x_pi,SAR,y_SAR)。

14)SAR影像灯柱PS点异常值剔除：潜在灯柱PS点候选对象SCR应足够高，但由于存在误差，进一步通过计算灯柱PI点的相位噪声σ_φ，进行异常值剔除。

15)SAR影像灯柱PS点坐标误差补偿：对提取的灯柱点p_i影像坐标(τ_pi,SAR,t_pi,SAR)进行系统定位误差和大气延迟误差的校正，得到改正后的时间量纲坐标(τ_pi,cor,SAR,t_pi,cor,SAR)。

16)立体平差处理生成灯柱控制点三维坐标：针对待生成控制点区域内可用的N张(N＞2)立体SAR影像，按照上述步骤11)-15)提取同名灯柱点对集合，然后进行立体平差处理，生成灯柱控制点的三维坐标。

所述基于高分光学影像上提取灯柱点包括以下步骤：

21)获取对应区域的光学影像I₀后，为了使灯柱的阴影更加突出，对光学影像I₀进行高频提升滤波处理：

I＝I₀+aI_m，

式中，I₀为原始图像，a是锐化因子标量，I_m是使用非锐化蒙板后的图像，为I₀与其对应模糊图像之差，其中a值越高，表示锐化程度越高。

22)在锐化后的光学图像I中提取任意包含灯柱阴影的模板图像T，模板大小设置为256×256，然后将模板与图像I进行相关运算，得到图像I在像方(u,v)坐标处的相似度：

式中，I(x,y)和T(x,y)分别表示图像I和模板图像T在(x,y)坐标处的像素值，N1,N2表示模板图像T的大小，

和/>

分别表示图像I和模板图像T的平均像素值，ρ(u,v)为归一化互相关值，ρ(u,v)不受图像亮度或对比度的变化影响，因此可提升模板匹配的效果。

23)光学图像I中经模板匹配后得到每个灯柱对象存在多个像素，进一步对模板提取到的像素进行聚类处理，得到表示灯柱唯一的像素点p_i，处理如下：

式中，p_i表示计算m(p_i)漂移向量的空间坐标点，p_j代表p_i附近的点，g为具有带宽h的高斯核函数，而|| ||表示欧几里得距离算子。

24)根据光学图像I上灯柱聚类点p_i的像素坐标(x_pi,I,y_pi,I)，获取其平面地理坐标(Lat_pi,I,Lon_pi,I)，并利用平面坐标在开源SRTM DEM数据上内插出对应高程H_pi,I，构成完整的三维坐标(Lat_pi,I,Lon_pi,I,H_pi,I)。

所述SAR影像匹配灯柱PS点包括以下步骤：

31)根据立体SAR影像的元数据信息构建距离多普勒RD几何定位方程：

式中，f_R、f_A分别表示距离方程和多普勒方程，τ、t分别表示距离向时间和方位向时间，C_S(t)和

分别表示成像t时刻SAR天线相位中心在WGS84坐标系下的位置矢量和速度矢量，C_T为观测目标在WGS84坐标系下的位置矢量，fdc(τ)表示在距离双程延时为τ处的多普中心频率；

32)逐个将灯柱聚类点p_i的三维坐标(Lat_pi,I,Lon_pi,I,H_pi,I)利用RD定位几何方程反投影至SAR影像的像方，得到灯柱点集。

33)根据一定阈值，对SAR影像上亮点进行二值掩膜处理。

34)利用迭代最近点法，实现反算在SAR像方的灯柱点集与SAR影像上亮点集匹配，SAR影像上只保留匹配到路灯点p_i的亮点，从而完成SAR影像灯柱点的匹配，并得到粗略匹配坐标(x_0,pi,SAR,y_0,pi,SAR)。

所述SAR影像灯柱PS点坐标精确提取包括以下步骤：

41)以近似坐标(x_0,pi,SAR,y_0,pi,SAR)为中心在SAR影像上截取32*32的影像块，并寻找最大峰值位置；

42)以最大峰值处为中心截取32*32样本窗做二维离散傅里叶变换，32倍采样因子频域补零后进行二维逆离散傅里叶变换，得到升采样后的点目标影像块；

43)对过采样后的影像块，以最大峰值处截取3×3影像子块，进行二维抛物面插值，椭圆抛物面插值模型如下：

f(x,y)＝a₅x²+a₄y²+a₃xy+a₂x+a₁y+a₀，

式中x,y为在3×3影像子块中的像素坐标，为a_i(1≤i≤6)为抛物面插值模型f(x,y)的参数，将3×3子块像素值代入上式中并通过最小二乘求解a_i。参数a_i求解后，进而获得灯柱p_i点精确坐标(x_pi,SAR,y_SAR)。

所述SAR影像灯柱PS点异常值剔除包括以下步骤：

51)计算灯柱p_i点的峰值功率P_peak和杂波功率P_clutter。峰值功率P_peak是以峰值(x_pi,SAR,y_SAR)为中心为2×2分辨率像元的主瓣区域的积分信号能量，杂波功率P_clutter是主瓣和旁瓣区域外的积分信号能量。

52)计算灯柱p_i点的信杂比SCR并得到相位噪声

计算公式如下：

/>

53)如果灯柱p_i点计算出的σ_φ大于0.5rad，则予以剔除，否则保留。并将灯柱影像p_i在SAR影像上的(x_pi,SAR,y_SAR)坐标转为时间量纲的坐标(τ_pi,SAR,t_pi,SAR)。转化公式如下：

式中，x、y为待转换的像素坐标，τ、t为转换后的时间量纲的坐标，分别表示距离向双程时延和方位向时间，PRF为方位向脉冲重复频率，F_s为距离向采样频率，τ_min为SAR影像近距端双程时延，t_start为方位向起始时间。

所述SAR影像灯柱PS点坐标误差补偿包括以下步骤：

61)从SAR影像元数据中读取系统误差定标参数(Δτ^cal、Δt^cal)，补偿至灯柱p_i点影像坐标中，公式如下：

上式中，τ_pi,1和t_pi,1分别为对灯柱p_i点影像初始时间量纲坐标(τ_pi,SAR,t_pi,SAR)进行系统误差补偿后的距离向和方位向时间坐标。

62)进一步对距离向时间τ_pi,1完成大气对流层延迟误差Δτ^Tro以及电离层的延迟Δτ^Ino误差改正，公式如下：

τ_pi,cor＝τ_pi,1-Δτ^Tro-Δτ^Ino，

经上述处理完后，灯柱p_i点得到误差补偿后的时间量纲坐标为(τ_pi,cor,SAR,t_pi,cor,SAR)。

所述立体平差处理生成灯柱控制点三维坐标包括以下步骤：

71)按照步骤11)-15)，逐个对N张SAR影像提取同名灯柱点对，针对某一灯柱点p_i得到同名点对集合

m表示同名灯柱点p_i所在SAR影像的标识，/>

和t_pi,cor,SARi分别表示在影像SARi上经过误差改正的距离向和方位向时间坐标。

72)对可用的立体SAR影像构建RD几何定位方程；

73)对某一灯柱点p_i的同名点对，构建立体观测方程如下：

式中,C_T＝[Lat_pi,Lon_pi,H_pi]^T为待求的三维坐标，其初始值

/>

74)对立体观测方程线性化，得到误差方程如下：

上式可简记为v＝AdC_T-l，。

75)求解立体观测误差方程，得到待求柱点p_i点的三维坐标改正值[dLat_pi,dLon_pi,dH_pi]^T。

76)将坐标改正值dC_T＝[dLat_pi,dLon_pi,dH_pi]^T补偿至初始值

中，更新/>

77)重复72)-76)，直到改正值小于阈值10^-5为止；

78)重复72)-77)，逐个计算提取出的灯柱点的三维坐标，直到所有灯柱点计算完为止。

有益效果

本发明的一种光学影像辅助的星载立体SAR影像控制点自动生成方法，与现有技术相比可实现不同升降轨SAR影像中灯柱PS点自动匹配，子像素级精确提取灯柱点在SAR影像上的坐标并完成异常值剔除，立体平差生成大量高精度可用的控制点，本方法可高效率、自动化的完成立体SAR影像上灯柱这一类PS点的匹配、提取，提取出的控制点可更好的服务基础测绘。

附图说明

图1为本发明方法顺序图；

图2为本发明所涉及的方法实施流程图；

图3为光学影像计算灯柱模板相关系数示意图；

图4为SAR影像灯柱PS点匹配示意图；

图5为SAR影像灯柱PS点坐标精确提取示意图。

具体实施方式

为使对本发明的结构特征及所达成的功效有更进一步的了解与认识，用以较佳的实施例及附图配合详细的说明，说明如下：

如图1和图2所示，本发明所述的一种光学影像辅助的星载立体SAR影像控制点自动生成方法，包括以下步骤：

第一步，基于高分光学影像上提取立体SAR影像上同名灯柱点的近似坐标：获取立体SAR影像覆盖区域所对应的高分光学遥感影像I₀，基于I₀检测到灯柱PS点后并得到其三维坐标，利用SAR影像的距离多普勒几何定位模型将灯柱的三维坐标反算至SAR影像像方，提取立体SAR影像上同名灯柱PS点的近似像点坐标(τ₀,t₀)。其具体步骤如下：

(1)获取对应区域的光学影像I₀后，为了使灯柱的阴影更加突出，对光学影像I₀进行高频提升滤波处理：

I＝I₀+aI_m，

式中，I₀为原始图像，a是锐化因子标量，I_m是使用非锐化蒙板后的图像，为I₀与其对应模糊图像之差，其中a值越高，表示锐化程度越高。

(2)如图3所示，在锐化后的光学图像I中提取任意包含灯柱阴影的模板图像T，模板大小设置为256×256，然后将模板T与图像I进行相关运算，得到图像I在像方(u,v)坐标处的相似度：

式中，I(x,y)和T(x,y)分别表示图像I和模板图像T在(x,y)坐标处的像素值，N₁,N₂表示模板图像T的大小，

和/>

分别表示图像I和模板图像T的平均像素值，ρ(u,v)为归一化互相关值因而不受图像亮度或对比度的变化影响，可提升模板匹配的效果。

(3)光学图像I中经模板匹配后每个灯柱对象存在多个像素，进一步对模板提取到的像素进行聚类处理，得到表示灯柱唯一的像素点p_i，处理如下：

式中，p_i表示计算m(p_i)漂移向量的空间坐标点，p_j代表p_i附近的点，g为具有带宽h的高斯核函数，而|| ||表示欧几里得距离算子。

(4)根据光学图像I上灯柱聚类点p_i的像素坐标(x_pi,I,y_pi,I)，获取其平面地理坐标(Lat_pi,I,Lon_pi,I)，并利用平面坐标在开源SRTM DEM数据上内插出对应高程H_pi,I，构成完整的三维坐标(Lat_pi,I,Lon_pi,I,H_pi,I)。

第二步，SAR影像灯柱PS点像方坐标粗略提取：对SAR影像进行阈值化处理，得到二值图，将光学影像中的路灯点坐标利用SAR的几何定位模型反算至其像方，应用迭代最近点法实现灯柱点集与SAR阈值处理后的亮点点集匹配，从而实现SAR影像灯柱PS点的匹配，并得到灯柱点对应SAR影像上的粗略坐标。其具体步骤如下：

(1)根据立体SAR影像的元数据信息构建距离多普勒RD几何定位方程：

式中，f_R、f_A分别表示距离方程和多普勒方程，τ、t分别表示距离向时间和方位向时间，C_S(t)和

分别表示成像t时刻SAR天线相位中心在WGS84坐标系下的位置矢量和速度矢量，C_T为观测目标在WGS84坐标系下的位置矢量，fdc(τ)表示在距离双程延时为τ处的多普中心频率；

(2)逐个将灯柱聚类点p_i的三维坐标(Lat_pi,I,Lon_pi,I,H_pi,I)利用RD定位几何方程反投影至SAR影像的像方，得到灯柱点集。

(3)根据一定阈值，本发明该阈值对灯柱点统计后自适应获取，SAR影像上亮点进行二值掩膜处理。

(4)如图4所示，利用迭代最近点法，实现反算在SAR像方的灯柱点集与SAR影像上亮点集匹配，SAR影像上只保留匹配到路灯点p_i的亮点，从而完成SAR影像灯柱点的匹配，并得到粗略匹配坐标(x_0,pi,SAR,y_0,pi,SAR)。

第三步，SAR影像灯柱PS点坐标精确提取：如图5所示，以近似坐标(x_0,pi,SAR,y_0,pi,SAR)为中心截取SAR影像块，进行点目标分析，精确提取灯柱p_i点的方位和距离坐标(x_pi,SAR,y_SAR)。其具体步骤如下：

(1)以近似坐标(x_0,pi,SAR,y_0,pi,SAR)为中心在SAR影像上截取32*32的影像块，并寻找最大峰值位置；

(2)以最大峰值处为中心截取32*32样本窗做二维离散傅里叶变换，32倍采样因子频域补零后进行二维逆离散傅里叶变换，得到升采样后的点目标影像块；

(3)对过采样后的影像块，以最大峰值处截取3×3影像子块，进行二维抛物面插值，椭圆抛物面插值模型如下：

f(x,y)＝a₅x²+a₄y²+a₃xy+a₂x+a₁y+a₀，

式中x,y为在3×3影像子块中的像素坐标，为a_i(1≤i≤6)为抛物面插值模型f(x,y)的参数，将3×3子块像素值代入上式中并通过最小二乘求解a_i。参数a_i求解后，进而获得灯柱p_i点精确坐标(x_pi,SAR,y_SAR)。

第四步，SAR影像灯柱PS点异常值剔除：理想灯柱PS点候选对象SCR应足够高，但由于存在误差，进一步通过计算灯柱PS点的相位噪声σ_φ，进行异常值剔除。其具体步骤如下：

(1)计算灯柱p_i点的峰值功率P_peak和杂波功率P_clutter。峰值功率P_peak是以峰值(x_pi,SAR,y_SAR)为中心为2×2分辨率像元的主瓣区域的积分信号能量，杂波功率_Pclutter是主瓣和旁瓣区域外的积分信号能量。

(2)计算灯柱p_i点的信杂比SCR并得到相位噪声

计算公式如下：

(3)如果灯柱p_i点计算出的σ_φ大于0.5rad，则予以剔除，否则保留。并将灯柱影像p_i在SAR影像上的(x_pi,SAR,y_SAR)坐标转为时间量纲的坐标(τ_pi,SAR,t_pi,SAR)，转化公式如下：

t＝t_start+y/PRF

τ＝τ_min+x/F_s

式中x、y为待转换的像素坐标，τ、t为转换后的时间量纲的坐标，分别表示距离向双程时延和方位向时间，PRF为方位向脉冲重复频率，F_s为距离向采样频率，τ_min为SAR影像近距端双程时延，t_start为方位向起始时间。

第五步，SAR影像灯柱PS点坐标误差补偿：对提取的灯柱点p_i影像坐标(τ_pi,SAR,t_pi,SAR)进行系统定位误差和大气延迟误差的校正，得到改正后的时间量纲坐标(τ_pi,cor,SAR,t_pi,cor,SAR)。其具体步骤如下：

(1)从SAR影像元数据中读取系统误差定标参数(Δτ^cal、Δt^cal)，补偿至灯柱p_i点影像坐标中，公式如下：

上式中τ_pi,1和t_pi,1分别为对灯柱p_i点影像初始时间量纲坐标(τ_pi,SAR,t_pi,SAR)进行系统误差补偿后的距离向和方位向时间坐标。

(2)进一步对距离向时间τ_pi,1完成大气对流层延迟误差Δτ^Tro以及电离层的延迟Δτ^Ino误差改正，公式如下：

τ_pi,cor＝τ_pi,1-Δτ^Tro-Δτ^Ino，

经上述处理完后，灯柱p_i点得到误差补偿后的时间量纲坐标为(τ_pi,cor,SAR,t_pi,cor,SAR)。

第六步，立体平差处理生成灯柱控制点三维坐标：针对待生成控制点区域内可用的N张(N＞2)立体SAR影像，按照上述第一步到第五步提取同名灯柱PS点对集合，然后进行立体平差处理，生成灯柱控制点的三维坐标。其具体步骤如下：

(1)按照第一步到第五步，逐个对N张SAR影像提取同名灯柱点对，针对某一灯柱点p_i得到同名点对集合

m表示同名灯柱点p_i所在SAR影像的标识，/>

和t_pi,cor,SARi分别表示在影像SARi上经过误差改正的距离向和方位向时间坐标。

(2)对可用的立体SAR影像构建RD几何定位方程；

(3)对某一灯柱点p_i的同名点对，构建立体观测方程如下：

式中C_T＝[Lat_pi,Lon_pi,H_pi]^T为待求的三维坐标，其初始值

/>

(4)对立体观测方程线性化，得到误差方程如下：

上式可简记为v＝AdC_T-l，。

(5)求解立体观测误差方程，得到待求柱点p_i点的三维坐标改正值[dLat_pi,dLon_pi,dH_pi]^T。

(6)将坐标改正值dC_T＝[dLat_pi,dLon_pi,dH_pi]^T补偿至初始值

中，更新/>

(7)重复(2)-(6)，直到改正值小于阈值10^-5为止；

(8)重复(2)-(7)，逐个计算提取出的灯柱点，直到所有灯柱点计算完为止。

下面以国产GF-3SAR卫星为例对本发明提出的方法进行说明：选取覆盖嵩山区域5景滑聚模式数据，数据中涵盖升降轨、左右视，不同组合立体角范围为70°到110°。光学影像采用谷歌地图嵩山区域最高分辨率图层的影像，其分辨率约为0.5m。为了使光学影像中灯柱更加突出及后续提高模板匹配精度，首先对光学影像进行滤波处理，选取一类灯柱作为图像模板进行模板匹配后，对匹配出的像素做聚类处理，得到45个灯柱点，通过SRTM DEM获取45个灯柱点的近似三维坐标。通过5景滑聚模式影像的几何定位模型，将45个灯柱点反算至SAR影像像方，并进一步提取精确的像方坐标，再经误差补偿后，进行立体平差处理，解算出灯柱点的三维坐标，经外业GPS测量检查点验证，平面精度优于1m，高程精度优于0.5m。

与现有技术相比，本发明具有如下优点与有益效果：

(1)基于光学影像辅助自动提取灯柱这一类PS点，而不是所有PS点，进而保证了后续控制点的可用性；

(2)可自动从不同升降轨的立体SAR影像中准确匹配灯柱PS点，实现异常值的剔除，从而保证控制点的精度。并且生成的控制点可为基础测绘、其他卫星几何校正提供较好的控制数据补充。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是本发明的原理，在不脱离本发明精神和范围的前提下本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明的范围内。本发明要求的保护范围由所附的权利要求书及其等同物界定。

Claims (7)

1.一种光学影像辅助的星载立体SAR影像控制点自动生成方法，其特征在于，包括以下步骤：

11)基于高分光学影像上提取立体SAR影像上同名灯柱点的近似坐标：获取立体SAR影像覆盖区域所对应的高分光学遥感影像I₀，基于I₀检测到灯柱点后并得到其近似三维坐标，利用SAR影像的距离多普勒几何定位模型将灯柱的近似三维坐标反算至SAR影像像方，完成立体SAR影像上同名灯柱点的近似像点坐标(τ₀,t₀)提取；

12)SAR影像灯柱PS点像方坐标粗略提取：对SAR影像进行阈值化处理，得到二值图，将光学影像中的路灯点坐标利用SAR的几何定位模型反算至其像方，应用迭代最近点法实现灯柱点集与SAR阈值处理后的亮点集进行匹配，从而实现SAR影像灯柱PS点的匹配，并得到灯柱点对应SAR影像上的粗略坐标；

13)SAR影像灯柱PS像方坐标精确提取：以近似坐标(x_0,pi,SAR,y_0,pi,SAR)为中心截取SAR影像块，进行点目标分析，精确提取灯柱p_i点的方位和距离坐标(x_pi,SAR,y_SAR)；

14)SAR影像灯柱PS点异常值剔除：潜在灯柱PS点候选对象SCR应足够高，但由于存在误差，进一步通过计算灯柱p_i点的相位噪声σ_φ，进行异常值剔除；

15)SAR影像灯柱PS点坐标误差补偿：对提取的灯柱点p_i影像坐标(τ_pi,SAR,t_pi,SAR)进行系统定位误差和大气延迟误差的校正，得到改正后的时间量纲坐标(τ_pi,cor,SAR,t_pi,cor,SAR)；

16)立体平差处理生成灯柱控制点三维坐标：针对待生成控制点区域内可用的N张(N＞2)立体SAR影像，按照上述步骤11)-15)提取同名灯柱点对集合，然后进行立体平差处理，生成灯柱控制点的三维坐标。

2.根据权利要求1所述的一种光学影像辅助的星载立体SAR影像控制点自动生成方法，其特征在于，所述基于高分光学影像上提取立体SAR影像上同名灯柱点的近似坐标包括以下步骤：

21)获取对应区域的光学影像I₀后，为使灯柱地面阴影更突出，对光学影像I₀进行高频提升滤波处理：

I＝I₀+aI_m，

式中，I₀为原始图像，a是锐化因子标量，I_m是使用非锐化蒙板后的图像，为I₀与其对应模糊图像之差，其中a值越高，表示锐化程度越高。

22)在锐化后的光学图像I中提取任意包含灯柱阴影的模板图像T，模板大小设置为256×256，然后将模板与图像I进行相关运算，得到图像I在像方坐标(u,v)处的相似度：

式中，(u,v)表示图像I上的像方坐标，I(x,y)和T(x,y)分别表示图像I和模板图像T在(x,y)坐标处的像素值，

和/>

分别表示图像I和模板图像T的像素值的平均值，ρ(u,v)为归一化互相关系数，ρ(u,v)不受图像亮度或对比度的变化影响，因此可提升模板匹配的效果。

23)光学图像I中经模板匹配后得到每个灯柱对象均包含多个像素，进一步对模板提取到的像素进行聚类处理，得到表示灯柱唯一的像素点p_i，具体处理如下：

式中，p_i表示计算m(p_i)漂移向量的空间坐标点，p_j代表p_i附近的点，g为具有带宽h的高斯核函数，而|| ||表示欧几里得距离算子。

24)根据光学图像I上灯柱聚类点p_i的像素坐标(x_pi,I,y_pi,I)，获取其平面地理坐标(Lat_pi,I,Lon_pi,I)，并利用平面坐标在开源SRTM DEM数据上内插出对应高程H_pi,I，构成完整的三维坐标(Lat_pi,I,Lon_pi,I,H_pi,I)。

3.根据权利要求1所述的一种光学影像辅助的星载立体SAR影像控制点自动生成方法，其特征在于，所述SAR影像灯柱PS点像方坐标粗略提取包括以下步骤：

31)根据立体SAR影像的元数据信息构建距离多普勒RD几何定位方程：

式中，f_R、f_A分别表示距离方程和多普勒方程，τ、t分别表示距离向时间和方位向时间，C_S(t)和

分别表示成像t时刻SAR天线相位中心在WGS84坐标系下的位置矢量和速度矢量，C_T为观测目标在WGS84坐标系下的位置矢量，fdc(τ)表示在距离双程延时为τ处的多普中心频率；

32)逐个将灯柱聚类点p_i的三维坐标(Lat_pi,I,Lon_pi,I,H_pi,I)利用RD定位几何方程反投影至SAR影像的像方，得到灯柱点集。

33)根据一定阈值，该阈值由对灯柱点像素值统计后获取，对SAR影像上亮点进行二值掩膜处理。

34)利用迭代最近点法，实现反算在SAR像方的灯柱点集与SAR影像上亮点集匹配，SAR影像上只保留匹配到路灯点p_i的亮点，从而完成SAR影像灯柱点的匹配，并得到粗略匹配坐标(x_0,pi,SAR,y_0,pi,SAR)。

4.根据权利要求1所述的一种光学影像辅助的星载立体SAR影像控制点自动生成方法，其特征在于，所述SAR影像间灯柱PS像方坐标精确提取包括以下步骤：

41)以近似坐标(x_0,pi,SAR,y_0,pi,SAR)为中心在SAR影像上截取32*32的影像块，并寻找最大峰值位置；

42)以最大峰值处为中心截取32*32样本窗做二维离散傅里叶变换，32倍采样因子频域补零后进行二维逆离散傅里叶变换，得到升采样后的点目标影像块；

43)对过采样后的影像块，以最大峰值处截取3×3影像子块，进行二维抛物面插值，椭圆抛物面插值模型如下：

f(x,y)＝a₅x²+a₄y²+a₃xy+a₂x+a₁y+a₀，

式中，x,y为在3×3影像子块中的像素坐标，为a_i(1≤i≤6)为抛物面插值模型f(x,y)的参数，将3×3子块像素值代入上式中并通过最小二乘求解a_i。参数a_i求解后，进而获得灯柱p_i点精确坐标(x_pi,SAR,y_SAR)。

5.根据权利要求1所述的一种光学影像辅助的星载立体SAR影像控制点自动生成方法，其特征在于，所述SAR影像灯柱PS点异常值剔除包括以下步骤：

51)计算灯柱p_i点的峰值功率P_peak和杂波功率P_clutter。峰值功率P_peak是以峰值(x_pi,SAR,y_SAR)为中心为2×2分辨率像元的主瓣区域的积分信号能量，杂波功率P_clutter是主瓣和旁瓣区域外的积分信号能量。

52)计算灯柱p_i点的信杂比SCR并得到相位噪声

计算公式如下：

53)如果灯柱p_i点计算出的σ_φ大于0.5rad，则予以剔除，否则保留。并将灯柱影像p_i在SAR影像上的(x_pi,SAR,y_SAR)坐标转为时间量纲的坐标(τ_pi,SAR,t_pi,SAR)。转化公式如下：

式中，x、y为待转换的像素坐标，τ、t为转换后的时间量纲的坐标，分别表示距离向双程时延和方位向时间，PRF为方位向脉冲重复频率，F_s为距离向采样频率，τ_min为SAR影像近距端双程时延，t_start为方位向起始时间。

6.根据权利要求1所述的一种光学影像辅助的星载立体SAR影像控制点自动生成方法，其特征在于，所述SAR影像灯柱PS点坐标误差补偿包括以下步骤：

61)从SAR影像元数据中读取系统误差定标参数(Δτ^cal、Δt^cal)，补偿至灯柱p_i点影像坐标中，公式如下：

式中，τ_pi,1和t_pi,1分别为对灯柱p_i点影像初始时间量纲坐标(τ_pi,SAR,t_pi,SAR)进行系统误差补偿后的距离向和方位向时间坐标。

62)进一步对距离向时间τ_pi,1完成大气对流层延迟误差Δτ^Tro以及电离层的延迟Δτ^Ino误差改正，公式如下：

τ_pi,cor＝τ_pi,1-Δτ^Tro-Δτ^Ino，

经上述处理完后，灯柱p_i点得到误差补偿后的时间量纲坐标为(τ_pi,cor,SAR,t_pi,cor,SAR)。

7.根据权利要求1所述的一种光学影像辅助的星载立体SAR影像控制点自动生成方法，其特征在于，所述立体平差处理生成灯柱控制点三维坐标包括以下步骤：

71)按照步骤11)-15)，逐个对N张SAR影像提取同名灯柱点对，针对某一灯柱点p_i得到同名点对集合

m表示同名灯柱点p_i所在SAR影像的标识，/>

和t_pi,cor,SARi分别表示在影像SARi上经过误差改正的距离向和方位向时间坐标。

72)对可用的立体SAR影像构建距离多普勒几何定位方程；

73)对某一灯柱点p_i的同名点对，构建立体观测方程如下：

式中，C_T＝[Lat_pi,Lon_pi,H_pi]^T为待求的三维坐标，其初始值

74)对立体观测方程线性化，得到误差方程如下：

上式可简记为v＝AdC_T-l，。

75)求解立体观测误差方程，得到待求柱点p_i点的三维坐标改正值[dLat_pi,dLon_pi,dH_pi]^T。

76)将坐标改正值dC_T＝[dLat_pi,dLon_pi,dH_pi]^T补偿至初始值

中，更新/>

77)重复72)-76)，直到改正值小于阈值10^-5为止；

78)重复72)-77)，逐个计算提取出的灯柱点，直到所有灯柱点坐标计算完为止。

Images