[0196] 其次,确定fBm的其他参数,
[0197] 当H≤0.75时,
[0198]
[0199] 当H>0.75时,
[0200]
[0201] 其中,H为Hurst参数,β、α、c1、c2为中间参数,无实际的数学意义,只是为了方便进行后续的公式描述。在本实施例中,H=0.7,对应的β=0,α=0.7,c1=0.3,c2=0.7。这4个参数将陆续在步骤S1.3.1.3中得到具体使用。
[0202] 步骤S1.3.1.3,使用步骤S1.3.1.2中的参数计算单块矩阵。单块矩阵是循环分块矩阵的基本构成,是步骤S1.3.1.4的输入参数。
[0203] 计算单块矩阵的目的是提供具有距离自相似性的循环分块矩阵。具体过程描述如下。
[0204] 引入各向同性函数r,其表述为
[0205]
[0206] 其中,其中(x,y)为待仿真的分块矩阵中每个元素的坐标,(x0,y0)为待仿真的单块矩阵中同性函数的起算点坐标,即(0,0)。待仿真单块矩阵中每个元素的距离权重函数表述为
[0207]
[0208] 其中c1、c2、α以及β为步骤S1.3.1.2中确定的距离权重参数。
[0209] 步骤S1.3.1.4,计算循环分块矩阵。
[0210] 将距离权重函数进行左右镜像和上下镜像,扩展之后的结果即为循环分块矩阵。扩展结果见下式右侧。
[0211]
[0212] 其中,ρsl为距离权重矩阵的第s列第l行的元素,S为距离权重矩阵的总列数,L为距离权重矩阵的总行数。在本实施例中,S=4000,L=4000。最终的循环分块矩阵为8000×8000的矩阵。
[0213] 步骤S1.3.2,计算与所述循环分块矩阵有相同协方差的随机场。
[0214] 这一步骤是为了求解公式(19)中的单位矩阵Q,使得单位矩阵Q满足公式(17)的定义。
[0215] 步骤S1.3.2具体包括以下子步骤:
[0216] 步骤S1.3.2.1,计算循环分块矩阵的特征值。
[0217] 循环分块矩阵的特征值设定为λ,根据特征值的定义,λ满足如下公式[0218]
[0219] 其中E为单位矩阵,C为循环分块矩阵,解算上述方程中的λ,即可得到矩阵的特征值,对于本实施例来说,矩阵具备8000个特征值,表示为(λ1 λ2 …… λ8000),特征值为步骤S1.3.2中的公式(17)中的矩阵Q的创立提供前提条件。
[0220] 步骤S1.3.2.2,生成二维复随机场。
[0221] 二维复随机场Fcr中的每个元素的表达式为:
[0222] Fcr(j,k)=a+bi (27)
[0223] 其中,Fcr(j,k)为二维复随机场中第j列第k行的元素,a为元素的实部,b为元素的虚部, 在本实施例中,随机数矩阵大小为(8000×8000)。因此实部和虚部均可组成一个8000×8000的矩阵。
[0224] 步骤S1.3.2.3,将步骤S1.3.2.1获取的特征值构成特征值矩阵,并与步骤S1.3.2.2生成的所述二维复随机场相乘,计算得到与所述循环分块矩阵有相同协方差的随机场。
[0225] 将步骤S1.3.2.1获取的特征值构成特征值矩阵,并与所述二维复随机场相乘,得到具有循环分块矩阵特征的随机场,即
[0226]
[0227] 其中,Fr为与循环分块矩阵有相同协方差的随机场,Fcr为步骤S1.3.2.1中生成的二维复随机场,(λ1 λ2 …… λ8000)为循环分块矩阵的特征值。
[0228] 步骤S1.3.2.4,对S1.3.2.3获取的随机场Fr进行傅里叶变换,并抽取傅里叶变换后的随机场Fr的正频部分的数据。
[0229] 在复数域进行FFT变换,其表达式为
[0230] F(j,k)=FFT(Fr)(j,k),j=1,2,...,8000,k=1,2,...,8000 (29)
[0231] 其中,F(j,k)为傅里叶变换之后,矩阵第j列第k行的元素,FFT为快速傅里叶变换函数,Fr为与循环分块矩阵有相同协方差的随机场,来自于步骤S1.3.2.2。将以上公式的矩阵形式拆分为元素形式,那么每个元素(j,k)的变换结果可表达为:
[0232]
[0233] 其中,F(j,k)为傅里叶变换之后,矩阵第j列第k行的元素,Fr为与循环分块矩阵有相同协方差的随机场,来自于步骤S1.3.2.2。M为Fr的矩阵大小,m与n是积分自变量,exp[]为指数函数, 对于本实施例来说,M=8000,对于上面的等式来说,傅里叶变换前后的数值均为复数,变换之后的矩阵,左下角为正频部分,右上角为负频部分,两者的
频谱特征相同。我们只取左下角正频部分的数据进行后续的仿真分析。
[0234] 步骤S1.3.2.5,使用S1.3.2.4中得到的傅里叶变换的正频数据解算独立随机场,得到与所述正频数据的实部和虚部分别对应的两个独立随机场;
[0235] 在复数的FFT变换之后,实频空间与复频空间是
正交的,两者形成的随机场相互独立,因此S1.3.2.4中的正频部分数据的实部与虚部,分别构成两个独立随机场。
[0236] 步骤S1.3.2.6,将步骤S1.3.2.5中得到的所述两个独立随机场的数据进行拉伸变换,得到拉伸变换后的随机场。
[0237] 步骤S1.3.2.5获取独立随机场数据是
频率域数据,值域范围极大,在计算过程中易造成数据溢出。因此需要将数据拉伸到[0,1]之间。步骤S1.3.2.6具体步骤包括:
[0238] S1.3.2.6.1,减掉所述独立随机场原点处的值,确保所述独立随机场的原点值为0,以符合步骤S1.3中描述的fBm的第一个特性。
[0239] S1.3.2.6.2,对所述两个独立随机场的数据分别进行数据拉伸。
[0240] 拉伸公式为:
[0241]
[0242] 其中,F′(j,k)为拉伸之后第j列第k行的值,F(e,k)为拉伸之前第j列第k行的值,Fmin为S1.3.2.6.1中经修正的独立随机场矩阵的最小值,Fmax为矩阵的最大值。
[0243] 步骤S1.3.3,对在S1.3.2中获取的随机场进行改正,得到最终的SAR影像的外接DEM。
[0244] 步骤S1.3.3具体包括以下子步骤:
[0245] 步骤S1.3.3.1,为所述随机场添加增量平稳特性。
[0246] 随机场改正的目的是,使得fBm在噪声分形的基础上,具备增量平稳性。随机场改正的方法,是在随机场中的每个像素值中,添加具备增量平稳特性的数值,此数值可表达为:
[0247]
[0248] 其中,G为增量平稳矩阵,G(j,k)为增量平稳矩阵G的第j列第k行的元素,n1与n2分别为两个随机数,txj为步骤S1.3.1.1中待仿真的单块矩阵行列坐标归一化后的第j列的x坐标,tyk为步骤S1.3.1.1中待仿真的单块矩阵行列坐标归一化后的第k列的y坐标,c2为步骤S1.3.1.2中定义的距离权重矩阵参数。
[0249] 步骤S1.3.3.2,将步骤S1.3.3.1获得的具有增量平稳特性的随机场的数据进行拉伸变换。具体的,将所述数据拉伸到[0,1]。
[0250] 拉伸过程与步骤S1.3.2.6相同。此处不再赘述。
[0251] 步骤S1.3.3.3,将步骤S1.3.2.5获得的随机场与步骤S1.3.3.2获得的随机场相加。
[0252] 步骤S1.3.3.4,将步骤S1.3.3.3获得的随机场拉伸到指定高程范围。
[0253] 由步骤S1.2.5可知,DEM高程范围为3500–6000m。因此,将fBm结果拉伸到3500–6000m。即可得到主影像的外接DEM。
[0254] 优选的,选用虚部生成的fBm结果进行后续仿真。
[0255] 仿真结果如图2(a)和图2(b)所示。其中图2(a)为步骤S1.3.2.5中得到的实部随机场,经步骤S1.3.3改正之后得到的外接DEM结果,图2(b)为步骤S1.3.2.5中得到的虚部随机场,经步骤S1.3.3改正之后得到的外接DEM结果。
[0256] 步骤S2,使用在步骤S1仿真的所述SAR影像的外接DEM仿真所述SAR影像的相位。
[0257] 按照干涉的理论,主影像与从影像的
相位差构成干涉相位。因此,为了完成主从影像的相位仿真,需要首先完成主从影像的干涉相位仿真,仿真过程中需要用到步骤S1.1中的主、从影像的成像几何参数,以及步骤S1.3中的仿真的外接DEM。随后完成从影像相位的仿真。最终完成主影像相位的仿真。
[0258] 各步骤详细如下。
[0259] 步骤S2.1,使用所述SAR影像的外接DEM,仿真所述主影像和所述从影像的干涉相位。
[0260] 计算过程中用到的干涉几何如图3所示。后续的讨论将基于此干涉几何进行。其中,S1为主影像对应的天线相位中心位置,S2为从影像对应的天线相位中心位置。B为基线长度,α为基线倾角,θ为侧视角,ρ为基线矢量与LOS向的夹角,B⊥为垂直基线,r和r+△r分别为主从影像天线相位中心到地面点P的距离,H为卫星高度,RH为卫星到地心的距离,Rh为地面点到地心的距离,Re为地球
曲率半径,h为地面点高程。本图是雷达干涉的基本原理图。
[0261] 步骤S2.1具体包括以下子步骤:
[0262] 步骤S2.1.1,将步骤S1仿真得到的DEM转换到SAR影像空间。
[0263] 采用步骤S1.2.1、S1.2.2以及S1.2.3中的R-D模型,将步骤S2得到的外接DEM中每个坐标点投影到SAR影像像素空间中,并保留像素坐标(x,y)满足1≤x≤S,且1≤y≤L的点,其中S为仿真得到的DEM的总列数,L为总行数。投影之后的高程图如图4所示,其中图4(a)为步骤S1.3仿真得到的外接DEM,黑色框线为主影像对应的影像范围,图4(b)为转换到SAR影像空间的DEM。
[0264] 步骤S2.1.2,使用通过步骤S2.1.1转换后的DEM,计算SAR影像空间中每个像素点的垂直基线和平行基线。
[0265] 干涉相位主要包括两类,即平地效应相位,以及地形相位。其中平地效应相位由平行基线表达,地形相位由垂直基线表达。因此本步骤需解算这两个基线参数。本步骤的输入值为步骤S2.1.1转换的DEM,输出值用以进行步骤S2.1.3的相位解算。
[0266] 步骤S2.1.2具体包括以下子步骤:
[0267] 步骤S2.1.2.1,计算每个像素点的视线向向量。
[0268] 每个像素点的视线向的向量表达如下:
[0269]
[0270] 其中, 为像素坐标(x,y)处的地物对应的视线向向量, 为主影像的相位中心坐标矢量,可通过步骤S1.2.1中卫星的三维坐标的解算方法解算获得。 为每个像素点的坐标矢量,坐标矢量使用R-D模型解算获得,解算过程参见步骤S1.2.1-S1.2.4,其中的高程信息由步骤S2.1.1转换的DEM提供。
[0271] 步骤S2.1.2.2,计算每个像素点垂直视线向的向量
[0272]
[0273] 其中,为卫星飞行方向矢量,采用卫星位置矢量 与步骤S2.1.2.1中解算得到的像素点视线向向量 解算获得,即
[0274]
[0275] 卫星位置矢量的解算过程参见步骤S1.2.1。
[0276] 步骤S2.1.2.3,计算SAR影像空间中每个像素点的基线矢量。
[0277] 步骤S2.1.2.3具体包括以下子步骤:
[0278] 步骤S2.1.2.3.1,计算每个点的像素矢量 在从影像中的像素坐标。
[0279] 计算过程中,使用的是步骤S1.2.1中的斜距方程,以及步骤S1.2.2中的多普勒方程。即使用两个方程解算两个未知数。解算过程中使用最小二乘或者
牛顿
迭代法均可,此处不再赘述。
[0280] 步骤S2.1.2.3.2,计算从影像中像素坐标对应的位置矢量
[0281] 卫星位置矢量的解算过程参见步骤S1.2.1。
[0282] 步骤S2.1.2.3.3,计算每个像素点的基线矢量
[0283]
[0284] 其中, 为主影像中卫星的位置矢量, 为从影像中卫星的位置矢量。两卫星矢量均可通过步骤S1.2.1中卫星的三维坐标的解算方法解算获得。
[0285] 步骤S2.1.2.4,根据每个像素点的基线矢量,计算每个像素点的垂直基线和平行基线。
[0286] 每个像素点的垂直基线B⊥以及平行基线B||可通过下式进行计算:
[0287]
[0288] 其中, 为步骤S2.1.2.1计算得到的像素点(x,y)的视线向量, 为步骤S2.1.2.2计算得到的像素点垂直视线向的向量。
[0289] 步骤S2.1.3,使用步骤S2.1.2得到的垂直基线和平行基线计算干涉相位。
[0290] 步骤S2.1.3具体包括以下子步骤:
[0291] 步骤S2.1.3.1,计算从影像每个像素点的斜距R2。
[0292]
[0293] 其中,ρ为基线矢量与雷达视线向的夹角,表达为 cos-1为反余弦函数, 为基线向量的单位向量,即 为步骤S2.1.2.3中计算的基线向量,为基线向量的模长, 为像素点视线向量的单位向量,即 为
步骤步骤S2.1.2.1中计算的像素点的视线向量,(x,y)为像素点的像素坐标, 为视线向量的模长。R1为主影像中的斜距,表达为R1=R0+x△R,R0为主影像近地点斜距,x为像素点的x方向坐标,△R为影像的斜距分辨率,近地点斜距和斜距分辨率在表1中均有描述。
[0294] 步骤S2.1.3.2,计算每个像素点的总相位
[0295]
[0296] 其中,λ为雷达波长,本实施例采用的是X波段数据,其波长为3cm,其中R1为主影像中的斜距,表达为R1=R0+x△R,其中R0为主影像近地点斜距,x为像素点的x方向坐标,△R为影像的斜距分辨率,近地点斜距和斜距分辨率在表1中均有描述。R2为步骤S2.1.3.1中解算得到的从影像每个像素点的斜距。本实施例中,计算得到的总相位如图5(a)所示。
[0297] 步骤S2.1.3.3,计算每个像素点的地形相位
[0298] 地形相位的表达式如下:
[0299]
[0300] 其中, 为步骤S2.1.3.2中获得的每个像素点的总相位,λ为雷达波长,本实施例采用的是X波段数据,其波长为3cm,B||为步骤S2.1.2.4解算得到的每个像素点的平行基线。计算得到的地形相位如图5(b)所示。
[0301] 步骤S2.2,使用随机噪声模型仿真从影像的相位。
[0302] 从影像的相位完全由随机噪声组成。依据干涉相位噪声与相干性的关系可知[0303]
[0304] 其中, 为相位噪声 的方差,γ为相干性,地形测绘过程中,一般令γ>0.9。NL为从影像的多视数,多视数可降低影像噪声,在仿真过程中,为了保持相位的真实性,一般令多视数为1。这种情况下,可解得相位噪声为
[0305]
[0306] 依据误差传播定律,单景影像的相位噪声 表达为
[0307]
[0308] 其中, 为公式42得到的相位噪声。
[0309] 因此,生成均值为0,标准差为0.242的随机噪声,其大小为8000×8000,作为从影像相位。即
[0310]
[0311] 其中 为从影像相位, 为服从正态分布的相位,N(0,0.242)代表均值为0,方差为0.242的标准正态分布,结果如图6(a)所示。
[0312] 步骤S2.3,使用步骤S2.1得到的所述从影像的干涉相位和步骤S2.2得到的所述从影像的相位,计算所述主影像的相位。
[0313] 主影像相位的表达式为
[0314]
[0315] 其中, 为步骤S2.2获得的从影像相位, 为主影像像素空间中每个像素点的干涉相位,由步骤S2.1.3.2获得, 为噪声相位,由其生成过程与步骤S2.2相同,不再赘述。仿真得到的主影像相位如图6(b)所示。
[0316] 步骤S3,使用在步骤S1仿真的所述SAR影像的外接DEM仿真所述SAR影像的振幅。
[0317] SAR影像为侧视成像,影像中存在叠掩和阴影。因此影像由三类主要特征构成,第一类即平地和非平地像素构成的一般像素,第二类为叠掩区的像素,第三类为阴影区的像素。计算过程中,会使用到所述主影像以及从影像像素空间中每个像素点的本地法线向量、像平面法线、名义入射角、本地入射角等。具体步骤如下。
[0318] 步骤S3.1,仿真所述主影像的振幅。
[0319] 步骤S3.1.1,计算所述主影像像素空间每个像素点的本地法线向量。
[0320] 本地法线向量 使用如下公式计算:
[0321]
[0322] 其中, 为像素点(1,y)对应的视线向向量, 为像素点(0,y)对应的视线向向量, 为像素点(1,y-1)对应的视线向向量, 代表向量的叉乘。公式中的三个视线向向量均可由步骤S2.1.2.1获得。本步骤为步骤S3.1.2提供本地法线向量的单位向量,进而判断某一像素是否隶属于平地、非平地、叠掩、阴影区域。
[0323] 步骤S3.1.2,基于所述主影像像素空间中像素点的本地法线向量,判断所述像素点是否是平地像素。
[0324] 本地法线向量与本地坐标向量的夹角可表示为χ,其表达式为
[0325]
[0326] 其中,cos-1为反余弦函数, 为本地法线向量的单位向量,其表达式为[0327]
[0328] 为步骤S3.1.1计算得到的每个像素点的本地法线向量, 为本地法线向量的模长, 为像素点坐标矢量的单位向量,其表达式与上式相同,需采用每个限速点的坐标矢量解算得来,坐标矢量使用R-D模型解算获得,解算过程参见步骤S1.2.1-S1.2.4,其中的高程信息由步骤S2.1.1转换的DEM提供。当χ小于一定角度时,可认为是平地像素(如图7所示)。优选的,令χ<1°时为平地像素。
[0329] 步骤S3.1.3,计算每个像素点的名义入射角β。
[0330]
[0331] 式中,各参数的含义与步骤S2.1相同,cos-1为反余弦函数,RH为卫星到地心的距离,R1为主影像斜距,表达为R1=R0+x△R,其中R0为主影像近地点斜距,x为像素点的x方向坐标,△R为影像的斜距分辨率,近地点斜距和斜距分辨率在表1中均有描述。Re为地球曲率半径,h为地面点高程。
[0332] 步骤S3.1.4,基于对所述像素点是否为平地像素的判断和像素点的名义入射角,计算每个像素点的本地入射角θ。
[0333]
[0334] 其中,β为步骤S3.1.3中得到的名义入射角,cos-1为反余弦函数, 为卫星速度的单位矢量,由步骤S1.2.2.1获得, 为本地法线向量的单位矢量,由步骤S3.1.3获得。
[0335] 步骤S3.1.5,基于对所述像素点是否为平地像素的判断和像素点的本地入射角,计算每个像素点的振幅值Am。
[0336] 公式表达式如下:
[0337]
[0338] 其中,θ为步骤S3.1.4中得到的每个像素点的本地入射角,μ为像素尺寸归一化参数,表达式如下:
[0339]
[0340] 其中,χ为本地法线向量与本地坐标向量的夹角,θ为步骤S3.1.4中得到的每个像素点的本地入射角
[0341] 至此,本发明已经完成了主影像像素空间中大多数像素的振幅计算。然而,SAR影像为侧视成像,在成像过程中会生成特有的叠掩、阴影等现象。因此需要在下面的步骤中计算叠掩区和阴影区的振幅值。
[0342] 步骤S3.1.6,计算叠掩区振幅值。
[0343] 步骤S3.1.6.1,确定叠掩点。
[0344] 计算每个点的本地坡向角ψ,其表达式为
[0345]
[0346] 其中,cos-1为反余弦函数, 为本地法线向量的单位向量,通过公式(48)计算得来, 为像平面法线的单位向量,将ψ大于90°的定义为叠掩点。像平面法线的单位向量表达为
[0347]
[0348] 其中,|Nimg|为像平面法线的模长,Nimg为像平面法线,其表达式为[0349]
[0350] 为主影像成像时刻的卫星速度矢量,其解算方法参见步骤S1.2.2.1。 为像素点(1,y-1)对应的视线向向量。
[0351] 步骤S3.1.6.2,计算叠掩点振幅值。
[0352] 寻找叠掩点左侧的最近有效像素及其振幅值(xl,yl,Aml),右侧的最近有效像素及其振幅值(xr,yr,Amr),使用线性插值获取此点(x,y)的振幅Am,即
[0353]
[0354] 步骤S3.1.7,计算阴影区灰度值。
[0355] 步骤S3.1.7.1,确定阴影点。
[0356] 本地入射角大于90°的像素点定义为阴影点。本地入射角由步骤S3.1.4计算获得。
[0357] 步骤S3.1.7.2,计算阴影点振幅值。
[0358] 寻找阴影点左侧的最近有效像素及其振幅值(xl,yl,Aml),右侧的最近有效像素及其振幅值(xr,yr,Amr),使用线性插值获取此点(x,y)的振幅Am,即
[0359]
[0360] 步骤S3.2,仿真所述从影像的振幅。
[0361] 从影像振幅仿真过程中,用到了表2中提供的参数。具体步骤与S3.1类似,此处不再赘述。图8给出了使用主影像以及从影像的仿真结果。
[0362] 步骤S4,基于步骤S2仿真得到的SAR影像的相位和步骤S3仿真得到的SAR影像的振幅,仿真所述SAR影像。
[0363] 步骤S4.1,仿真所述主影像。
[0364] 对主影像每个像素来说,其像素值表达为一个复数,即
[0365] F(x,y,m)=am+bmi (58)
[0366] 其中,F(x,y,m)为主影像每个像素的像素值,(x,y)为主影像像素的坐标,m为主影像像素的复数值, Am为步骤S3.1仿真的主影像振幅, 为步骤2.3仿真的主影像的相位。
[0367] 步骤S4.2,仿真所述从影像。
[0368] 对从影像每个像素来说,其像素值表达为一个复数,即
[0369] F(x,y,s)=as+bsi (59)
[0370] 其中,F(x,y,s)为从影像每个像素的像素值,(x,y)为为从影像像素的坐标,s为从影像像素的复数值, As为步骤S3.2仿真的从影像振幅, 为步骤2.2仿真的从影像相位。
[0371] 图9给出了相位和振幅叠加的结果,其中图9(a)为主影像结果,图9(b)为从影像结果。
[0372] 现有技术多采用SRTM数据或已知的DEM数据进行仿真,仿真过程中无法从纯数学角度探讨地形测绘过程中的误差来源,本发明提出了基于fBm的仿真方法,一方面,fBm能够表达数据的自相似性,这与地形本身一定范围内的相似性一致;另一方面,fBm能够表达数据的随机性,这与地形本身极小范围内的随机起伏一致。因此,fBm的数学特征能够很好的描述地形信息,这使得从理论层面进行地形特性的研究成为可能。此外,传统仿真过程中,参数为仿真系统生成,而非真实参数,本发明仿真过程中采用的参数为真实的TanDEM-X卫星干涉参数,参数具有极好的工程应用性,能够确保在干涉处理过程中,所有的误差来源不会与干涉几何状态有太大关联,从而使得研究人员能够基于本发明得到更为精确的误差仿真结论。本发明为基于干涉SAR技术的卫星测绘提供了良好的参考。
[0373] 以上内容仅为本发明的较佳实施例,对于本领域的普通技术人员,依据本发明的思想,在具体实施方式及应用范围上均会有改变之处,本
说明书内容不应理解为对本发明的限制。