CN113900069A - 一种基于干涉成像高度计的垂线偏差计算方法及其系统 - Google Patents
一种基于干涉成像高度计的垂线偏差计算方法及其系统 Download PDFInfo
- Publication number
- CN113900069A CN113900069A CN202111153467.6A CN202111153467A CN113900069A CN 113900069 A CN113900069 A CN 113900069A CN 202111153467 A CN202111153467 A CN 202111153467A CN 113900069 A CN113900069 A CN 113900069A
- Authority
- CN
- China
- Prior art keywords
- point
- sea surface
- data
- longitude
- resampling
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 48
- 238000003384 imaging method Methods 0.000 title claims abstract description 26
- 238000012952 Resampling Methods 0.000 claims abstract description 157
- 238000000034 method Methods 0.000 claims abstract description 53
- 230000002159 abnormal effect Effects 0.000 claims abstract description 25
- 238000005070 sampling Methods 0.000 claims abstract description 18
- 238000001914 filtration Methods 0.000 claims abstract description 12
- 230000008030 elimination Effects 0.000 claims abstract description 4
- 238000003379 elimination reaction Methods 0.000 claims abstract description 4
- 239000000126 substance Substances 0.000 claims description 27
- 238000004590 computer program Methods 0.000 claims description 9
- 238000009499 grossing Methods 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 3
- 238000013507 mapping Methods 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 7
- 238000005259 measurement Methods 0.000 description 7
- 230000005484 gravity Effects 0.000 description 5
- 238000007796 conventional method Methods 0.000 description 2
- OMPJBNCRMGITSC-UHFFFAOYSA-N Benzoylperoxide Chemical compound C=1C=CC=CC=1C(=O)OOC(=O)C1=CC=CC=C1 OMPJBNCRMGITSC-UHFFFAOYSA-N 0.000 description 1
- 244000178993 Brassica juncea Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 1
- 239000010931 gold Substances 0.000 description 1
- 229910052737 gold Inorganic materials 0.000 description 1
- 208000028755 loss of height Diseases 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 229940060894 topex Drugs 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/40—Means for monitoring or calibrating
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C13/00—Surveying specially adapted to open water, e.g. sea, lake, river or canal
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C5/00—Measuring height; Measuring distances transverse to line of sight; Levelling between separated points; Surveyors' levels
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/02—Systems using reflection of radio waves, e.g. primary radar systems; Analogous systems
- G01S13/06—Systems determining position data of a target
- G01S13/08—Systems for measuring distance only
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/882—Radar or analogous systems specially adapted for specific applications for altimeters
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/9021—SAR image post-processing techniques
- G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
Abstract
本发明属于海洋测绘技术领域，具体地说，涉及一种基于干涉成像高度计的垂线偏差计算方法，该方法包括：步骤1)对原始数据中的每个原始采样点进行重采样，同时计算每个重采样点的海面高度误差；步骤2)对重采样数据中的每个重采样点的海面高度进行异常值剔除和低通滤波，得到多个重采样点的平滑海面高度；步骤3)计算的每个重采样点的平滑数据；步骤4)将每个重采样点的平滑数据作为一个数据点，合并组成数据集；步骤5)对数据集进行空间划分，得到多个经纬度网格点；步骤6)对每个经纬度网格点，采用小邻域窗口方法，得到每个经纬度网格点的挑选的所有数据点；步骤7)利用最小二乘最优法则，计算每个经纬网格的垂线偏差分量。
Description
技术领域
本发明属于海洋测绘、卫星高度计测高数据处理、海洋大地测量、遥感信息提取及海洋参数反演技术领域，具体地说，涉及一种基于干涉成像高度计的垂线偏差计算方法及其系统。
背景技术
卫星高度计是测量海表面高度的仪器，通常搭载在卫星平台上。干涉成像高度计测高的原理是利用双天线的相位差转化成距离差，从而测得海表面高度的具体数值。并且宽刈幅的特性使干涉成像高度计能够同时测量大范围的海域，通常可达到百公里级，大大改善了传统高度计在交轨方向上测高数据的缺失，对于海洋高度测量的空间分辨率提升到5-10km，精度保证在3-5cm，对海洋中尺度和亚中尺度过程海洋现象的反演提供了可靠的数据支撑。
海洋垂线偏差的定义是大地水准面上某点的重力方向与通过该点到正常椭球面的法线方向之间的夹角。海洋垂线偏差是计算海洋重力异常乃至计算海洋深度的中间量，对于物理海洋学、地球物理学、地球动力学、海洋动力学、海洋地质学等基础研究至关重要，同时在资源勘探、海啸预警、海底电缆铺设、惯性导航等民用或军用领域也发挥着作用。
海洋垂线偏差通常描述为北向(子午向)分量ξ和东向(卯酉向)分量η，星载雷达高度计的主要观测量包括：海表面高度、有效波高和海面风速。其中，海表面高度包含大地水准面高和海面地形两部分信息；其中，大地水准面高不仅可以直接反演海洋重力异常，还可用于获取垂线偏差信息。
早在20世纪80年代，国外学者最早提出基于原始海面高观测值计算沿轨垂线偏差的想法，相邻数据的差分过程能够削弱径向轨道误差等长波长误差项的影响，从而有效地揭示海洋大地水准面的短波起伏(文献[1]：D.T.Sandwell,“A Detailed View of theSouth-Pacific Geoid from Satellite Altimetry,”J.Geophys.Res.,vol.89,no.B2,pp.1089–1104,Feb.1984,doi:10.1029/JB089iB02p01089.)。1992年，Sandwell学者提出了一种交叉点处垂线偏差分量的严密计算方法，需要计算大地水准面对时间的导数以及星下点轨迹处的纬向和径向速率，该方法主要受限于交叉点的分布(文献[2]： D.T.Sandwell,“Antarctic marine gravity field from high-density satellite altimetry,”Geophysical Journal International,vol.109,no.2,pp.437–448,May 1992)。1995年，Olgiati等学者在交叉点处计算与地面轨迹正交方向的垂线偏差，然后采用Akima样条插值法，推算每个数据点处的轨迹正交方向的垂线偏差，最后，联立沿轨垂线偏差信息，计算每个数据点处的垂线偏差方向分量，从而解除了对于交叉点分布的限制，但精度严重依赖于插值算法(文献[3]：A.Olgiati,G.Balmino,M.Sarrailh,and C.M. Green,“Gravityanomalies from satellite altimetry:comparison between computation via geoidheights and via deflections of the vertical,”Bulletin Géodésique,vol.69,no.4,pp. 252–260,Dec.1995)。1998年，中国台湾学者黄金维提出一种基于沿轨垂线偏差和方位角信息最小二乘解算垂线偏差分量的网格化方法，称为沿轨加权最小二乘法，该方法可以有效地融合使用多代卫星测高资料(文献[4]：C.Hwang,E.-C.Kao,and B.Parsons,“Global derivation of marine gravity anomalies from Seasat,Geosat,ERS-1andTOPEX/POSEIDON altimeter data,”Geophysical Journal International,vol.134,no.2,pp. 449–459,Aug.1998)。
然而，目前的卫星高度计垂线偏差计算方法都是针对星下点观测方式的传统型高度计，计算方法仅适用于单条测线。对于干涉成像高度计宽刈幅的相位测量模式，现有的方法仅从单一方向坡度计算垂线偏差，会大大增加误差。同时，也没有考虑精确度低、误差噪声大的远幅侧和近幅侧对垂线偏差精度的影响。
发明内容
为解决现有技术存在的上述缺陷，本发明提出了一种基于干涉成像高度计的垂线偏差计算方法，该方法包括：
步骤1)对原始数据中的每个原始采样点进行重采样，得到由多个重采样点组成的重采样数据，同时计算每个重采样点的海面高度误差；
步骤2)对重采样数据中的每个重采样点的海面高度进行异常值剔除和低通滤波，得到多个重采样点的平滑海面高度；
步骤3)分别沿卫星轨道飞行方向、垂直于卫星轨道方向和倾斜于卫星轨道方向，根据步骤1)得到的每个重采样点的海面高度误差和步骤2)得到的每个重采样点的平滑海面高度，计算的每个重采样点的平滑数据；其中，上述平滑数据包括：每个重采样点的海面高度的海面坡度、海面方位角和海面坡度误差；
步骤4)将每个重采样点的平滑数据作为一个数据点，合并组成数据集；
步骤5)根据实际需求，确定均匀分布的经纬度网格的尺寸，对数据集进行空间划分，得到多个经纬度网格点；其中，每个经纬度网格点内包括多个数据点；
步骤6)对每个经纬度网格点，采用小邻域窗口方法，挑选出小邻域窗口内的所有数据点，得到每个经纬度网格点的挑选的所有数据点；
步骤7)利用最小二乘最优法则，使用每个经纬度网格点的挑选的所有数据点，计算每个经纬网格的垂线偏差分量。
作为上述技术方案的改进之一，所述步骤1)具体包括：
将原始数据中的每个50m×50m空间分辨率的原始采样点进行重采样，得到由多个2km×2km空间分辨率的重采样点组成的重采样数据，包括经度、纬度、海面高度，并在重采样的过程中，计算每个重采样点的海面高度误差：
重采样点邻域窗口尺寸范围为1km×1km，在重采样点邻域窗口内计算邻域窗口尺寸范围内的所有原始采样点的海面高度的标准差，作为该采样点的海面高度误差σSSH：
σSSH i＝std(SSHnear)
其中，σSSH i表示第i个重采样点的海面高度误差，SSHnear表示该重采样点对应的邻域窗口尺寸范围内的所有原始采样点的海面高度，符号std表示计算标准差；
将第i个重采样点的经度、纬度、海面高度、海面高度误差组成卫星高度计重采样后的第i个重采样数据。
作为上述技术方案的改进之一，所述步骤2)具体包括：
根据卫星高度计重采样后的所有重采样数据中的每个重采样点的海面高度SSH的经度和纬度组成的位置信息，在已发布的平均海面高数据集中，插值出对应点的平均海面高MSS；重采样后的所有海面高度与对应位置的平均海面高度做差，得到二者的差值，并计算该差值的标准差，记为σall；
判断每个重采样点的海面高度是否为异常值：
对第i个重采样点的海面高度SSHi与对应位置的平均海面高度MSSi做差，得到第i个重采样点的差值，若该第i个重采样点的差值的绝对值大于3σall，则认为该重采样点的海面高度是异常值，并删除该重采样点的海面高度；
|SSHi-MSSi|＞3×σall
若该第i个重采样点的差值的绝对值小于或等于3σall，则认为该重采样点的海面高度是正常值，保留该重采样点的海面高度；
根据判断结果，得到剔除异常值后的重采样数据点；
重复上述过程，对重采样数据中的每个重采样点的海面高度进行异常值剔除，得到由多个剔除异常值后的重采样数据点组成的剔除后的数据集；
采用高斯加权滤波器，对剔除异常值后的数据集中的重采样数据点中的海面高度进行低通滤波，降低高频噪声影响，得到滤波后的海面高度，记为多个重采样点的平滑海面高度。
作为上述技术方案的改进之一，所述步骤3)具体包括：
针对沿卫星轨道飞行方向：
假设步骤2)中得到的每个重采样点的平滑海面高度中对应的重采样点为数据点，将沿卫星轨道飞行方向的两个数据点的平滑海面高度做差，将该差值除以两个相邻数据点之间的真实距离，得到两个相邻数据点之间的沿卫星轨道飞行方向的海面坡度SSSn：
其中，Hi,Hj分别表示两个相邻的数据点中的第i个数据点和第j个数据点的沿卫星轨道飞行方向的平滑海面高度；sij表示两个相邻数据点间的真实距离；
在计算两个相邻数据点之间的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向的方位角α，记为沿卫星轨道飞行方向的海面方位角：
其中，
是两个相邻的数据点中的第i个数据点的纬度；λi是两个相邻的数据点中的第i个数据点的经度；
是两个相邻的数据点中的第j个数据点的纬度；λj是两个相邻的数据点中的第j个数据点的经度；
此外，计算沿卫星轨道飞行方向的海面坡度误差σSSS：
其中，
是两个相邻数据点中的第i个数据点的海面高度误差；
分别是两个相邻数据点中的第j个数据点的海面高度误差；
将沿卫星轨道飞行方向的海面坡度SSSn，沿卫星轨道飞行方向的海面方位角α和沿卫星轨道飞行方向的海面坡度误差σSSS作为沿卫星轨道飞行方向的每个重采样点的平滑数据；
针对垂直于卫星轨道方向：
计算两个相邻数据点之间的垂轨向的海面坡度SSSN：
其中，HI,HJ分别表示两个相邻的数据点中的第I个数据点和第J个数据点的垂直于卫星轨道方向的平滑海面高度；sIJ表示两个相邻数据点间的垂直于卫星轨道方向的真实距离；
在计算两个相邻数据点之间的垂轨向的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向垂轨向的海面方位角αN：
其中，
是两个相邻的数据点中的第I个数据点的垂轨向的纬度；λI是两个相邻的数据点中的第I个数据点的垂轨向的经度；
是两个相邻的数据点中的第J个数据点的垂轨向的纬度；λJ是两个相邻的数据点中的第J个数据点的垂轨向的经度；
此外，计算垂轨向海面坡度误差σSSSN：
其中，
是两个相邻数据点中的第I个数据点的垂轨向的海面高度误差；
分别是两个相邻数据点中的第J个数据点的垂轨向的海面高度误差；
将垂轨向的海面坡度SSSN，垂轨向的海面方位角αN和垂轨向海面坡度误差σSSSN作为沿垂直于卫星轨道方向的每个重采样点的平滑数据；
针对倾斜于卫星轨道135°方向：
计算两个相邻数据点之间的倾斜于卫星轨道135°方向的海面坡度SSSM：
其中，HP,HQ分别表示两个相邻的数据点中的第P个数据点和第Q个数据点的倾斜于卫星轨道135°方向的平滑海面高度；sPQ表示两个相邻数据点间的倾斜于卫星轨道135°方向的真实距离；
在计算两个相邻数据点之间的倾斜于卫星轨道135°方向的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向倾斜于卫星轨道135°方向的海面方位角αM：
其中，
是两个相邻的数据点中的第P个数据点的倾斜于卫星轨道135°方向的纬度；λP是两个相邻的数据点中的第P个数据点的倾斜于卫星轨道135°方向的经度；
是两个相邻的数据点中的第Q个数据点的倾斜于卫星轨道135°方向的纬度；λQ是两个相邻的数据点中的第Q个数据点的倾斜于卫星轨道135°方向的经度；
此外，计算倾斜于卫星轨道135°方向海面坡度误差σSSSM：
其中，
是两个相邻数据点中的第P个数据点的倾斜于卫星轨道135°方向的海面高度误差；
分别是两个相邻数据点中的第Q个数据点的倾斜于卫星轨道 135°方向的海面高度误差；
将倾斜于卫星轨道135°方向的海面坡度SSSM，倾斜于卫星轨道135°方向的海面方位角αM和倾斜于卫星轨道135°方向海面坡度误差σSSSM作为沿垂直于卫星轨道135°方向的每个重采样点的平滑数据；
针对倾斜于卫星轨道45°方向：
计算两个相邻数据点之间的倾斜于卫星轨道45°方向的海面坡度SSSO：
其中，Hg,Hh分别表示两个相邻的数据点中的第g个数据点和第h个数据点的倾斜于卫星轨道45°方向的平滑海面高度；sgh表示两个相邻数据点间的倾斜于卫星轨道45°方向的真实距离；
在计算两个相邻数据点之间的倾斜于卫星轨道45°方向的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向倾斜于卫星轨道45°方向的海面方位角αO：
其中，
是两个相邻的数据点中的第g个数据点的倾斜于卫星轨道45°方向的纬度；λg是两个相邻的数据点中的第g个数据点的倾斜于卫星轨道45°方向的经度；
是两个相邻的数据点中的第h个数据点的倾斜于卫星轨道45°方向的纬度；λh是两个相邻的数据点中的第h个数据点的倾斜于卫星轨道45°方向的经度；
此外，计算倾斜于卫星轨道45°方向海面坡度误差σSSSO：
其中，
是两个相邻数据点中的第g个数据点的倾斜于卫星轨道45°方向的海面高度误差；
分别是两个相邻数据点中的第h个数据点的倾斜于卫星轨道45°方向的海面高度误差；
将倾斜于卫星轨道45°方向的海面坡度SSSO，倾斜于卫星轨道45°方向的海面方位角αO和倾斜于卫星轨道45°方向海面坡度误差σSSSO作为沿垂直于卫星轨道 45°方向的每个重采样点的平滑数据。
作为上述技术方案的改进之一，所述步骤5)具体包括：
根据实际需求，确定均匀分布的经纬度网格的尺寸为2′×2′或1′×1′的经纬度网格，按照上述确定的经纬度网格尺寸，对数据集进行空间划分，得到多个经纬度网格点。
作为上述技术方案的改进之一，所述步骤6)具体包括：
挑选出小邻域窗口内的所有数据点，得到挑选的所有数据点；
对于某个经纬度网格点，选择小邻域窗口的尺寸为12′×12′，即选定经纬度范围为网格点的经纬度加减6′，对步骤4)得到的数据集中的每个数据点进行挑选，将经纬度处于小于小邻域窗口尺寸范围内的数据点挑选出来，得到每个经纬度网格的挑选的所有数据点。
作为上述技术方案的改进之一，所述步骤7)具体包括：
利用最小二乘最优法则，使用每个经纬网格的挑选的所有数据点，计算每个经纬网格的垂线偏差分量；其中，每个经纬网格的垂线偏差分量包括：每个经纬网格的垂线偏差北向分量ξ和每个经纬网格的垂线偏差南向分量η；
SSSi+vi＝ξcosαi+ηsinαi,i＝1,2,…,n
其中，SSSi表示第i个数据点的海面坡度；vi表示第i个数据点的残差项，αi是第 i个数据点的海面坡度SSSi的方位角；
基于最小二乘准则和间接平差方法，直接计算该经纬网格处的垂线偏差北向分量ξ和东向分量η；
将上式改写成矩阵形式：
V＝AX-L
其中，
因此，V的表现形式改写为：
假定卫星对地面的观测点都是独立观测，即不同观测点之间不存在互相影响，那么权重P的形式为：
其中，si是两个数据点之间的距离；σSSSZ是第i个数据点的海面坡度误差；该坡度误差是通过步骤3)的四个方向的坡度误差整合得到的；其中，四个方向指的是沿卫星轨道飞行方向、垂轨方向、倾斜于卫星轨道45°方向和倾斜于卫星轨道135°方向。
最终的第m个经纬度网格点的垂线偏差北向分量ξ和东向分量η为：
其中，ξm为第m个经纬度网格点的垂线偏差北向分量；ηm为第m个经纬度网格点的垂线偏差东向分量；
使用第m个经纬度网格点的挑选出的所有数据点，计算得到该经纬度网格点处的垂线偏差分量；
重复上述过程，计算每个经纬度网格点处的垂线偏差分量。
本发明还提供了一种基于干涉成像高度计的垂线偏差计算系统，该系统包括：
高度误差获取模块，用于对原始数据中的每个原始采样点进行重采样，得到由多个重采样点组成的重采样数据，同时计算每个重采样点的海面高度误差；
海面高度获取模块，用于对重采样数据中的每个重采样点的海面高度进行异常值剔除和低通滤波，得到多个重采样点的平滑海面高度；
平滑数据获取模块，用于分别沿卫星轨道飞行方向、垂直于卫星轨道方向和倾斜于卫星轨道方向，根据得到的每个重采样点的海面高度误差和得到的每个重采样点的平滑海面高度，计算的每个重采样点的平滑数据；其中，上述平滑数据包括：每个重采样点的海面高度的坡度、方位角和坡度误差；将每个重采样点的平滑数据作为一个数据点，合并组成数据集；
经纬划分模块，用于根据实际需求，确定均匀分布的经纬度网格的尺寸，对数据集进行空间划分，得到多个经纬度网格点；其中，每个经纬度网格点内包括多个数据点；
挑选模块，用于对每个经纬度网格点，采用小邻域窗口方法，挑选出小邻域窗口内的所有数据点，得到每个经纬度网格点的挑选的所有数据点；和
垂线偏差计算模块，用于利用最小二乘最优法则，使用每个经纬度网格点的挑选的所有数据点，计算每个经纬网格的垂线偏差分量。
本发明还提供了一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现所述的方法。
本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行所述的方法。
本发明与现有技术相比的有益效果是：
1、本发明的方法，在计算坡度时综合考虑了各个方向，降低了单一方向坡度造成的误差，使得垂线偏差计算结果精度，要高于传统星下点高度计；同时考虑了沿轨和垂轨方向误差分布特性，通过不同的权重来指导最小二乘计算，使得结果更加精确；
2、对使用干涉成像机理进行测高的卫星高度计数据，本方法具有普适性；
3、本发明的方法减轻误差噪声大的远幅侧和近幅侧对整体计算结果的影响。
附图说明
图1是本发明的一种基于干涉成像高度计的垂线偏差计算方法的流程图；
图2是图1的一种基于干涉成像高度计的垂线偏差计算方法的实施例的具体方法流程图；
图3是图1的方法转中的步骤3)的沿轨方向的坡度计算示意图；
图4是图1的方法转中的步骤3)的垂轨方向的坡度计算示意图；
图5是图1的方法转中的步骤3)的倾斜轨道135°方向的坡度计算示意图；
图6是图1的方法转中的步骤3)的倾斜轨道45°方向的坡度计算示意图。
具体实施方式
现结合附图和实例对本发明作进一步的描述。
如图1所示，本发明提供了一种基于干涉成像高度计的垂线偏差计算方法，具体涉及一种考虑误差分布的干涉成像高度计垂线偏差计算方法，对于干涉式测高宽刈幅的特性，在计算坡度时综合考虑了各个方向，降低了单一方向坡度造成的误差，使得垂线偏差计算结果精度要高于传统星下点高度计。同时，在计算坡度时，考虑了干涉成像高度计沿轨方向误差分布和垂轨方向误差分布特性，针对上述误差分布特性，采用不同权重，用于后续垂线偏差的精确计算，减轻误差噪声大的远幅侧和近幅侧对整体计算结果的影响。生成的垂线偏差产品可用于重力异常、海洋深度等其他海洋参数反演应用。
如图1和2所示，该方法具体包括：
步骤1)对原始数据中的每个原始采样点进行重采样，得到由多个重采样点组成的重采样数据，同时计算每个重采样点的海面高度误差；如图2所示，所示原始数据的格式为nc格式文件；其中，海面高度就是从读取的nc文件里的原始数据中提取的；
具体地，将原始数据中的每个50m×50m空间分辨率的原始采样点进行重采样，得到由多个2km×2km空间分辨率的重采样点组成的重采样数据，包括经度、纬度、海面高度，并在重采样的过程中，计算每个重采样点的海面高度误差：
重采样点邻域窗口尺寸范围为1km×1km，在重采样点邻域窗口内计算邻域窗口尺寸范围内的所有原始采样点的海面高度的标准差，作为该邻域窗口内的第i个重采样点的海面高度误差σi：
σi＝std(SSHnear)
其中，σi表示第i个重采样点的海面高度误差，SSHnear表示该重采样点对应的邻域窗口尺寸范围内的所有原始采样点的海面高度，符号std表示计算标准差。
重复上述计算过程，计算每个重采样点的海面高度误差；将第i个重采样点的经度、纬度、海面高度、海面高度误差组成卫星高度计重采样后的第i个重采样数据；并将其另存文件，便于后续步骤处理。
干涉成像高度计的重采样点的海面高度误差主要有基线倾角误差、基线长度误差、随机相位误差和海况偏差等。这几种误差项对垂轨方向的测高数据具有不同的影响，在近刈幅端，即靠近卫星地面轨迹的一侧误差较大；在远刈幅端，即远离卫星地面轨迹的一侧误差也较大；但是在观测刈幅的中间位置，误差较小。
步骤2)对重采样数据中的每个重采样点的海面高度进行异常值剔除和低通滤波，得到每个重采样点的平滑海面高度；
根据卫星高度计重采样后的所有重采样数据中的每个重采样点的海面高度SSH的经度和纬度组成的位置信息，在已发布的平均海面高数据集中，插值出对应点的平均海面高MSS；重采样后的所有海面高度与对应位置的平均海面高度做差，得到二者的差值，并计算该差值的标准差，记为σall；
判断每个重采样点的海面高度是否为异常值：
对第i个重采样点的海面高度SSHi与对应位置的平均海面高度MSSi做差，得到第i个重采样点的差值，若该第i个重采样点的差值的绝对值大于3σall，则认为该重采样点的海面高度是异常值，并删除该重采样点的海面高度；
|SSHi-MSSi|＞3×σall
若该第i个重采样点的差值的绝对值小于或等于3σall，则认为该重采样点的海面高度是正常值，保留该重采样点的海面高度；
根据判断结果，得到剔除异常值后的重采样数据点；
重复上述过程，对重采样数据中的每个重采样点的海面高度进行异常值剔除(即图2中的粗差剔除)，得到由多个剔除异常值后的重采样数据点组成的剔除后的数据集；
采用高斯加权滤波器，对剔除异常值后的数据集中的重采样数据点中的海面高度进行低通滤波，降低高频噪声影响，得到滤波后的海面高度，记为多个重采样点的平滑海面高度。
采用高斯加权滤波器，对保留的重采样点进行低通滤波，得到滤波后的平滑海面高度，记为每个重采样点的平滑海面高度；参数设计为滤波尺寸30km，半功率点对应7.5km波长。其中，进行低通滤波是为了降低噪声，平滑海面高度，提高信噪比；其中，采用高斯加权滤波器，能很好的保留第i个重采样点的平滑海面高度的高频成分。
步骤3)如图3所示，分别沿卫星轨道飞行方向、垂直于卫星轨道方向和倾斜于卫星轨道方向，根据步骤1)得到的每个重采样点的海面高度误差和步骤2)得到的每个重采样点的平滑海面高度，计算的每个重采样点的平滑数据；其中，上述平滑数据包括：每个重采样点的海面高度的海面坡度、海面方位角和海面坡度误差；
具体地，针对沿卫星轨道飞行方向：
如图3所示，沿卫星轨道飞行方向计算坡度信息示意，假设步骤2)中得到的每个重采样点的平滑海面高度中对应的重采样点为数据点，步骤1)中得到的重采样点的海面高度误差为该步骤中的数据点的海面高度误差；
将沿卫星轨道飞行方向的两个数据点的平滑海面高度做差，将该差值除以两个相邻数据点之间的真实距离，得到两个相邻数据点之间的沿卫星轨道飞行方向的海面坡度SSSn，具体公式如下：
其中，Hi,Hj分别表示两个相邻的数据点中的第i个数据点和第j个数据点的沿卫星轨道飞行方向的平滑海面高度，即步骤2)中得到的平滑海面高度；sij表示两个相邻数据点间的真实距离；
在计算两个相邻数据点之间的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向的沿卫星轨道飞行方向的海面方位角α：
其中，
是两个相邻的数据点中的第i个数据点的纬度；λi是两个相邻的数据点中的第i个数据点的经度；
是两个相邻的数据点中的第j个数据点的纬度；λj是两个相邻的数据点中的第j个数据点的经度；
此外，计算沿卫星轨道飞行方向的海面坡度误差σSSS：
其中，
是两个相邻数据点中的第i个数据点的海面高度误差；
分别是两个相邻数据点中的第j个数据点的海面高度误差；即该两个值
和
是通过步骤1)计算得到的；通过这个操作，将海面高度的误差分布特性传递到坡度误差分布。
将沿卫星轨道飞行方向的海面坡度SSSn，坡度方向的沿卫星轨道飞行方向的海面方位角α和沿卫星轨道飞行方向的海面坡度误差σSSS作为沿卫星轨道飞行方向的每个重采样点的平滑数据；
图4为垂轨方向坡度计算示意图，传统星下点观测高度计在垂轨向只有一列数据，导致坡度只能在单一沿轨方向计算，使得垂线偏差计算结果精度不高。而干涉成像高度计宽刈幅的特性，使得可以在垂轨方向(垂直于卫星轨道方向)和倾斜轨道方向(即倾斜于卫星轨道方向)计算坡度，大大提高了计算精度。计算方法与沿轨向一致，将垂轨方向的两个海面高度数据点做差，差值除以两个数据点之间的真实距离即为两个数据点之间的海面坡度。垂轨向的方位角计算和坡度误差计算同样按照沿轨向的计算方法进行。
针对垂直于卫星轨道方向：
图4为垂轨方向坡度计算示意，如图4所示，计算两个相邻数据点之间的垂轨向的海面坡度SSSN，具体公式如下：
其中，HI,HJ分别表示两个相邻的数据点中的第I个数据点和第J个数据点的垂直于卫星轨道方向的平滑海面高度；sIJ表示两个相邻数据点间的垂直于卫星轨道方向的真实距离；
在计算两个相邻数据点之间的垂轨向的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向垂轨向的海面方位角αN：
其中，
是两个相邻的数据点中的第I个数据点的垂轨向的纬度；λI是两个相邻的数据点中的第I个数据点的垂轨向的经度；
是两个相邻的数据点中的第J个数据点的垂轨向的纬度；λJ是两个相邻的数据点中的第J个数据点的垂轨向的经度；
此外，计算垂轨向海面坡度误差σSSSN：
其中，
是两个相邻数据点中的第I个数据点的垂轨向的海面高度误差；
分别是两个相邻数据点中的第J个数据点的垂轨向的海面高度误差；即该两个值
和
是通过步骤1)计算得到的；通过这个操作，将海面高度的误差分布特性传递到坡度误差分布。
将垂轨向的海面坡度SSSN，坡度方向垂轨向的海面方位角αN和垂轨向海面坡度误差σSSSN作为沿垂直于卫星轨道方向的每个重采样点的平滑数据；
图5为倾斜轨道135°方向的坡度计算示意图，图6为倾斜轨道45°方向的坡度计算示意图。倾斜于卫星飞行轨迹方向的坡度计算，及相应的方位角和坡度误差计算，能够按照沿轨向计算方法进行。需要注意的是两个数据点之间的距离s略大于沿轨向和垂轨向的距离s。
针对倾斜于卫星轨道135°方向，即图2中的斜轨坡度计算：
图5为倾斜轨道135°方向的坡度计算示意图，如图5所示，计算两个相邻数据点之间的倾斜于卫星轨道135°方向的海面坡度SSSM：
其中，HP,HQ分别表示两个相邻的数据点中的第P个数据点和第Q个数据点的倾斜于卫星轨道135°方向的平滑海面高度；sPQ表示两个相邻数据点间的倾斜于卫星轨道135°方向的真实距离；
在计算两个相邻数据点之间的倾斜于卫星轨道135°方向的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向倾斜于卫星轨道135°方向的海面方位角αM：
其中，
是两个相邻的数据点中的第P个数据点的倾斜于卫星轨道135°方向的纬度；λP是两个相邻的数据点中的第P个数据点的倾斜于卫星轨道135°方向的经度；
是两个相邻的数据点中的第Q个数据点的倾斜于卫星轨道135°方向的纬度；λQ是两个相邻的数据点中的第Q个数据点的倾斜于卫星轨道135°方向的经度；
此外，计算倾斜于卫星轨道135°方向海面坡度误差σSSSM：
其中，
是两个相邻数据点中的第P个数据点的倾斜于卫星轨道135°方向的海面高度误差；
分别是两个相邻数据点中的第Q个数据点的倾斜于卫星轨道 135°方向的海面高度误差；即该两个值
和
是通过步骤1)计算得到的；
将倾斜于卫星轨道135°方向的海面坡度SSSM，倾斜于卫星轨道135°方向的海面方位角αM和倾斜于卫星轨道135°方向海面坡度误差σSSSM作为沿垂直于卫星轨道135°方向的每个重采样点的平滑数据；
针对倾斜于卫星轨道45°方向，即图2中的斜轨坡度计算：
图6为倾斜轨道45°方向的坡度计算示意图，如图6所示，计算两个相邻数据点之间的倾斜于卫星轨道45°方向的海面坡度SSSO：
其中，Hg,Hh分别表示两个相邻的数据点中的第g个数据点和第h个数据点的倾斜于卫星轨道45°方向的平滑海面高度；sgh表示两个相邻数据点间的倾斜于卫星轨道45°方向的真实距离；
在计算两个相邻数据点之间的倾斜于卫星轨道45°方向的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向倾斜于卫星轨道45°方向的海面方位角αO：
其中，
是两个相邻的数据点中的第g个数据点的倾斜于卫星轨道45°方向的纬度；λg是两个相邻的数据点中的第g个数据点的倾斜于卫星轨道45°方向的经度；
是两个相邻的数据点中的第h个数据点的倾斜于卫星轨道45°方向的纬度；λh是两个相邻的数据点中的第h个数据点的倾斜于卫星轨道45°方向的经度；
此外，计算倾斜于卫星轨道45°方向海面坡度误差σSSSO：
其中，
是两个相邻数据点中的第g个数据点的倾斜于卫星轨道45°方向的海面高度误差；
分别是两个相邻数据点中的第h个数据点的倾斜于卫星轨道45°方向的海面高度误差；即该两个值
和
是通过步骤1)计算得到的；
将倾斜于卫星轨道45°方向的海面坡度SSSO，倾斜于卫星轨道45°方向的海面方位角αO和倾斜于卫星轨道45°方向海面坡度误差σSSSO作为沿垂直于卫星轨道 45°方向的每个重采样点的平滑数据。
步骤4)将每个重采样点的平滑数据作为一个数据点，合并组成数据集；
具体地，将沿卫星轨道飞行方向的海面坡度SSSn，沿卫星轨道飞行方向的海面方位角α和沿卫星轨道飞行方向的海面坡度误差σSSS作为沿卫星轨道飞行方向的每个重采样点的平滑数据；
将垂轨向的海面坡度SSSN，坡度方向垂轨向的海面方位角αN和垂轨向海面坡度误差σSSSN作为沿垂直于卫星轨道方向的每个重采样点的平滑数据；
将倾斜于卫星轨道135°方向的海面坡度SSSM，倾斜于卫星轨道135°方向的海面方位角αM和倾斜于卫星轨道135°方向海面坡度误差σSSSM作为沿垂直于卫星轨道135°方向的每个重采样点的平滑数据；
将倾斜于卫星轨道45°方向的海面坡度SSSO，倾斜于卫星轨道45°方向的海面方位角αO和倾斜于卫星轨道45°方向海面坡度误差σSSSO作为沿垂直于卫星轨道 45°方向的每个重采样点的平滑数据；
将上述每个方向的平滑数据作为一个数据点，合并组成数据集。
步骤5)根据实际需求，确定均匀分布的经纬度网格的尺寸，对数据集进行空间划分，得到多个经纬度网格点；其中，每个经纬度网格点内包括多个数据点；
具体地，根据实际需求，确定均匀分布的经纬度网格的尺寸为2′×2′或1′×1′的经纬度网格，按照上述确定的经纬度网格尺寸，对数据集进行空间划分，得到多个经纬度网格点；其中，每个经纬度网格点内包括多个数据点。
步骤6)对每个经纬度网格点，采用小邻域窗口方法，挑选出小邻域窗口内的所有数据点，得到每个经纬度网格点的挑选的所有数据点；
具体地，挑选出小邻域窗口内的所有数据点，得到挑选的所有数据点；
对于某个经纬度网格的第一个数据点，选择小邻域窗口的尺寸为12′×12′，即选定经纬度范围为网格数据点的经纬度加减6′，挑选出处于小于小邻域窗口尺寸范围内的所有数据点，得到挑选的所有数据点；凡是数据点经纬度落在选定经纬度范围内的，就将该数据点挑出。这样从整体数据集中能够挑出具有不同方向的坡度数据，且数据经纬度都满足选定的经纬度范围。
对每个经纬度网格点，采用小邻域窗口方法，基于上述挑选过程，对步骤4)得到的数据集中的每个数据点进行挑选，将经纬度处于小于小邻域窗口尺寸范围内的数据点挑选出来，得到每个经纬度网格点的挑选的所有数据点。
步骤7)利用最小二乘最优法则，使用每个经纬度网格点挑选的所有数据点，计算每个经纬度网格点的垂线偏差分量。
具体地，利用最小二乘最优法则，使用每个经纬网格的挑选的所有数据点，计算每个经纬网格的垂线偏差分量；其中，每个经纬网格的垂线偏差分量包括：每个经纬网格的垂线偏差北向分量ξ和每个经纬网格的垂线偏差南向分量η；
SSSi+vi＝ξcosαi+ηsinαi,i＝1,2,…,n
其中，SSSi表示第i个数据点的海面坡度；vi表示第i个数据点的残差项，αi是第 i个数据点的海面坡度SSSi的方位角；
基于最小二乘准则和间接平差方法，直接计算该经纬网格处的垂线偏差北向分量ξ和东向分量η；
将上式改写成矩阵形式：
V＝AX-L
其中，
因此，V的表现形式改写为：
假定卫星对地面的观测点都是独立观测，即不同观测点之间不存在互相影响，那么权重P的形式为：
其中，si是两个数据点之间的距离；σSSSZ是第i个数据点的海面坡度误差；该坡度误差是通过步骤3)的四个方向的坡度误差整合得到的；其中，四个方向指的是沿卫星轨道飞行方向、垂轨方向、倾斜于卫星轨道45°方向和倾斜于卫星轨道135°方向。
最终的第m个经纬度网格点的垂线偏差北向分量ξ和东向分量η为：
其中，ξm为第m个经纬度网格点的垂线偏差北向分量；ηm为第m个经纬度网格点的垂线偏差东向分量；
使用第m个经纬度网格点的挑选出的所有数据点，计算得到该经纬度网格点处的垂线偏差分量；
重复上述过程，计算每个经纬度网格点处的垂线偏差分量。
本发明还提供了一种基于干涉成像高度计的垂线偏差计算系统，该系统包括：
高度误差获取模块，用于对原始数据中的每个原始采样点进行重采样，得到由多个重采样点组成的重采样数据，同时计算每个重采样点的海面高度误差；
海面高度获取模块，用于对重采样数据中的每个重采样点的海面高度进行异常值剔除和低通滤波，得到多个重采样点的平滑海面高度；
平滑数据获取模块，用于分别沿卫星轨道飞行方向、垂直于卫星轨道方向和倾斜于卫星轨道方向，根据得到的每个重采样点的海面高度误差和得到的每个重采样点的平滑海面高度，计算的每个重采样点的平滑数据；其中，上述平滑数据包括：每个重采样点的海面高度的坡度、方位角和坡度误差；将每个重采样点的平滑数据作为一个数据点，合并组成数据集；
经纬划分模块，用于根据实际需求，确定均匀分布的经纬度网格的尺寸，对数据集进行空间划分，得到多个经纬度网格点；其中，每个经纬度网格点内包括多个数据点；
挑选模块，用于对每个经纬度网格点，采用小邻域窗口方法，挑选出小邻域窗口内的所有数据点，得到每个经纬度网格点的挑选的所有数据点；和
垂线偏差计算模块，用于利用最小二乘最优法则，使用每个经纬度网格点的挑选的所有数据点，计算每个经纬网格的垂线偏差分量。
本发明还提供了一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现所述的方法。
本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行所述的方法。
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。
Claims (10)
1.一种基于干涉成像高度计的垂线偏差计算方法，该方法包括：
步骤1)对原始数据中的每个原始采样点进行重采样，得到由多个重采样点组成的重采样数据，同时计算每个重采样点的海面高度误差；
步骤2)对重采样数据中的每个重采样点的海面高度进行异常值剔除和低通滤波，得到多个重采样点的平滑海面高度；
步骤3)分别沿卫星轨道飞行方向、垂直于卫星轨道方向和倾斜于卫星轨道方向，根据步骤1)得到的每个重采样点的海面高度误差和步骤2)得到的每个重采样点的平滑海面高度，计算的每个重采样点的平滑数据；其中，上述平滑数据包括：每个重采样点的海面高度的海面坡度、海面方位角和海面坡度误差；
步骤4)将每个重采样点的平滑数据作为一个数据点，合并组成数据集；
步骤5)根据实际需求，确定均匀分布的经纬度网格的尺寸，对数据集进行空间划分，得到多个经纬度网格点；其中，每个经纬度网格点内包括多个数据点；
步骤6)对每个经纬度网格点，采用小邻域窗口方法，挑选出小邻域窗口内的所有数据点，得到每个经纬度网格点的挑选的所有数据点；
步骤7)利用最小二乘最优法则，使用每个经纬度网格点的挑选的所有数据点，计算每个经纬网格的垂线偏差分量。
2.根据权利要求1所述的基于干涉成像高度计的垂线偏差计算方法，其特征在于，所述步骤1)具体包括：
将原始数据中的每个50m×50m空间分辨率的原始采样点进行重采样，得到由多个2km×2km空间分辨率的重采样点组成的重采样数据，包括经度、纬度、海面高度，并在重采样的过程中，计算每个重采样点的海面高度误差：
重采样点邻域窗口尺寸范围为1km×1km，在重采样点邻域窗口内计算邻域窗口尺寸范围内的所有原始采样点的海面高度的标准差，作为该采样点的海面高度误差σSSH：
σSSHi＝std(SSHnear)
其中，σSSHi表示第i个重采样点的海面高度误差，SSHnear表示该重采样点对应的邻域窗口尺寸范围内的所有原始采样点的海面高度，符号std表示计算标准差；
将第i个重采样点的经度、纬度、海面高度、海面高度误差组成卫星高度计重采样后的第i个重采样数据。
3.根据权利要求1所述的基于干涉成像高度计的垂线偏差计算方法，其特征在于，所述步骤2)具体包括：
根据卫星高度计重采样后的所有重采样数据中的每个重采样点的海面高度SSH的经度和纬度组成的位置信息，在已发布的平均海面高数据集中，插值出对应点的平均海面高MSS；重采样后的所有海面高度与对应位置的平均海面高度做差，得到二者的差值，并计算该差值的标准差，记为σall；
判断每个重采样点的海面高度是否为异常值：
对第i个重采样点的海面高度SSHi与对应位置的平均海面高度MSSi做差，得到第i个重采样点的差值，若该第i个重采样点的差值的绝对值大于3σall，则认为该重采样点的海面高度是异常值，并删除该重采样点的海面高度；
|SSHi-MSSi|＞3×σall
若该第i个重采样点的差值的绝对值小于或等于3σall，则认为该重采样点的海面高度是正常值，保留该重采样点的海面高度；
根据判断结果，得到剔除异常值后的重采样数据点；
重复上述过程，对重采样数据中的每个重采样点的海面高度进行异常值剔除，得到由多个剔除异常值后的重采样数据点组成的剔除后的数据集；
采用高斯加权滤波器，对剔除异常值后的数据集中的重采样数据点中的海面高度进行低通滤波，降低高频噪声影响，得到滤波后的海面高度，记为多个重采样点的平滑海面高度。
4.根据权利要求1所述的基于干涉成像高度计的垂线偏差计算方法，其特征在于，所述步骤3)具体包括：
针对沿卫星轨道飞行方向：
假设步骤2)中得到的每个重采样点的平滑海面高度中对应的重采样点为数据点，将沿卫星轨道飞行方向的两个数据点的平滑海面高度做差，将该差值除以两个相邻数据点之间的真实距离，得到两个相邻数据点之间的沿卫星轨道飞行方向的海面坡度SSSn：
其中，Hi,Hj分别表示两个相邻的数据点中的第i个数据点和第j个数据点的沿卫星轨道飞行方向的平滑海面高度；sij表示两个相邻数据点间的真实距离；
在计算两个相邻数据点之间的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向的方位角α，记为沿卫星轨道飞行方向的海面方位角：
其中，
是两个相邻的数据点中的第i个数据点的纬度；λi是两个相邻的数据点中的第i个数据点的经度；
是两个相邻的数据点中的第j个数据点的纬度；λj是两个相邻的数据点中的第j个数据点的经度；
此外，计算沿卫星轨道飞行方向的海面坡度误差σSSS：
其中，
是两个相邻数据点中的第i个数据点的海面高度误差；
分别是两个相邻数据点中的第j个数据点的海面高度误差；
将沿卫星轨道飞行方向的海面坡度SSSn，沿卫星轨道飞行方向的海面方位角α和沿卫星轨道飞行方向的海面坡度误差σSSS作为沿卫星轨道飞行方向的每个重采样点的平滑数据；
针对垂直于卫星轨道方向：
计算两个相邻数据点之间的垂轨向的海面坡度SSSN：
其中，HI,HJ分别表示两个相邻的数据点中的第I个数据点和第J个数据点的垂直于卫星轨道方向的平滑海面高度；sIJ表示两个相邻数据点间的垂直于卫星轨道方向的真实距离；
在计算两个相邻数据点之间的垂轨向的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向垂轨向的海面方位角αN：
其中，
是两个相邻的数据点中的第I个数据点的垂轨向的纬度；λI是两个相邻的数据点中的第I个数据点的垂轨向的经度；
是两个相邻的数据点中的第J个数据点的垂轨向的纬度；λJ是两个相邻的数据点中的第J个数据点的垂轨向的经度；
此外，计算垂轨向海面坡度误差σSSSN：
其中，
是两个相邻数据点中的第I个数据点的垂轨向的海面高度误差；
分别是两个相邻数据点中的第J个数据点的垂轨向的海面高度误差；
将垂轨向的海面坡度SSSN，垂轨向的海面方位角αN和垂轨向海面坡度误差σSSSN作为沿垂直于卫星轨道方向的每个重采样点的平滑数据；
针对倾斜于卫星轨道135°方向：
计算两个相邻数据点之间的倾斜于卫星轨道135°方向的海面坡度SSSM：
其中，HP,HQ分别表示两个相邻的数据点中的第P个数据点和第Q个数据点的倾斜于卫星轨道135°方向的平滑海面高度；sPQ表示两个相邻数据点间的倾斜于卫星轨道135°方向的真实距离；
在计算两个相邻数据点之间的倾斜于卫星轨道135°方向的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向倾斜于卫星轨道135°方向的海面方位角αM：
其中，
是两个相邻的数据点中的第P个数据点的倾斜于卫星轨道135°方向的纬度；λP是两个相邻的数据点中的第P个数据点的倾斜于卫星轨道135°方向的经度；
是两个相邻的数据点中的第Q个数据点的倾斜于卫星轨道135°方向的纬度；λQ是两个相邻的数据点中的第Q个数据点的倾斜于卫星轨道135°方向的经度；
此外，计算倾斜于卫星轨道135°方向海面坡度误差σSSSM：
其中，
是两个相邻数据点中的第P个数据点的倾斜于卫星轨道135°方向的海面高度误差；
分别是两个相邻数据点中的第Q个数据点的倾斜于卫星轨道135°方向的海面高度误差；
将倾斜于卫星轨道135°方向的海面坡度SSSM，倾斜于卫星轨道135°方向的海面方位角αM和倾斜于卫星轨道135°方向海面坡度误差σSSSM作为沿垂直于卫星轨道135°方向的每个重采样点的平滑数据；
针对倾斜于卫星轨道45°方向：
计算两个相邻数据点之间的倾斜于卫星轨道45°方向的海面坡度SSSO：
其中，Hg,Hh分别表示两个相邻的数据点中的第g个数据点和第h个数据点的倾斜于卫星轨道45°方向的平滑海面高度；sgh表示两个相邻数据点间的倾斜于卫星轨道45°方向的真实距离；
在计算两个相邻数据点之间的倾斜于卫星轨道45°方向的海面坡度的同时，计算该两个相邻数据点之间的海面坡度的坡度方向倾斜于卫星轨道45°方向的海面方位角αO：
其中，
是两个相邻的数据点中的第g个数据点的倾斜于卫星轨道45°方向的纬度；λg是两个相邻的数据点中的第g个数据点的倾斜于卫星轨道45°方向的经度；
是两个相邻的数据点中的第h个数据点的倾斜于卫星轨道45°方向的纬度；λh是两个相邻的数据点中的第h个数据点的倾斜于卫星轨道45°方向的经度；
此外，计算倾斜于卫星轨道45°方向海面坡度误差σSSSO：
其中，
是两个相邻数据点中的第g个数据点的倾斜于卫星轨道45°方向的海面高度误差；
分别是两个相邻数据点中的第h个数据点的倾斜于卫星轨道45°方向的海面高度误差；
将倾斜于卫星轨道45°方向的海面坡度SSSO，倾斜于卫星轨道45°方向的海面方位角αO和倾斜于卫星轨道45°方向海面坡度误差σSSSO作为沿垂直于卫星轨道45°方向的每个重采样点的平滑数据。
5.根据权利要求1所述的基于干涉成像高度计的垂线偏差计算方法，其特征在于，所述步骤5)具体包括：
根据实际需求，确定均匀分布的经纬度网格的尺寸为2′×2′或1′×1′的经纬度网格，按照上述确定的经纬度网格尺寸，对数据集进行空间划分，得到多个经纬度网格点。
6.根据权利要求1所述的基于干涉成像高度计的垂线偏差计算方法，其特征在于，所述步骤6)具体包括：
挑选出小邻域窗口内的所有数据点，得到挑选的所有数据点；
对于某个经纬度网格点，选择小邻域窗口的尺寸为12′×12′，即选定经纬度范围为网格点的经纬度加减6′，对步骤4)得到的数据集中的每个数据点进行挑选，将经纬度处于小于小邻域窗口尺寸范围内的数据点挑选出来，得到每个经纬度网格的挑选的所有数据点。
7.根据权利要求1所述的基于干涉成像高度计的垂线偏差计算方法，其特征在于，所述步骤7)具体包括：
利用最小二乘最优法则，使用每个经纬网格的挑选的所有数据点，计算每个经纬网格的垂线偏差分量；其中，每个经纬网格的垂线偏差分量包括：每个经纬网格的垂线偏差北向分量ξ和每个经纬网格的垂线偏差南向分量η；
SSSi+vi＝ξcosαi+ηsinαi,i＝1,2,…,n
其中，SSSi表示第i个数据点的海面坡度；vi表示第i个数据点的残差项，αi是第i个数据点的海面坡度SSSi的方位角；
基于最小二乘准则和间接平差方法，直接计算该经纬网格处的垂线偏差北向分量ξ和东向分量η；
将上式改写成矩阵形式：
V＝AX-L
其中，
因此，V的表现形式改写为：
假定卫星对地面的观测点都是独立观测，即不同观测点之间不存在互相影响，那么权重P的形式为：
其中，si是两个数据点之间的距离；σSSSZ是第i个数据点的海面坡度误差；该坡度误差是通过步骤3)的四个方向的坡度误差整合得到的；其中，四个方向指的是沿卫星轨道飞行方向、垂轨方向、倾斜于卫星轨道45°方向和倾斜于卫星轨道135°方向；
最终的第m个经纬度网格点的垂线偏差北向分量ξ和东向分量η为：
其中，ξm为第m个经纬度网格点的垂线偏差北向分量；ηm为第m个经纬度网格点的垂线偏差东向分量；
使用第m个经纬度网格点的挑选出的所有数据点，计算得到该经纬度网格点处的垂线偏差分量；
重复上述过程，计算每个经纬度网格点处的垂线偏差分量。
8.一种基于干涉成像高度计的垂线偏差计算系统，其特征在于，该系统包括：
高度误差获取模块，用于对原始数据中的每个原始采样点进行重采样，得到由多个重采样点组成的重采样数据，同时计算每个重采样点的海面高度误差；
海面高度获取模块，用于对重采样数据中的每个重采样点的海面高度进行异常值剔除和低通滤波，得到多个重采样点的平滑海面高度；
平滑数据获取模块，用于分别沿卫星轨道飞行方向、垂直于卫星轨道方向和倾斜于卫星轨道方向，根据得到的每个重采样点的海面高度误差和得到的每个重采样点的平滑海面高度，计算的每个重采样点的平滑数据；其中，上述平滑数据包括：每个重采样点的海面高度的坡度、方位角和坡度误差；将每个重采样点的平滑数据作为一个数据点，合并组成数据集；
经纬划分模块，用于根据实际需求，确定均匀分布的经纬度网格的尺寸，对数据集进行空间划分，得到多个经纬度网格点；其中，每个经纬度网格点内包括多个数据点；
挑选模块，用于对每个经纬度网格点，采用小邻域窗口方法，挑选出小邻域窗口内的所有数据点，得到每个经纬度网格点的挑选的所有数据点；和
垂线偏差计算模块，用于利用最小二乘最优法则，使用每个经纬度网格点的挑选的所有数据点，计算每个经纬网格的垂线偏差分量。
9.一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1-7中任一项所述的方法。
10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行权利要求1-7中任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111153467.6A CN113900069A (zh) | 2021-09-29 | 2021-09-29 | 一种基于干涉成像高度计的垂线偏差计算方法及其系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111153467.6A CN113900069A (zh) | 2021-09-29 | 2021-09-29 | 一种基于干涉成像高度计的垂线偏差计算方法及其系统 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN113900069A true CN113900069A (zh) | 2022-01-07 |
Family
ID=79189390
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111153467.6A Pending CN113900069A (zh) | 2021-09-29 | 2021-09-29 | 一种基于干涉成像高度计的垂线偏差计算方法及其系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113900069A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113406628A (zh) * | 2021-05-12 | 2021-09-17 | 中国科学院国家空间科学中心 | 一种用于干涉成像高度计的数据重采样方法 |
| CN115629364A (zh) * | 2022-12-22 | 2023-01-20 | 中国海洋大学 | 一种面向动态海面的星载小角度sar海况偏差仿真方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20200059085A (ko) * | 2018-11-20 | 2020-05-28 | 서울대학교산학협력단 | 해양 표층 표류부이 관측 수온자료와 인공위성 적외영상자료를 활용한 해양 표층 해수면온도 산출 방법 |
| RU2727584C1 (ru) * | 2019-12-20 | 2020-07-22 | Федеральное государственное бюджетное образовательное учреждение высшего образования «Государственный университет морского и речного флота имени адмирала С.О. Макарова» | Альтиметрический способ выполнения морской высотометрии с определением уклонения отвесной линии и устройство для его осуществления |
| CN111505619A (zh) * | 2020-03-03 | 2020-08-07 | 自然资源部第一海洋研究所 | 一种时空分布不规则不均匀卫星高度计测高数据网格化处理方法 |
| CN112556660A (zh) * | 2021-02-20 | 2021-03-26 | 中国测绘科学研究院 | 基于卫星测高数据的海域重力异常反演方法及系统 |
-
2021
- 2021-09-29 CN CN202111153467.6A patent/CN113900069A/zh active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20200059085A (ko) * | 2018-11-20 | 2020-05-28 | 서울대학교산학협력단 | 해양 표층 표류부이 관측 수온자료와 인공위성 적외영상자료를 활용한 해양 표층 해수면온도 산출 방법 |
| RU2727584C1 (ru) * | 2019-12-20 | 2020-07-22 | Федеральное государственное бюджетное образовательное учреждение высшего образования «Государственный университет морского и речного флота имени адмирала С.О. Макарова» | Альтиметрический способ выполнения морской высотометрии с определением уклонения отвесной линии и устройство для его осуществления |
| CN111505619A (zh) * | 2020-03-03 | 2020-08-07 | 自然资源部第一海洋研究所 | 一种时空分布不规则不均匀卫星高度计测高数据网格化处理方法 |
| CN112556660A (zh) * | 2021-02-20 | 2021-03-26 | 中国测绘科学研究院 | 基于卫星测高数据的海域重力异常反演方法及系统 |
Non-Patent Citations (2)
| Title |
|---|
| 万剑华;王莉娟;范陈清;杨俊钢;: "联合多卫星测高数据确定中国近海及其邻域垂线偏差", 海洋学研究, no. 03, 15 September 2012 (2012-09-15) * |
| 孙杨;姜效典;王D;李朝阳;: "卫星测高反演重力场在西太平洋海域的应用", 中国海洋大学学报(自然科学版), no. 06, 15 June 2017 (2017-06-15) * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113406628A (zh) * | 2021-05-12 | 2021-09-17 | 中国科学院国家空间科学中心 | 一种用于干涉成像高度计的数据重采样方法 |
| CN115629364A (zh) * | 2022-12-22 | 2023-01-20 | 中国海洋大学 | 一种面向动态海面的星载小角度sar海况偏差仿真方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110058236B (zh) | 一种面向三维地表形变估计的InSAR和GNSS定权方法 | |
| Baltsavias et al. | Digital surface modelling by airborne laser scanning and digital photogrammetry for glacier monitoring | |
| Catalão et al. | Merging GPS and atmospherically corrected InSAR data to map 3-D terrain displacement velocity | |
| CN106970398B (zh) | 顾及卫星遮挡条件的卫星可见性分析及星历预报方法 | |
| EP3084350B1 (en) | Method and system for generating a geoid via three computation spaces and airborne-acquired gravity data | |
| CN108896040B (zh) | 天空海一体化水下潜器惯性/重力组合导航方法和系统 | |
| CN104597471A (zh) | 面向时钟同步多天线gnss接收机的定向测姿方法 | |
| CN113900069A (zh) | 一种基于干涉成像高度计的垂线偏差计算方法及其系统 | |
| CN109085655B (zh) | 一种水下平台重力测量方案与验证方法 | |
| CN113189559B (zh) | 一种星载成像高度计遥感数据海底地形反演方法 | |
| CN114689015B (zh) | 一种提高光学卫星立体影像dsm高程精度的方法 | |
| CN112146663B (zh) | 基于偏置星载加速度计的gnss卫星精密定轨增强方法及系统 | |
| CN104567802B (zh) | 集成船载重力和gnss的测线式陆海高程传递方法 | |
| CN111854699A (zh) | 一种基于无人机航测河道崩岸过程的监测方法 | |
| Zeybek et al. | Investigation of landslide detection using radial basis functions: a case study of the Taşkent landslide, Turkey | |
| Woodroffe et al. | Reference water level and tidal datum | |
| Fei et al. | Dem Development and Precision Analysis for Antarctic Ice Sheet Using Cryosat‐2 Altimetry Data | |
| CN115165019A (zh) | 一种在航gnss潮位测量方法 | |
| CN113238228B (zh) | 基于水准约束的三维地表形变获取方法、系统及装置 | |
| CN116861647A (zh) | 联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法 | |
| CN114046774B (zh) | 综合cors网和多源数据的地面形变连续监测方法 | |
| CN114252875B (zh) | 一种成像高度计数据的高精度网格化方法 | |
| Saray et al. | Evaluation of Global Gravitational Models Based on DGPS/leveling Data over Baghdad University (IRAQ) | |
| CN115018973A (zh) | 一种低空无人机点云建模精度的无靶标评估方法 | |
| Feng et al. | Automatic selection of permanent scatterers-based GCPs for refinement and reflattening in InSAR DEM generation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |