一种基于偏振三维成像的遥感地形三维测绘方法

1.本发明属于遥感成像技术领域，具体涉及一种基于偏振三维成像的遥感地形三维测绘方法。


背景技术：

2.航天光学遥感技术通常指利用卫星、航天飞机等在高空对地面上的目标进行远距离的探测，通常卫星平台搭载着空间相机或成像光谱仪等设备来收集目标所反射或发射出的辐射信息，利用所捕获到的信息进行光电转换、信息处理等来实现目标的监测。同时，利用卫星遥感影像的优势，在国土资源统计与监测、灾情监测与城市变迁等领域有着重要应用。随着科技水平的飞速发展，原始的二维遥感影像已无法满足应用需求，因此，需要寻求一种三维成像技术来对地形建立三维重建模型，而三维重建的主要目标是通过计算获取目标区域地形的数字高程信息。
3.目前，常见的三维重建方法是利用航空航天影像结合摄影测量理论，通过获取立体像对进行三维重建，其中，遥感卫星的立体像对是指对相同目标区域从不同的角度进行拍摄，获取到两幅或者多幅图像，再利用模拟人眼的立体视觉，最终实现目标地形区域的三维重建。通过确定立体像对的内、外方位元素，结合共线条件以及两幅或多幅匹配图像(通常是同一卫星对同一目标区域不同角度所拍摄得到的多幅，具有一定重叠区域的图像)上的像点坐标进行影像匹配，得到立体像对上的同名像点，再进一步构建出两条同名射线，通过两条同名射线的交点即可计算出对应地面点的三维坐标，最终实现基于立体像对的三维重建。
4.由于卫星遥感图像地形多为山地且起伏较大，影像匹配得到的同名点数量往往会存在极大的限制，导致最终重建效果受到影响；另外，虽然影像匹配方法较多，但其往往存在对噪声与光照影响精度与算法运行速度上无法同时兼顾的缺陷；再者，由于基于立体像对遥感影像三维重建技术匹配时往往需要多幅图片中包含相同的重叠区域，因此对于大范围大场景的目标区域拼接往往较为繁琐复杂，无法高效精确实现。


技术实现要素：

5.为了解决现有技术中存在的上述问题，本发明提供了一种基于偏振三维成像的遥感地形三维测绘方法。本发明要解决的技术问题通过以下技术方案实现：
6.本发明提供一种基于偏振三维成像的遥感地形三维测绘方法，包括：
7.采集预设角度下不同目标地形区域的偏振图像；
8.针对每个目标地形区域，计算目标地形区域表面的偏振度，并根据所述偏振度计算所述目标地形区域表面的入射光的入射角θ和方位角
9.根据所述入射角θ和所述方位角构建所述目标地形区域表面的法向量与所述入射角θ及所述方位角之间的第一映射关系，获得第一法向量信息；
10.计算太阳的天顶角θz与方位角θs，构建光源矢量与所述天顶角θz及所述方位角θs之间的第二映射关系，获得所述目标地形区域表面的光源矢量；
11.构建所述目标地形区域表面的辐射模型，并结合所述目标地形区域表面的光源矢量计算先验法向量信息；
12.基于所述先验法向量信息对所述第一法向量信息进行修正，得到所述目标地形区域表面的修正法向量信息；
13.根据所述修正法向量信息和所述预设角度下不同目标地形区域的偏振图像，进行目标区域的三维场景重建。
14.在本发明的一个实施例中，所述预设角度包括0
°
、45
°
、90
°
和135
°
；
15.针对每个目标地形区域，计算目标地形区域表面的偏振度，并根据所述偏振度计算所述目标地形区域表面的入射光的入射角θ和方位角的步骤，包括：
16.针对每个目标地形区域，获取采集预设角度下的偏振图像时通过0
°
偏光片的光强i
00
、通过45
°
偏光片的光强i
+450
、通过90
°
偏光片的光强i
900
和通过135
°
方向线偏振器的光强i-45
°
，并计算斯托克斯参数s0、s1和s2：
17.s0＝i0°
+i
90
°
；
18.s1＝i0°-i
90
°
；
19.s2＝i
+45
°-i-45
°
；
20.根据所述斯托克斯参数s0、s1和s2，计算目标地形区域表面的偏振度：
[0021][0022]
根据所述斯托克斯参数s0、s1及所述偏振度，计算目标地形区域表面的入射光的入射角θ和方位角
[0023]
在本发明的一个实施例中，根据所述斯托克斯参数s0、s1及所述偏振度，计算目标地形区域表面的入射光的入射角θ和方位角的步骤，包括：
[0024]
根据所述斯托克斯参数s0、s1计算目标地形区域表面的入射光的方位角
[0025][0026]
根据所述偏振度dolp计算入射光的入射角θ：
[0027][0028]
式中，n表示目标地形区域表面的折射率，n＝1.5。
[0029]
在本发明的一个实施例中，根据所述入射角θ和所述方位角构建所述目标地形区域表面的法向量与所述入射角θ及所述方位角之间的第一映射关系，获得第一法向量
信息的步骤，包括：
[0030]
根据所述入射角θ和所述方位角构建所述目标地形区域表面的法向量与所述入射角θ及所述方位角之间的第一映射关系：
[0031][0032]
式中，n
x
、ny、nz分别表示所述目标地形区域表面的法向量在x方向、y方向和z方向上的分量，其中，z方向与所述目标地形区域表面垂直，x方向与y方向垂直且二者所在平面与所述z方向垂直；
[0033]
根据所述第一映射关系计算第一法向量信息d。
[0034]
在本发明的一个实施例中，计算太阳的天顶角θz与方位角θs，构建光源矢量与所述天顶角θz及所述方位角θs之间的第二映射关系，获得所述目标地形区域表面的光源矢量的步骤，包括：
[0035]
根据所述目标地形区域的被拍摄时间，计算太阳赤纬角：
[0036][0037]
式中，n表示积日数；
[0038]
根据所述太阳赤纬角δ、目标地形区域的纬度φ和太阳时角ω计算太阳的高度角：
[0039]
sinhs＝sinδ
·
sinφ+cosδ
·
cosφ
·
cosω；
[0040]
式中，ω＝(12-t)
×
15
°
，t表示目标地形区域被拍摄时的当地时间；
[0041]
根据太阳赤纬角δ、目标地形区域的经纬度信息和太阳高度角hs计算太阳的方位角：
[0042][0043]
式中，φ表示目标地形区域的纬度；
[0044]
根据太阳的高度角hs计算太阳的天顶角θz：
[0045]
θz＝90
°‑hs
；
[0046]
构建光源矢量与所述天顶角θz及所述方位角θs之间的第二映射关系：
[0047][0048]
根据所述第二映射关系，获得所述目标地形区域表面的光源矢量(ps,qs,-1)。
[0049]
在本发明的一个实施例中，构建所述目标地形区域表面的辐射模型，并结合所述目标地形区域表面的光源矢量计算先验法向量信息的步骤，包括：
[0050]
根据朗伯反射模型、所述目标地形区域表面的光源矢量(ps,qs,-1)和光滑性约束，构建所述目标地形区域表面的辐射模型：
[0051][0052]
所述光滑性约束为：
[0053][0054]
式中，i(x',y')表示目标地形区域(x',y')处在0
°
和90
°
下的偏振图像的光源强度之和，ρ表示预设的偏振图像反照率，p
x
,py,q
x
,qy分别为p(x',y')、q(x',y')对应于x方向和y方向上的导数；
[0055]
根据所述辐射模型及光滑性约束，计算得到所述目标地形区域表面(x',y')处的先验法向量信息y＝(p(x',y'),q(x',y'),-1)。
[0056]
在本发明的一个实施例中，基于所述先验法向量信息对所述第一法向量信息进行修正，得到所述目标地形区域表面修正后的法向量的步骤，包括：
[0057]
按照如下公式计算二元算子
[0058][0059]
式中，t表示非线性算子；
[0060]
基于所述二元算子对第一法向量信息d进行修正：
[0061][0062]
式中，n
cor
表示所述目标地形区域表面的修正法向量信息。
[0063]
在本发明的一个实施例中，所述不同目标区域包括第一目标区域和第二目标区域，所述第一目标区域与所述第二目标区域至少存在部分重叠区域；
[0064]
根据所述修正法向量信息和所述预设角度下不同目标地形区域的偏振图像，进行目标区域的三维场景重建的步骤，包括：
[0065]
从所述第一目标区域的任意一幅偏振图像中选取一特征点，所述特征点位于重叠区域内；
[0066]
获取所述特征点在所述第二目标区域中的位置信息，并建立所述第一目标区域的偏振图像与所述第二目标区域的偏振图像之间的变换关系；
[0067]
根据所述变换关系拼接所述第一目标区域的偏振图像与所述第二目标区域的偏振图像；
[0068]
根据拼接后的图像、第一目标区域表面的修正法向量信息和第二目标区域表面的修正法向量信息进行三维重建。
[0069]
在本发明的一个实施例中，根据拼接后的图像、第一目标区域表面的修正法向量信息和第二目标区域表面的修正法向量信息，按照如下公式进行三维重建：
[0070][0071]
式中，f{}和f-1
{}分别表示离散傅里叶变换和傅里叶逆变换，m、n分别表示拼接后的图像中像素点的行数及列数，u、ν分别表示离散微算子在x和y方向上的傅里叶系数。
[0072]
与现有技术相比，本发明的有益效果在于：
[0073]
本发明提供一种基于偏振三维成像的遥感地形三维测绘方法，通过利用偏振图像中的信息进行三维重建，避免了现有技术中基于立体像对遥感影像三维重建时无法避开的图像匹配的过程，提高了运算速度；与此同时，本发明利用偏振三维成像进行逐像素拼接，可以快速高效地实现大场景下的卫星遥感地形三维重建，解决了现有基于立体像对遥感影像三维重建过程中需要多幅图像且只能对重叠区域进行三维重建的问题。另外，基于偏振三维重建的遥感影像地形恢复技术具有被动式、高精度、远距离、干扰较小的优势，使得利用卫星遥感图像进行三维重建的精度进一步提升成为可能。
[0074]
以下将结合附图及实施例对本发明做进一步详细说明。
附图说明
[0075]
图1是本发明实施例提供的基于偏振三维成像的遥感地形三维测绘方法的一种流程图；
[0076]
图2是本发明实施例提供的基于偏振三维成像的遥感地形三维测绘方法的另一种流程图；
[0077]
图3是本发明实施例提供的目标地形区域表面法向量与入射光的方位角和入射角的示意图；
[0078]
图4是本发明实施例提供的目标地形区域表面法向量与太阳天顶角和方位角的示意图；
[0079]
图5是本发明实施例提供的逐像素拼接的示意图。
具体实施方式
[0080]
下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。
[0081]
图1是本发明实施例提供的基于偏振三维成像的遥感地形三维测绘方法的一种流程图，图2是本发明实施例提供的基于偏振三维成像的遥感地形三维测绘方法的另一种流程图。如图1-2所示，本发明提供一种基于偏振三维成像的遥感地形三维测绘方法，包括：
[0082]
s1、采集预设角度下不同目标地形区域的偏振图像；
[0083]
s2、针对每个目标地形区域，计算目标地形区域表面的偏振度，并根据偏振度计算目标地形区域表面的入射光的入射角θ和方位角
[0084]
s3、根据入射角θ和方位角构建目标地形区域表面的法向量与入射角θ及方位角之间的第一映射关系，获得第一法向量信息；
[0085]
s4、计算太阳的天顶角θz与方位角θs，构建光源矢量与天顶角θz及方位角θs之间的第二映射关系，获得目标地形区域表面的光源矢量；
[0086]
s5、构建目标地形区域表面的辐射模型，并结合目标地形区域表面的光源矢量计算先验法向量信息；
[0087]
s6、基于先验法向量信息对第一法向量信息进行修正，得到目标地形区域表面的修正法向量信息；
[0088]
s7、根据修正法向量信息和预设角度下不同目标地形区域的偏振图像，进行目标区域的三维场景重建。
[0089]
本实施例中，预设角度包括0
°
、45
°
、90
°
和135
°
。具体地，可以利用卫星上所搭载的成像探测器cmos相机采集目标地形区域的反射光，得到每个目标地形区域在0
°
、45
°
、90
°
和135
°
的偏振图像。
[0090]
上述步骤s2中，针对每个目标地形区域，计算目标地形区域表面的偏振度，并根据偏振度计算目标地形区域表面的入射光的入射角θ和方位角的步骤，包括：
[0091]
s201、针对每个目标地形区域，获取采集预设角度下的偏振图像时通过0
°
偏光片的光强i0°
、通过45
°
偏光片的光强i
+45
°
、通过90
°
偏光片的光强i
90
°
和通过135
°
方向线偏振器的光强i-45
°
，并计算斯托克斯参数s0、s1和s2：
[0092]
s0＝i0°
+i
90
°
；
[0093]
s1＝i0°-i
90
°
；
[0094]
s2＝i
+45
°-i-45
°
；
[0095]
s202、根据斯托克斯参数s0、s1和s2，计算目标地形区域表面的偏振度：
[0096][0097]
s203、根据斯托克斯参数s0、s1及所述偏振度，计算目标地形区域表面的入射光的入射角θ和方位角
[0098]
应当理解，利用两个或两个以上具有不同偏振态的光子可以描述光的偏振状态，偏振状态及其和由斯托克参数描述，其分量为四个斯托克斯参数：
[0099][0100]
其中：
[0101][0102]
在实际应用中，斯托克斯参数可以通过偏振器测量获得：
[0103][0104]
上式中，i0°
表示通过0
°
偏光片的光强，i
90
°
是通过90
°
偏光片的光强，i
+45
°
是与水平正方向成45
°
偏光片的光强，i-45
°
是135
°
方向的线偏振器得到的光强，i
rh
和i
lh
分别表示右旋圆偏振光分量和左旋圆偏振光分量。
[0105]
进一步地，利用斯托克斯参数可以直接确定任何光信号的偏振度(dop)，偏振度表征偏振分量的比例：
[0106][0107]
由于本实施例中采集偏振图像时的光束是线性偏振，因而右旋圆偏振光分量与左旋圆偏振光分量均为零，则目标地形区域表面的偏振度可表示为：
[0108][0109]
式中，dolp一般称为线性偏振度，表征了光束的线性偏振度。
[0110]
上述步骤s203中，根据斯托克斯参数s0、s1及偏振度，计算目标地形区域表面的入射光的入射角θ和方位角的步骤，包括：
[0111]
根据斯托克斯参数s0、s1计算目标地形区域表面的入射光的方位角
[0112][0113]
根据偏振度dolp计算入射光的入射角θ：
[0114][0115]
式中，n表示目标地形区域表面的折射率，n＝1.5。
[0116]
图3是本发明实施例提供的目标地形区域表面法向量与入射光的方位角和入射角的示意图。请参见图3，上述步骤s3中，根据入射角θ和所述方位角构建目标地形区域表面的法向量与入射角θ及方位角之间的第一映射关系，获得第一法向量信息的步骤，包括：
[0117]
根据入射角θ和方位角构建目标地形区域表面的法向量与入射角θ及方位角之间的第一映射关系：
[0118][0119]
式中，n
x
、ny、nz分别表示目标地形区域表面的法向量在x方向、y方向和z方向上的分量，其中，z方向与目标地形区域表面垂直，x方向与y方向垂直且二者所在平面与z方向垂直；
[0120]
根据第一映射关系计算第一法向量信息d。
[0121]
需要说明的是，在光强一定时存在对应多个偏振器旋转角度的情况，因此通过光强来计算入射光的方位角时会存在一个180
°
的多值性问题，也就无法准确获取目标地形区域表面的法向量信息。
[0122]
基于上述分析，本实施例还需要进一步对根据第一映射关系计算获得的第一法向量信息d进行修正。
[0123]
图4是本发明实施例提供的目标地形区域表面法向量与太阳天顶角和方位角的示意图。如图4所示，上述步骤s4中，计算太阳的天顶角θz与方位角θs，构建光源矢量与天顶角θz及方位角θs之间的第二映射关系，获得目标地形区域表面的光源矢量的步骤，包括：
[0124]
s401、根据目标地形区域的被拍摄时间，计算太阳赤纬角：
[0125][0126]
式中，n表示积日数；
[0127]
s402、根据太阳赤纬角δ、目标地形区域的纬度φ和太阳时角ω计算太阳的高度角：
[0128]
sinhs＝sinδ
·
sinφ+cosδ
·
cosφ
·
cosω；
[0129]
式中，ω＝(12-t)
×
15
°
，t表示目标地形区域被拍摄时的当地时间；
[0130]
s403、根据太阳赤纬角δ、目标地形区域的经纬度信息和太阳高度角hs计算太阳的方位角：
[0131][0132]
式中，φ表示目标地形区域的纬度；
[0133]
s404、根据太阳的高度角hs计算太阳的天顶角θz：
[0134]
θz＝90
°‑hs
；
[0135]
s405、构建光源矢量与天顶角θz及方位角θs之间的第二映射关系：
[0136][0137]
s406、根据第二映射关系，获得目标地形区域表面的光源矢量(ps,qs,-1)。
[0138]
具体而言，首先根据目标地形区域的经纬度信息与被拍摄时间，计算太阳的方位
角θs与高度角hs为：
[0139][0140]
sinhs＝sinδ
·
sinφ+cosδ
·
cosφ
·
cosω；
[0141]
式中，δ为太阳中心点和地球中心点的连线与赤道所在平面之间所夹的锐角，又称为太阳赤纬角，φ表示目标地形区域的纬度，ω表示太阳时角。
[0142]
其中，赤纬角δ可由下式解得：
[0143][0144]
n表示从拍摄时间所在年份的1月1日到拍摄当天的日期之间的天数，又称积日数。
[0145]
太阳时角ω的计算公式为：
[0146]
ω＝(12-t)
×
15
°
[0147]
t表示目标地形区域被拍摄时的当地时间。
[0148]
请参见图4，太阳的高度角hs与天顶角θz互为余角，因此太阳的天顶角可由高度角hs确定：
[0149]
θz＝90
°‑hs
。
[0150]
上述步骤s5中，构建目标地形区域表面的辐射模型，并结合目标地形区域表面的光源矢量计算先验法向量信息的步骤，包括：
[0151]
s501、根据朗伯反射模型、目标地形区域表面的光源矢量(ps,qs,-1)和光滑性约束，构建目标地形区域表面的辐射模型：
[0152][0153]
光滑性约束为：
[0154][0155]
式中，光滑性约束表示相邻像素点法向量方向接近，i(x',y')表示目标地形区域(x',y')处在0
°
和90
°
下的偏振图像的光源强度之和，ρ表示预设的偏振图像反照率，p
x
,py,q
x
,qy分别为p(x',y')、q(x',y')对应于x方向和y方向上的导数；
[0156]
s502、根据辐射模型及光滑性约束，计算得到目标地形区域表面(x',y')处的先验法向量信息y＝(p(x',y'),q(x',y'),-1)。
[0157]
上述步骤s6中，基于先验法向量信息对第一法向量信息进行修正，得到目标地形区域表面的修正法向量信息的步骤，包括：
[0158]
按照如下公式计算二元算子
[0159][0160]
式中，t表示非线性算子；
[0161]
基于所述二元算子对第一法向量信息d进行修正：
[0162][0163]
式中，n
cor
表示目标地形区域表面的修正法向量信息。
[0164]
本实施例利用法向梯度场的先验信息所提供的目标地形区域表面在校正过程中的变化趋势，解决了方位角模糊所导致表面法向量不准确的问题，得到修正后的目标表面法线，进一步获得目标表面梯度场，然后可利用表面梯度场重建方法可以恢复目标表面形状。
[0165]
可选地，不同目标区域包括第一目标区域和第二目标区域，第一目标区域与第二目标区域至少存在部分重叠区域；
[0166]
上述步骤s7中，根据修正法向量信息和预设角度下不同目标地形区域的偏振图像，进行目标区域的三维场景重建的步骤，包括：
[0167]
s701、从第一目标区域的任意一幅偏振图像中选取一特征点，特征点位于重叠区域内；
[0168]
s702、获取特征点在第二目标区域中的位置信息，并建立第一目标区域的偏振图像与第二目标区域的偏振图像之间的变换关系；
[0169]
s703、根据变换关系拼接第一目标区域的偏振图像与第二目标区域的偏振图像；
[0170]
s704、根据拼接后的图像、第一目标区域表面的修正法向量信息和第二目标区域表面的修正法向量信息进行三维重建。
[0171]
对于实际中大场景大区域的目标区域地形三维重建时，由于相机视场有限，无法使一幅图像中包含所有目标区域，因此需要对多区域的图像进行拼接。步骤s701～s702中，首先从第一目标区域的任意一幅偏振图像中选取中一个特征点，然后结合这一特征点在第二目标区域的任意一幅偏振图像中的位置信息，建立两幅偏振图像之间的变换关系，从而实现图像的无缝拼接。同时，由于偏振三维成像是基于每一个独立的像素点进行的，因此可以同时实现大场景的逐像素拼接与重建。
[0172]
图5是本发明实施例提供的逐像素拼接的示意图。如图5所示，以方格代表偏振图像中的像素、实心方格表示选取的特征点，其中，i(m,n)表示在第一目标区域的偏振图像中位于(m,n)处的像素点，i(i,j)表示第二目标区域的偏振图像中位于(i,j)处的像素点，结合两幅图像之间的变换关系以及修正法向量信息即可实现大场景下遥感地形的逐像素重建。
[0173]
具体而言，将目标曲面记作z(x,y)，利用可积函数的梯度表达形式来表示目标区域的表面法线：
[0174][0175]
结合法向量与入射角θ及方位角之间的映射关系，联立可得：
[0176]
[0177][0178]
定义函数u(z)表示利用梯度场求解目标曲面形状与实际形状之间的误差：
[0179]
u(z)＝(z
x-p)2+(z
y-q)2[0180]
其中，z
x
与zy分别表示目标曲面z(x,y)在x方向和y方向上的偏导数。
[0181]
为了得到最真实的目标曲面z(x,y)，需要使得目标曲面上所有点的误差最小化，因此可定义目标代价函数f(z)来进行优化处理：
[0182]
f(z)＝∫∫u(z)dxdy。
[0183]
对于全局连续的法线梯度场，通常采用全局或局部积分的方法。局部积分方法容易被待处理数据中的噪声所干扰，而实际采集初始数据时往往会因为环境杂散光以及本身探测器的影响，存在不可避免的噪声影响，而全局积分方法对噪声具有更好的鲁棒性，重建的表面更光滑。因此，本实施例采用全局积分方法来求解目标区域曲面使得误差最小化，并利用frankot chellappa算法获得最终重建结果。
[0184]
frankot chellappa算法主要是通过在频率域空间中来重建出目标区域的表面函数，再通过傅里叶逆变换来实现目标的三维重建模型。可选地，按照如下公式进行三维重建：
[0185][0186]
式中，f{}和f-1
{}分别表示离散傅里叶变换和傅里叶逆变换，m、n分别表示拼接后的图像中像素点的行数及列数，u、ν分别表示离散微算子在x和y方向上的傅里叶系数。
[0187]
通过上述各实施例可知，本发明的有益效果在于：
[0188]
本发明提供一种基于偏振三维成像的遥感地形三维测绘方法，通过利用偏振图像中的信息进行三维重建，避免了现有技术中基于立体像对遥感影像三维重建时无法避开的图像匹配的过程，提高了运算速度；与此同时，本发明利用偏振三维成像进行逐像素拼接，可以快速高效地实现大场景下的卫星遥感地形三维重建，解决了现有基于立体像对遥感影像三维重建过程中需要多幅图像且只能对重叠区域进行三维重建的问题。另外，基于偏振三维重建的遥感影像地形恢复技术具有被动式、高精度、远距离、干扰较小的优势，使得利用卫星遥感图像进行三维重建的精度进一步提升成为可能。
[0189]
在本发明的描述中，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。
[0190]
参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例
或示例中以合适的方式结合。此外，本领域的技术人员可以将本说明书中描述的不同实施例或示例进行接合和组合。
[0191]
尽管在此结合各实施例对本技术进行了描述，然而，在实施所要求保护的本技术过程中，本领域技术人员通过查看所述附图、公开内容、以及所附权利要求书，可理解并实现所述公开实施例的其他变化。
[0192]
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

技术特征：
1.一种基于偏振三维成像的遥感地形三维测绘方法，其特征在于，包括：采集预设角度下不同目标地形区域的偏振图像；针对每个目标地形区域，计算目标地形区域表面的偏振度，并根据所述偏振度计算所述目标地形区域表面的入射光的入射角θ和方位角根据所述入射角θ和所述方位角构建所述目标地形区域表面的法向量与所述入射角θ及所述方位角之间的第一映射关系，获得第一法向量信息；计算太阳的天顶角θ
z
与方位角θ
s
，构建光源矢量与所述天顶角θ
z
及所述方位角θ
s
之间的第二映射关系，获得所述目标地形区域表面的光源矢量；构建所述目标地形区域表面的辐射模型，并结合所述目标地形区域表面的光源矢量计算先验法向量信息；基于所述先验法向量信息对所述第一法向量信息进行修正，得到所述目标地形区域表面的修正法向量信息；根据所述修正法向量信息和所述预设角度下不同目标地形区域的偏振图像，进行目标区域的三维场景重建。2.根据权利要求1所述的基于偏振三维成像的遥感地形三维测绘方法，其特征在于，所述预设角度包括0
°
、45
°
、90
°
和135
°
；针对每个目标地形区域，计算目标地形区域表面的偏振度，并根据所述偏振度计算所述目标地形区域表面的入射光的入射角θ和方位角的步骤，包括：针对每个目标地形区域，获取采集预设角度下的偏振图像时通过0
°
偏光片的光强i0°
、通过45
°
偏光片的光强i
+45
°
、通过90
°
偏光片的光强i
90
°
和通过135
°
方向线偏振器的光强i-45
°
，并计算斯托克斯参数s0、s1和s2：s0＝i0°
+i
90
°
；s1＝i0°-i
90
°
；s2＝i
+45
°-i-45
°
；根据所述斯托克斯参数s0、s1和s2，计算目标地形区域表面的偏振度：根据所述斯托克斯参数s0、s1及所述偏振度，计算目标地形区域表面的入射光的入射角θ和方位角3.根据权利要求2所述的基于偏振三维成像的遥感地形三维测绘方法，其特征在于，根据所述斯托克斯参数s0、s1及所述偏振度，计算目标地形区域表面的入射光的入射角θ和方位角的步骤，包括：根据所述斯托克斯参数s0、s1计算目标地形区域表面的入射光的方位角
根据所述偏振度dolp计算入射光的入射角θ：式中，n表示目标地形区域表面的折射率，n＝1.5。4.根据权利要求2所述的基于偏振三维成像的遥感地形三维测绘方法，其特征在于，根据所述入射角θ和所述方位角构建所述目标地形区域表面的法向量与所述入射角θ及所述方位角之间的第一映射关系，获得第一法向量信息的步骤，包括：根据所述入射角θ和所述方位角构建所述目标地形区域表面的法向量与所述入射角θ及所述方位角之间的第一映射关系：式中，n
x
、n
y
、n
z
分别表示所述目标地形区域表面的法向量在x方向、y方向和z方向上的分量，其中，z方向与所述目标地形区域表面垂直，x方向与y方向垂直且二者所在平面与所述z方向垂直；根据所述第一映射关系计算第一法向量信息d。5.根据权利要求4所述的基于偏振三维成像的遥感地形三维测绘方法，其特征在于，计算太阳的天顶角θ
z
与方位角θ
s
，构建光源矢量与所述天顶角θ
z
及所述方位角θ
s
之间的第二映射关系，获得所述目标地形区域表面的光源矢量的步骤，包括：根据所述目标地形区域的被拍摄时间，计算太阳赤纬角：式中，n表示积日数；根据所述太阳赤纬角δ、目标地形区域的纬度φ和太阳时角ω计算太阳的高度角：sinh
s
＝sinδ
·
sinφ+cosδ
·
cosφ
·
cosω；式中，ω＝(12-t)
×
15
°
，t表示目标地形区域被拍摄时的当地时间；根据太阳赤纬角δ、目标地形区域的经纬度信息和太阳高度角h
s
计算太阳的方位角：
式中，φ表示目标地形区域的纬度；根据太阳的高度角h
s
计算太阳的天顶角θ
z
：θ
z
＝90
°‑
h
s
；构建光源矢量与所述天顶角θ
z
及所述方位角θ
s
之间的第二映射关系：根据所述第二映射关系，获得所述目标地形区域表面的光源矢量(p
s
,q
s
,-1)。6.根据权利要求5所述的基于偏振三维成像的遥感地形三维测绘方法，其特征在于，构建所述目标地形区域表面的辐射模型，并结合所述目标地形区域表面的光源矢量计算先验法向量信息的步骤，包括：根据朗伯反射模型、所述目标地形区域表面的光源矢量(p
s
,q
s
,-1)和光滑性约束，构建所述目标地形区域表面的辐射模型：所述光滑性约束为：式中，i(x',y')表示目标地形区域(x',y')处在0
°
和90
°
下的偏振图像的光源强度之和，ρ表示预设的偏振图像反照率，p
x
,p
y
,q
x
,q
y
分别为p(x',y')、q(x',y')对应于x方向和y方向上的导数；根据所述辐射模型及光滑性约束，计算得到所述目标地形区域表面(x',y')处的先验法向量信息y＝(p(x',y'),q(x',y'),-1)。7.根据权利要求6所述的基于偏振三维成像的遥感地形三维测绘方法，其特征在于，基于所述先验法向量信息对所述第一法向量信息进行修正，得到所述目标地形区域表面修正后的法向量的步骤，包括：按照如下公式计算二元算子按照如下公式计算二元算子式中，t表示非线性算子；基于所述二元算子对第一法向量信息d进行修正：式中，n
cor
表示所述目标地形区域表面的修正法向量信息。8.根据权利要求1所述的基于偏振三维成像的遥感地形三维测绘方法，其特征在于，所
述不同目标区域包括第一目标区域和第二目标区域，所述第一目标区域与所述第二目标区域至少存在部分重叠区域；根据所述修正法向量信息和所述预设角度下不同目标地形区域的偏振图像，进行目标区域的三维场景重建的步骤，包括：从所述第一目标区域的任意一幅偏振图像中选取一特征点，所述特征点位于重叠区域内；获取所述特征点在所述第二目标区域中的位置信息，并建立所述第一目标区域的偏振图像与所述第二目标区域的偏振图像之间的变换关系；根据所述变换关系拼接所述第一目标区域的偏振图像与所述第二目标区域的偏振图像；根据拼接后的图像、第一目标区域表面的修正法向量信息和第二目标区域表面的修正法向量信息进行三维重建。9.根据权利要求8所述的基于偏振三维成像的遥感地形三维测绘方法，其特征在于，根据拼接后的图像、第一目标区域表面的修正法向量信息和第二目标区域表面的修正法向量信息，按照如下公式进行三维重建：式中，f{}和f-1
{}分别表示离散傅里叶变换和傅里叶逆变换，m、n分别表示拼接后的图像中像素点的行数及列数，u、ν分别表示离散微算子在x和y方向上的傅里叶系数。

技术总结
本发明公开了一种基于偏振三维成像的遥感地形三维测绘方法，包括：采集不同目标地形区域的偏振图像；计算目标地形区域表面的偏振度以及入射光的入射角、方位角；构建法向量与入射光的入射角、方位角之间的第一映射关系，获得第一法向量信息；计算太阳的天顶角与方位角，构建光源矢量与天顶角、方位角之间的第二映射关系，获得目标地形区域表面的光源矢量；构建目标地形区域表面的辐射模型，并计算先验法向量信息；基于先验法向量信息对第一法向量信息进行修正，得到修正法向量信息；根据修正法向量信息和预设角度下不同目标地形区域的偏振图像，进行目标区域的三维场景重建。本发明可以实现大场景下的卫星遥感地形三维重建，提高了运算速度。提高了运算速度。提高了运算速度。


技术研发人员：李轩 祝锦涛 邵晓鹏 蔡玉栋 刘志强 吴秉松 宋家伟
受保护的技术使用者：西安电子科技大学杭州研究院
技术研发日：2023.04.28
技术公布日：2023/8/13