WO2021000361A1 - 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法 - Google Patents

一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法 Download PDF

Info

Publication number
WO2021000361A1
WO2021000361A1 PCT/CN2019/097374 CN2019097374W WO2021000361A1 WO 2021000361 A1 WO2021000361 A1 WO 2021000361A1 CN 2019097374 W CN2019097374 W CN 2019097374W WO 2021000361 A1 WO2021000361 A1 WO 2021000361A1
Authority
WO
WIPO (PCT)
Prior art keywords
data
time
matrix
missing
value
Prior art date
Application number
PCT/CN2019/097374
Other languages
English (en)
French (fr)
Inventor
张丰
陈奕君
刘仁义
叶伟文
Original Assignee
浙江大学
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by 浙江大学 filed Critical 浙江大学
Priority to JP2020556955A priority Critical patent/JP7004844B2/ja
Priority to US17/042,214 priority patent/US11790580B2/en
Publication of WO2021000361A1 publication Critical patent/WO2021000361A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/77—Retouching; Inpainting; Scratch removal
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/40—Filling a planar surface by adding surface attributes, e.g. colour or texture
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10024—Color image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20024—Filtering details
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation

Definitions

  • the invention belongs to the field of ocean remote sensing, and specifically relates to a method for reconstructing stationary ocean water color satellite data based on an empirical orthogonal function decomposition method.
  • Ocean remote sensing satellites can carry out large-scale and long-term observations of the global ocean, providing humans with a data source that cannot be replaced by other observation methods for understanding and understanding the ocean.
  • the emergence of geostationary orbiting ocean water color satellites provides a new perspective for ocean remote sensing research-monitoring ocean water color changes from time to time.
  • the world's first geostationary ocean color sensor-GOCI (Geostationary Ocean Color Imager, GOCI) provides for the first time ocean color data with a time resolution of up to an hour, making it possible to monitor changes in the ocean environment over time.
  • Missing data is a common problem in ocean remote sensing data.
  • data processing methods have been developed at home and abroad for data reconstruction, such as optimal interpolation, kriging interpolation, empirical orthogonal function Decomposition method, among which Empirical Orthogonal Function Decomposition (Data INterpolating Empirical Orthogonal Functions, DINEOF) is currently a more effective ocean remote sensing data reconstruction method, which has been applied to reconstruct a variety of ocean remote sensing data, such as sea surface temperature, chlorophyll a Concentration, total suspended solids concentration, sea surface salinity, etc. Due to its own limitations, the empirical orthogonal function decomposition method has many adaptability problems for the reconstruction of stationary ocean aqua satellite data.
  • the reconstruction research of Liu Xiaoming and Wang Menghua has a shorter time span than quarterly, and focuses on the data missing rate in summer and autumn. Relatively low time period.
  • the purpose of the present invention is to overcome the existing shortcomings and provide a method for reconstructing stationary ocean aqua satellite data based on the empirical orthogonal function decomposition method.
  • S2 Construct a two-dimensional data matrix based on the spatio-temporal field of the preprocessed ocean color remote sensing data set, perform empirical orthogonal function decomposition and reconstruction, and use Laplace smoothing filter to perform the time covariance matrix before each iteration of decomposition Filter; After reconstruction, get the optimal number of retained modes, and the corresponding temporal and spatial characteristic modes;
  • each step is preferably implemented in the following specific manner.
  • the preferred implementation manners of each step can be combined correspondingly without conflict, which does not constitute a limitation.
  • an abnormal pixel identification process is performed according to S11 to S15:
  • O i is the current pixel value of the concentration to be calculated
  • O m is the median of the sample data set O
  • MAD(O) is the median absolute deviation of the sample data set O;
  • w prox is the edge detection weight
  • w conc is the concentration detection weight
  • w prox + w conc 1;
  • the method for removing the extremely high missing rate data is:
  • the missing rate of remote sensing images of geostationary orbit ocean color satellites at any time is greater than the upper limit of the missing rate, the corresponding images at that time will be eliminated; in space, if any of all geostationary orbit ocean color satellite remote sensing images
  • the time series data missing rate of a location point is greater than the upper limit of the missing rate, and the pixel value of the location point in all images is eliminated.
  • the upper limit of the deletion rate is 95%.
  • step S2 is as follows:
  • S21 Set the spatiotemporal field of the preprocessed ocean water color remote sensing data set as an m ⁇ n two-dimensional data matrix X 0 , where m is the number of spatial location points in the geostationary orbit ocean water color satellite remote sensing image, and n is the ocean water color remote sensing The total number of moments represented by the images in the data set, and m>n, where m rows of matrix X 0 represent the time series values of m different spatial locations, and n columns represent the values of all spatial locations at n different moments;
  • Steps S231 to S234 are executed in each iteration:
  • u i is the i th column of the spatial eigenmode of the U P
  • v i is the i th column time of eigenmodes V P
  • X T X is the time covariance matrix of n ⁇ n;
  • C l before data representing Laplacian smoothing C l data representing the Laplacian smoothing filter
  • [alpha] is a smoothing intensity, 0 ⁇ ⁇ 1;
  • (u p ) i is the i-th element in the p-th column of the spatial eigenmode U P , Is the j-th element in the p- th column of the time characteristic mode V P , and ⁇ p is the corresponding singular value;
  • N is the total number of data in the cross-validation set; Are the reconstructed value and the true value of the p-th element in the cross-validation set;
  • root mean square error R If the root mean square error R is greater than the error threshold, update the reconstruction values of all missing data points in the matrix X obtained in S233 to the matrix X, and then repeat steps S231 to S233 with the updated matrix X until the root mean square When the error R converges below the threshold, stop repeating, and record the converged root mean square error R as the root mean square error R P under the current retained mode number P ;
  • the iteration range of the reserved modal number P is [1, k MAX ], and the maximum number of iterations k XAX ⁇ min(m, n).
  • step S23 if the number of times of repeating steps S231 to S233 for singular value decomposition and Laplace smoothing filtering reaches the maximum number of times, stop repeating, and record the final root mean square error R as the current reservation The root mean square error R P under the number of modes P.
  • step S3 is as follows:
  • v i (t) ⁇ v i1 (d1),v i2 (dk), ⁇ ,v iK (dK) ⁇
  • K is the total number of frequency intraday data Geostationary ocean color satellite remote sensing image;
  • dk represents the date included in the data in v ik (dk), and dk does not include the missing date corresponding to the extremely high missing rate data;
  • n is the total number of moments represented by the images in the ocean color remote sensing data set
  • c j (dk) is the j-th intrinsic modal function component
  • r n (dk) is the decomposition margin
  • the total daily data frequency K of the geostationary orbit ocean color satellite remote sensing image is 8.
  • step S4 is as follows:
  • Is a spatial characteristic mode The i-th element in the p-th column of Is the jth element in the pth column of the time characteristic mode V T , Is the corresponding singular value;
  • S42 Update the reconstructed values of all missing points in the matrix X obtained in S41 to the matrix X, and the remaining element values of the matrix X remain unchanged to obtain the spatiotemporal field of the ocean color remote sensing data set that complements the missing values, and then obtain A geostationary ocean aqua satellite data set covering all time and space.
  • the present invention has the following beneficial effects:
  • the present invention fully considers its characteristics of high time resolution, as well as the continuous high missing rate and complete missing of data missing.
  • concentration detection and edge detection are used to detect abnormal pixels in remote sensing data.
  • the obtained temporal eigenmodes are decomposed and interpolated twice to realize the reconstruction of high missing rate data and continuous missing data, and finally obtain a water color remote sensing reanalysis data set that maintains its original accuracy and full coverage of time and space.
  • the present invention has improved mode number retention (Mode), reconstruction accuracy (MSE, RMSE, MAE), signal interpretation rate and reconstruction rate, and the operation is simple and easy, and the reconstruction result is intuitive.
  • the present invention provides demonstration and reference for the reconstruction of other remote sensing data with high time resolution characteristics or water color satellite data with large time interval changes.
  • Figure 1 is a flowchart of a method for reconstructing stationary ocean aqua satellite data based on the empirical orthogonal function decomposition method
  • Figure 2 shows the comparison between the original and reconstructed images of GOCI total suspended solids concentration on June 2, 2017.
  • Figure 3 shows the results of GOCI total suspended solids concentration reconstruction accuracy on June 2, 2017.
  • FIG. 1 it is a flowchart of a method for reconstructing stationary ocean aqua satellite data based on the empirical orthogonal function decomposition method provided in a preferred embodiment of the present invention.
  • the main steps include 4 steps, S1 ⁇ S4:
  • S2 Construct a two-dimensional data matrix based on the spatio-temporal field of the preprocessed ocean color remote sensing data set, perform empirical orthogonal function decomposition and reconstruction, and use Laplace smoothing filter to perform the time covariance matrix before each iteration of decomposition Filter; After reconstruction, get the optimal number of retained modes, and the corresponding temporal and spatial characteristic modes.
  • S3 Decompose the time eigenmodes obtained in S2 into eigenmode vectors in columns, and then decompose each column of eigenmode vectors into multiple sub-time eigenmodes according to the time the data belongs to, and then divide each sub-time eigenmode Carry out empirical mode decomposition, use cubic spline interpolation to complement the intrinsic mode function components after decomposition and the data corresponding to the position of the extremely high missing rate data in the margin, and re-stack and merge to obtain the continuous and complete time characteristic mode of the time series .
  • S4 Use the spatial characteristic mode obtained in S2 and the complete temporal characteristic mode obtained in S3 to reconstruct a stationary ocean aqua satellite data set with full coverage in time and space.
  • step S1 is to process these abnormal pixels and prevent noise from being introduced in the reconstruction process.
  • concentration detection is to judge whether the data is "outlier" based on the average and standard deviation of all data.
  • the calculation of this method is relatively simple, but if there are a large number of outliers in the data set to be tested, it will greatly affect the description of the overall characteristics of the data by the mean and standard deviation, resulting in greater uncertainty in the detection of outliers.
  • the use of median and median absolute deviation (Median Absolute Deviation, MAD) to replace the mean and standard deviation is a more robust detection method.
  • the weighted results of concentration detection and edge detection are used as the standard in the preprocessing process, and abnormal pixel identification and processing are performed according to S11 to S15:
  • O i is the current pixel value of the concentration to be calculated
  • O m is the median of the sample data set O
  • MAD(O) is the median absolute deviation of the sample data set O.
  • pixels with O conc greater than 3 can be regarded as abnormal pixels.
  • w prox is the edge detection weight
  • w conc is the concentration detection weight
  • w prox + w conc 1 is satisfied.
  • w prox and w conc can be adjusted according to actual conditions.
  • weighted results to detect abnormal pixels on the one hand, it improves the stability of detection and avoids excessive punishment by a single method, which leads to more data missing and loss of original change characteristics of the data. On the other hand, it can also target the given original data.
  • the source and characteristics of abnormal pixels adjust the weight to improve the effect of abnormal value detection.
  • the detection threshold O threshold determines the number of abnormal pixels to be eliminated.
  • the detection threshold O threshold is set to 3.
  • the empirical orthogonal function decomposition method (hereinafter referred to as DINEOF) cannot reconstruct these extremely high missing rate data, and it needs to be eliminated before reconstruction, otherwise the reconstruction result of the high missing rate image is the average of the original remote sensing data set Value substitution.
  • the method for removing extremely high missing rate data is: in time, if the missing rate of the geostationary orbit ocean color satellite remote sensing image at any time is greater than the missing rate upper limit, then the image corresponding to that time is removed; In space, if the time series data missing rate of any point in all geostationary orbit ocean color satellite remote sensing images is greater than the upper limit of the missing rate, then the pixel value of that point in all images will be removed, and the removed pixel will not As data missing points.
  • the upper limit of the missing rate can generally be set to 95%.
  • the original data reconstructed by the DINEOF method should be continuous in time, and the data interval should be the same.
  • the data itself is continuous within the same day and the time interval is consistent, but after the high-frequency daytime data, there is a nighttime gap without data for a long time.
  • the nighttime gap usually exceeds 12 hours, and it is more frequent in winter. long.
  • the cloud cover time is usually longer than the monitoring period of the stationary ocean aqua satellite, there are many cases of continuous high or even complete lack of regional data, which further reduces the continuity of the data.
  • the present invention considers a filter to smooth the original remote sensing data, which is different from general smoothing methods (such as median filter, mean filter, etc.).
  • the filter eliminates sudden changes and noise while retaining the original change characteristics of the data.
  • the present invention introduces the Laplacian smoothing method, constructs a filter suitable for time series data smoothing, and applies it to the reconstruction process of DINEOF. The following first introduces the theory of the filtering process.
  • Laplacian filtering to m time series of the data matrix X(m ⁇ n) can be understood as extending the Laplacian smoothing filter (n ⁇ 1) acting on one-dimensional time series to two-dimensional, constructing Laplacian filter L is an n ⁇ n matrix, and the matrix X is filtered in the time domain:
  • v i is the ith column of the time characteristic mode V
  • the filter L can be combined with Combined with the SVD decomposition, the original matrix filtering is changed to the time covariance matrix filtering.
  • Laplace smoothing filter Carrying out eigenvector decomposition expression, we can know the time covariance matrix after Laplace smoothing filter:
  • the Laplacian filter L filters the columns of the time covariance matrix A first, and then filters the rows of the time covariance matrix A.
  • the root mean square error of DINEOF decomposition can determine the optimal value of ⁇ and P. Generally, the value of ⁇ is kept constant and the value of P is optimized.
  • step S2 the empirical orthogonal function decomposition and reconstruction process in this embodiment is mainly implemented through step S2, and the specific method is described in detail below:
  • the m rows of the matrix X 0 respectively represent the time series values of m different spatial location points
  • the n columns respectively represent the values of all spatial location points at n different moments.
  • Steps S231 to S234 are executed in each iteration:
  • u i is the i th column of the spatial eigenmode of the U P
  • v i is the i th column time of eigenmodes V P
  • X T X is the time covariance matrix of n ⁇ n
  • the superscript T represents transpose (the same below).
  • C l before data representing Laplacian smoothing C l data representing the Laplacian smoothing filter
  • [alpha] is a smoothing intensity, 0 ⁇ ⁇ 1.
  • the superscript s represents the s-th Laplacian smoothing filter, and the Laplacian smoothing filter number s is 1,2,...,s max .
  • the smoothing effect of filtering depends on the value of ⁇ and the maximum number of iterations s max . This filtering realizes the smoothing of the time series data while retaining the change characteristics of the time series, and considers the influence of time relativity.
  • i, j are the space and time subscripts of X
  • (u p ) i is the i-th element in the p-th column of the spatial eigenmode U P
  • ⁇ p is the corresponding singular value.
  • the reconstruction values of all missing points of the data in the matrix X can be obtained. Since part of the data is taken from the matrix X in the foregoing steps as the cross-validation set, it can be based on the true value of the cross-validation set and The predicted reconstruction value is used to calculate the error.
  • N is the total number of data in the cross-validation set; Are the reconstructed value and the true value of the p-th element in the cross-validation set;
  • the updated matrix X re can be expressed as:
  • I a matrix of reconstructed values of missing points.
  • steps S231 to S233 are repeated.
  • Each repetition needs to calculate the root mean square error R, until the root mean square error R converges below the threshold, stop repeating, and record the converged root mean square error R as the root mean square error under the current retained mode number P R P.
  • an upper limit of the number of repetitions needs to be set, that is, the aforementioned s max . If repeating steps S231 ⁇ S233 for singular value decomposition and Laplace smoothing and filtering reaches the maximum number of times s max , do not continue to repeat, directly record the final root mean square error R as the current retained mode number P Root mean square error R P.
  • the optimal number of retained modes can be determined Obtain the optimal singular value matrix Time eigenmode And optimal spatial eigenmode
  • DINEOF relies on the statistical information of the original data for reconstruction. For extremely high missing rate data with a missing rate greater than 95%, it cannot be accurately reconstructed. Therefore, it needs to be eliminated before reconstruction. This part of the eliminated data does not participate in the S2 step.
  • the Laplacian filtering process needs to rely on step S3 for interpolation and reconstruction.
  • the time characteristic modal reconstructed by the DINEOF method itself reflects the time series change characteristics of the data set, and the data of adjacent times have the correlation and continuity of the time characteristic modal value corresponding to the same modal. Therefore, it is a feasible method to interpolate the time eigenmodes to obtain the time eigenmode vector corresponding to the extremely high missing rate. But for stationary ocean aqua satellite data, this method has certain flaws and cannot be directly used to reconstruct data with extremely high missing rates. The main reason is that the continuity of the temporal eigenmodes obtained from the reconstruction of the stationary ocean aqua satellite data is limited to the day on the time resolution scale.
  • the use of neighboring data interpolation on different days needs to correspond the temporal characteristic modal to the time resolution, and perform non-equidistant interpolation.
  • the traditional method has poor interpolation effect.
  • the data at all times during the day will be completely missing data continuously.
  • the interpolation results are difficult to match the actual situation.
  • the present invention adopts a method of reducing time resolution, using high-time-resolved remote sensing data that still has good periodicity at low time-resolution, and proposes to first reduce the time characteristic modal obtained by reconstruction.
  • Time resolution decomposition construct multiple sub-time eigenmodes with high continuity and consistent time intervals, and then use empirical mode decomposition method to decompose the sub-time eigenmodes, and use cubic spline method for missing signals in the decomposition results Interpolation fills, and finally merges into sub-time eigenmodes with continuous time and consistent intervals, and then obtains a new time eigenmode that includes all times.
  • This method decomposes the time eigenmodes twice and then interpolates.
  • the invention calls it "time mode quadratic decomposition interpolation".
  • the core of this method is that the data acquisition time using stationary ocean aqua satellite data is at equal time intervals at low time resolution, that is, the data acquired at the same time for multiple consecutive days have the same time interval (24 hours). There is periodicity.
  • the choice to interpolate the intrinsic modal function components obtained through empirical mode decomposition instead of directly interpolating the sub-time characteristic mode vector is because the decomposed signal volatility is lower than the original signal, which is beneficial to Interpolation and filling, and the intrinsic modal function component contains more details that the original vector does not have, which can achieve higher accuracy.
  • data obtained at the same time on different dates are referred to as data at the same time.
  • the data interval of the same time on two adjacent days is a fixed 24h.
  • step S3 the specific method for performing "time mode quadratic decomposition interpolation" in step S3 is as follows:
  • v i (t) ⁇ v i1 (d1),v i2 (dk), ⁇ ,v iK (dK) ⁇
  • K is the total daily data frequency of the geostationary ocean water color satellite remote sensing image.
  • the GOCI data used in the present invention scans the area 8 times a day from 8:30 to 15:30 Beijing time, and obtains 8 images with a time interval of 1 hour. , So the value of K is 8.
  • v ik (dk) of the v i (t) of the k-th time times corresponding to the sub-time eigenmodes which comprises a k-th time views of data for all dates of v i (t);
  • dk represents v ik (dk The date included in the data in ), and because the extremely high missing rate data in the data is eliminated, the date corresponding to these data in the decomposed sub-time characteristic mode is missing, so it is represented by dk, that is, dk does not include the high rate of missing data corresponding to the missing date and v i t (t) is also used to denote missing data. The next steps will perform interpolation completion on the data at these missing date positions.
  • n is the total number of moments represented by the images in the ocean color remote sensing data set
  • c j (dk) is the j-th intrinsic modal function component
  • r n (dk) is the decomposition margin
  • step S4 a complete spatial characteristic mode and a complete temporal characteristic mode have been obtained, and a stationary ocean aqua satellite data set with full coverage in time and space can be reconstructed according to the DINEOF method.
  • the superscript T all means transposition, which is different from "the time T corresponding to the extremely high missing rate data";
  • S42 Update the reconstruction values of all missing data points in the matrix X obtained in S41 to the matrix X.
  • the updated matrix of is the spatio-temporal field of the ocean water color remote sensing data set with missing values.
  • the updated matrix can be re-decomposed and restored to obtain a geostationary ocean aqua satellite data set covering time and space. All geostationary ocean aqua satellite remote sensing images in this dataset are continuous and complete in time and space.
  • the missing rate (time of missing) of the time series data of valid locations in the study area is the highest 99.93% (2918) and the lowest 53.94% (1575).
  • the spatial missing rate of 2920 GOCI TSM data is the highest 100%, and the lowest is 0.89%.
  • a 9 ⁇ 9 spatial window is used for concentration detection.
  • the degree detection weight w conc and the edge detection weight w prox are both Set the detection threshold O threshold to 3, and remove the original data from abnormal pixels as data missing points. Because of the limitation of DINEOF on the original data missing rate, all pixels where the original data corresponding to the time when the image missing rate is greater than 95% and the time series missing rate greater than 95% are further eliminated.
  • the data removal operation reduced the number of available GOCI TSM data from 2920 to 1794, a reduction of 38.56%.
  • step S4 using the spatial characteristic mode obtained in S2 and the complete temporal characteristic mode obtained in S3, reconstruct a stationary ocean aqua satellite data set with full coverage in time and space.
  • the reconstruction method of the present invention is named DINEOF-G, which is compared with the original classic DINEOF method ( BARTH A,RIXEN M,et al.Reconstruction of incomplete oceanographic data sets using empirical orthogonal functions:application to the Adriatic Sea surface temperature[J].Ocean Modelling.2005,9(4):325-346.)
  • the number (Mode), reconstruction accuracy (MSE, RMSE, MAE), signal interpretation rate and reconstruction rate have been improved, as shown in Table 1 below.
  • the reconstruction accuracy is increased by 8% (calculated by MSE), and the data reconstruction rate is increased from 58% to 94%, which is an increase of 36%.
  • test two are consistent with the classic DINEOF
  • the parameters of test five are consistent with the DINEOF-G proposed by the present invention.
  • the sequence of the test parameters is consistent with the sequence of the operating procedures in the test. From the analysis of specific test parameters and test results, it can be seen that by comparing test two and test three, abnormal pixel detection can improve the reconstruction accuracy. Although some data is lost, it has less impact on the reconstruction rate. The impact of data on the rate of spatio-temporal missing is limited. Comparing Experiment 3 and Experiment 4, Laplacian smoothing filter further improves the reconstruction accuracy while significantly increasing the number of retained modes and correspondingly improving the signal interpretation rate.
  • test five uses the time mode quadratic decomposition and interpolation to reconstruct the data when the spatial missing rate is greater than 95% based on the reconstruction result of test four.
  • the reconstruction accuracy of Experiment 5 is slightly lower than that of Experiment 4, mainly due to the introduction of additional sample observations (data with a missing rate higher than 95% and less than 100%), but it is still better than Experiment 1 using the same sample. This also reflects the high accuracy of the interpolation results of the time mode quadratic decomposition.

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法。其步骤如下：1）对原始海洋水色遥感区域利用集中度检测和边缘检测剔除异常像元；2）在经验正交函数分解法中使用拉普拉斯平滑滤波对时间协方差矩阵进行平滑过滤；3）对经验正交函数分解法获得时间特征模态值先分解为子时间特征模态，再采用经验模态分解，最后进行插值，获得保持其原始精度和时空间全覆盖的水色遥感再分析数据集。本方法的优点在于充分考虑静止海洋水色卫星数据高时间分辨率的特点，以及数据缺失存在连续高缺失率和完全缺失情况，适用于静止海洋水色卫星数据的重构，其准确性好、解释率高。对于静止海洋水色卫星数据的重构具有十分重要的实际应用价值。

Description

一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法 技术领域
本发明属于海洋遥感领域，具体涉及一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法。
背景技术
海洋遥感卫星能够对全球海洋开展大范围、长时期的观测，为人类了解海洋和认识海洋提供了其他观测方式都无法替代的数据源。静止轨道海洋水色卫星的出现为海洋遥感研究提供了全新的视角——逐时监测海洋水色变化。世界上第一个静止轨道海洋水色传感器——GOCI(Geostationary Ocean Color Imager，GOCI)首次提供了时间分辨率达小时级的海洋水色数据，使得对海洋环境的逐时变化监测成为可能。然而受到海洋上空云、雾和霾的影响，数据出现连续高缺失率甚至完全缺失的情况较多，对数据的时空连续性造成了严重的影响，使得反演出的海洋水色产品可用性大大降低，不能充分发挥其应有的价值。
数据缺失是海洋遥感数据中的常见问题，为了补全缺失的遥感数据，国内外已经发展了多种数据处理方法进行数据重构，如最优插值法、克里金插值法、经验正交函数分解法等，其中经验正交函数分解法(Data INterpolating Empirical Orthogonal Functions，DINEOF)是目前较为有效的海洋遥感数据重构方法，已经被应用于重构多种海洋遥感数据，如海表温度，叶绿素a浓度，总悬浮物浓度，海表盐度等。经验正交函数分解法由于自身的限制，对于重构静止海洋水色卫星数据存在诸多适应性问题，Liu Xiaoming和Wang Menghua的重构研究其时间跨度短于季度，且集中于夏、秋季数据缺失率相对较低的时段。
发明内容
本发明的目的是克服现有的不足，提供一种基于经验正交函数分解法的静止海洋水色卫 星数据重构方法。
为实现本发明目的，提供的技术方案如下：
一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法，它的步骤如下：
S1：对待重构的原始海洋水色遥感数据集中的每幅静止轨道海洋水色卫星遥感影像进行预处理将异常像元设为数据缺失点，同时剔除数据集中时间或空间上缺失率高于上限的极高缺失率数据；
S2：以预处理后的海洋水色遥感数据集的时空场构建二维数据矩阵，进行经验正交函数分解重构，且在每次迭代分解前使用拉普拉斯平滑滤波对时间协方差矩阵进行过滤；重构后，得到最优的保留模态数，以及对应的时间特征模态和空间特征模态；
S3：将所述的时间特征模态按列分解为特征模态向量，每一列特征模态向量再按照数据所属的时次分解成多个子时间特征模态，再将每个子时间特征模态进行经验模态分解，利用三次样条插值法补全分解后固有模态函数分量以及余量中的极高缺失率数据所对应位置的数据，重新叠加合并获得时间序列连续的完整时间特征模态；
S4：利用S2中得到的空间特征模态以及S3中得到的完整时间特征模态，重构获得时空间全覆盖的静止海洋水色卫星数据集。
基于上述技术方案，各步骤优选采用如下具体方式实现。其中各步骤的优选实现方式在没有冲突的情况下均可进行相应组合，不构成限制。
优选的，所述步骤S1中，针对每幅静止轨道海洋水色卫星遥感影像，按照S11～S15进行异常像元识别处理：
S11：针对待重构的原始海洋水色遥感数据集中的每幅静止轨道海洋水色卫星遥感影像，分别以影像中每个有值像元为中心设定一个边长为L的空间窗口；
S12：分别针对每个有值像元，以该像元对应的空间窗口内所有有值像元作为样本数据集O，计算该像元的集中度O conc：
Figure PCTCN2019097374-appb-000001
式中，O i为待计算集中度的当前像元值，O m为样本数据集O的中位数；标准差δ的计算 公式为：δ＝k*MAD(O)，k为比例因子常量，MAD(O)为样本数据集O的中位数绝对偏差；
S13：分别针对每个有值像元，计算该像元的边缘检测值O prox，若该像元周围存在无值像元则O prox取正数值，若该像元周围不存在无值像元则O prox取零值；
S14：分别针对每个有值像元，计算该像元的异常检测结果O final：
O final＝w proxO prox+w concO conc
式中，w prox为边缘检测权重，w conc为集中度检测权重，且满足w prox+w conc＝1；
S15：根据预设的异常像元检测门槛值O threshold，将O final大于门槛值O threshold的像元作为异常像元，并将异常像元设为数据缺失点。
优选的，所述步骤S1中，所述的极高缺失率数据的剔除方法为：
在时间上，若任一时刻的静止轨道海洋水色卫星遥感影像的缺失率大于缺失率上限，则将该时刻对应的影像进行剔除；在空间上，若所有静止轨道海洋水色卫星遥感影像中任一位置点的时间序列数据缺失率大于缺失率上限，则将该位置点在所有影像中的像元值进行剔除。
进一步优选的，所述的缺失率上限为95％。
优选的，所述步骤S2的具体方法如下：
S21：将经过预处理的海洋水色遥感数据集的时空场设为m×n的二维数据矩阵X 0，m为静止轨道海洋水色卫星遥感影像中的空间位置点个数，n为海洋水色遥感数据集中影像代表的时刻总数，且m>n，其中矩阵X 0的m行分别代表m个不同空间位置点的时间序列值，n列分别代表n个不同时刻所有空间位置点的值；
S22：计算矩阵X 0中所有有效数据的平均值
Figure PCTCN2019097374-appb-000002
令
Figure PCTCN2019097374-appb-000003
从矩阵X中随机抽取部分有效数据作为交叉验证集X c，将矩阵X中所有数据缺失点以及被抽取到交叉验证集X c的元素均赋值为0；
S23：对保留模态数P进行迭代，计算不同P值下的空间特征模态U P、时间特征模态V P，以及交叉验证集的重构值与原始值的均方根误差R；其中每一轮迭代均执行步骤S231～S234：
S231：对矩阵X进行奇异值分解，通过特征向量分解来计算P个最大奇异值和奇异向量，进而获得m×P的空间特征模态U P、P×P的奇异值矩阵S P和n×P的时间特征模态V P，计算 方程为：
Figure PCTCN2019097374-appb-000004
Figure PCTCN2019097374-appb-000005
式中，u i是空间特征模态U P的第i列，v i是时间特征模态V P的第i列，ρ i是与u i和v i对应的奇异值，i＝1,...,P；X TX为n×n的时间协方差矩阵；
S232：使用拉普拉斯平滑滤波先对所述时间协方差矩阵的列进行过滤，然后再对所述时间协方差矩阵的行进行过滤；对于时间协方差矩阵的行或列，均将其视为包含n个时刻数据的时间序列数据集合C＝{c 1,c 2,···,c n}，C对应的时间集合为T＝{t 1,t 2,···,t n}，c l为第l个时刻t l的时间序列数据，l＝1,2,…,n，且满足t 1<t 2<····<t n；拉普拉斯平滑滤波公式为：
Figure PCTCN2019097374-appb-000006
Figure PCTCN2019097374-appb-000007
t′ l＝t l-1+t l
式中，
Figure PCTCN2019097374-appb-000008
表示拉普拉斯平滑滤波前的数据c l，
Figure PCTCN2019097374-appb-000009
表示拉普拉斯平滑滤波后的数据c l；α为平滑强度，0<α<1；
S233：完成一次拉普拉斯平滑滤波后，计算矩阵X中所有数据缺失点的重构值，其中位于矩阵X第i行第j列的缺失点的重构值
Figure PCTCN2019097374-appb-000010
为：
Figure PCTCN2019097374-appb-000011
式中：(u p) i是空间特征模态U P的第p列中第i个元素，
Figure PCTCN2019097374-appb-000012
是时间特征模态V P的第p列中第j个元素，ρ p是相对应的奇异值；
S234：计算所述交叉验证集中所有元素的重构值与真值的均方根误差R：
Figure PCTCN2019097374-appb-000013
式中：N为所述交叉验证集中的数据总个数；
Figure PCTCN2019097374-appb-000014
分别为交叉验证集中第p个元素的 重构值、真值；
若均方根误差R大于误差阈值，则将S233中得到的矩阵X中所有数据缺失点的重构值更新至矩阵X中，然后以更新后的矩阵X重复步骤S231～S233，直至均方根误差R收敛至阈值以下时，停止重复，并将收敛的均方根误差R记为当前保留模态数P下的均方根误差R P；
S24：对保留模态数P＝1,2,…,k MAX下的均方根误差R P进行比较，以R P最小对应的P值作为最优的保留模态数
Figure PCTCN2019097374-appb-000015
获得最优的时间特征模态
Figure PCTCN2019097374-appb-000016
和最优的空间特征模态
Figure PCTCN2019097374-appb-000017
优选的，所述步骤S23中，保留模态数P的迭代范围为[1,k MAX]，最大迭代次数k XAX≤min(m,n)。
进一步优选的，所述步骤S23中，若重复步骤S231～S233进行奇异值分解和拉普拉斯平滑滤波的次数达到最大次数时，停止重复，并将最终的均方根误差R记为当前保留模态数P下的均方根误差R P。
优选的，所述步骤S3的具体方法如下：
S31：将最优的时间特征模态
Figure PCTCN2019097374-appb-000018
中的每一列特征模态向量进行分解，分解时以日为周期将不同日中属于同一时次的数据分解到一个子时间特征模态中，第i列特征模态向量v i(t)被分解为：
v i(t)＝{v i1(d1),v i2(dk),···,v iK(dK)}
式中：K为静止轨道海洋水色卫星遥感影像的日内数据频次总数；v ik(dk)为v i(t)中第k个时次对应的子时间特征模态，其中包含v i(t)中所有日期的第k个时次的数据；
Figure PCTCN2019097374-appb-000019
k＝1,2,…,K；dk表示v ik(dk)中的数据包括的日期，dk中不包括所述极高缺失率数据对应的缺失日期；
S32：对于每一列的特征模态向量中的每一个子时间特征模态v ik(dk)进行经验分解：
Figure PCTCN2019097374-appb-000020
式中：n为海洋水色遥感数据集中影像代表的时刻总数，c j(dk)是第j个固有模态函数分量，r n(dk)为分解的余量；
S33：分别针对每一个子时间特征模态v ik(dk)，利用三次样条插值法对所有固有模态函数分量c j(dk)以及余量信号r n(dk)中数据缺失日期位置进行插值填补，然后再重新叠加所有插值填补后的固有模态函数分量和余量信号，获得时间序列完整的子时间特征模态v ik(Dk)，Dk包含有第k个时次的所有连续日期；
S34：将同一列特征模态向量的K个子时间特征模态v ik(Dk)重新叠加合并，得到包含有所有时刻的完整特征模态向量v i(T)，T包含了极高缺失率数据对应的时刻；
S35：将所有的
Figure PCTCN2019097374-appb-000021
列完整特征模态向量v i(T)重新叠加合并，获得时间序列连续的完整时间特征模态V T。
进一步优选的，所述的静止轨道海洋水色卫星遥感影像的日内数据频次总数K为8。
优选的，所述步骤S4的具体方法如下：
S41：利用S2中得到的空间特征模态
Figure PCTCN2019097374-appb-000022
以及S3中得到的完整时间特征模态V T，重新计算矩阵X中所有数据缺失点的重构值，其中位于矩阵X第i行第j列的缺失点的重构值
Figure PCTCN2019097374-appb-000023
为：
Figure PCTCN2019097374-appb-000024
式中：
Figure PCTCN2019097374-appb-000025
是空间特征模态
Figure PCTCN2019097374-appb-000026
的第p列中第i个元素，
Figure PCTCN2019097374-appb-000027
是时间特征模态V T的第p列中第j个元素，
Figure PCTCN2019097374-appb-000028
是相对应的奇异值；
S42：将S41中得到的矩阵X中所有数据缺失点的重构值更新至矩阵X中，矩阵X的其余元素值不变，得到补全缺失值的海洋水色遥感数据集的时空场，进而获得时空间全覆盖的静止海洋水色卫星数据集。
本发明与现有技术相比具有的有益效果：
本发明针对DINEOF方法的限制，充分考虑了其具有高时间分辨率的特点，以及数据缺失存在连续高缺失率和完全缺失情况，先采用集中度检测和边缘检测对遥感数据中的异常像元进行检查，避免异常值对DINEOF方法中交叉验证的干扰，再引入拉普拉斯平滑滤波对DINEOF分解中的时间协方差矩阵进行过滤，减少因夜间间隔和数据连续缺失造成数据断层，然后对重构获得的时间特征模态二次分解后插值，实现高缺失率数据和连续缺失数据 的重构，最终获得保持其原始精度和时空间全覆盖的水色遥感再分析数据集。本发明在保留模态数(Mode)、重构精度(MSE、RMSE、MAE)、信号解释率和重构率方面均有提升，且操作简单易行，重构结果直观。本发明对其他具有高时间分辨率特点的遥感数据或时间间隔变化较大的水色卫星数据的重构提供示范和借鉴。
附图说明
图1为基于经验正交函数分解法的静止海洋水色卫星数据重构方法流程图；
图2为2017年6月2日的GOCI总悬浮物浓度原始与重构图像对比。
图3为2017年6月2日的GOCI总悬浮物浓度重构精度结果。
具体实施方式
下面结合附图和具体实施方式对本发明做进一步阐述和说明。
如图1所示，为本发明的一个较佳实施例中提供的一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法流程图。其主要步骤包括4步，分别为S1～S4：
S1：对待重构的原始海洋水色遥感数据集中的每幅静止轨道海洋水色卫星遥感影像(GOCI影像)进行预处理识别异常像元，并将异常像元设为数据缺失点，同时剔除数据集中时间或空间上缺失率高于上限的极高缺失率数据。
S2：以预处理后的海洋水色遥感数据集的时空场构建二维数据矩阵，进行经验正交函数分解重构，且在每次迭代分解前使用拉普拉斯平滑滤波对时间协方差矩阵进行过滤；重构后，得到最优的保留模态数，以及对应的时间特征模态和空间特征模态。
S3：将S2中得到的时间特征模态按列分解为特征模态向量，每一列特征模态向量再按照数据所属的时次分解成多个子时间特征模态，再将每个子时间特征模态进行经验模态分解，利用三次样条插值法补全分解后固有模态函数分量以及余量中的极高缺失率数据所对应位置的数据，重新叠加合并获得时间序列连续的完整时间特征模态。
S4：利用S2中得到的空间特征模态以及S3中得到的完整时间特征模态，重构获得时空间全覆盖的静止海洋水色卫星数据集。
下面对本实施例中S1～S4的具体实现方式以及其效果进行详细描述。
首先，GOCI数据虽然已经经过了数据质量控制，但仍可能存在合理值域范围内的异常值，因此步骤S1的目的之一是为了对这些异常像元进行处理，防止重构过程中引入噪声。最普遍的集中度检测方式是根据所有数据的平均值和标准差来判断数据是否“离群”。该方法计算较为简单，但如果待测数据集合中存在大量异常值，将极大地影响平均值和标准差对数据整体特征的描述，从而导致异常值检测存在较大的不确定性。为了避免上述情况，使用中位数和中位数绝对偏差(Median Absolute Deviation，MAD)替换平均值和标准差是一种更为稳健的检测方式。
在本实施例中，针对每幅静止轨道海洋水色卫星遥感影像，预处理过程中以集中度检测和边缘检测的加权结果为标准，按照S11～S15进行异常像元识别、处理：
S11：针对待重构的原始海洋水色遥感数据集中的每幅静止轨道海洋水色卫星遥感影像，分别以影像中每个有值像元为中心设定一个边长为L的空间窗口，L可以根据实际进行调整。
S12：分别针对每个有值像元，以该像元对应的空间窗口内所有有值像元作为样本数据集O，计算该像元的集中度O conc：
Figure PCTCN2019097374-appb-000029
式中，O i为待计算集中度的当前像元值，O m为样本数据集O的中位数；标准差δ的计算公式为：δ＝k*MAD(O)，k为比例因子常量，MAD(O)为样本数据集O的中位数绝对偏差。
对于正态分布数据，假设正常值落在中间50％的区域，则异常值落在两侧50％的区域里，可以换算得到k＝1.4826。根据3σ标准，O conc大于3的像元可以被认定为异常像元。
S13：为了增强异常像元监测的稳定性，本发明依据异常像元的产生原因，即异常像元在云层边缘发生概率较高，采用一种简单的边缘检测方法标记云层边缘位置的像元：分别针对每个有值像元，计算该像元的边缘检测值O prox，若该像元周围的8个像元存在无值像元则O prox取正数值M，若该像元周围的8个像元不存在无值像元则O prox取零值。本实施例中，为了与集中度检测的结果保持数值一致，则如果像元周围存在无值像元，标记为M＝3，不存在标记为0，即：
Figure PCTCN2019097374-appb-000030
S14：分别针对每个有值像元，计算该像元的异常检测结果O final：
O final＝w proxO prox+w concO conc
式中，w prox为边缘检测权重，w conc为集中度检测权重，且满足w prox+w conc＝1。w prox、w conc可以根据实际进行调整。
使用加权结果对异常像元进行检测，一方面提高了检测的稳定性，避免单一方法的过度惩罚，导致数据缺失较多，损失数据原有的变化特征，另一方面还可以针对给定原始数据异常像元的来源和特点调整权重提高异常值检测的效果。
S15：根据预设的异常像元检测门槛值O threshold，将O final大于门槛值O threshold的像元作为异常像元，并将异常像元设为数据缺失点。阈值O threshold的高低决定了被剔除的异常像元的数量，本实施例中设置检测阈值O threshold为3。
另外，如果原始数据在时间维或空间维缺失率过高(通常以超过95％计)，则这些数据不但不能提供有效的信息，反而可能影响最终结果。因此，经验正交函数分解法(以下简称DINEOF)无法重构这些极高缺失率的数据，需要在重构前进行剔除，否则被高缺失率图像的重构结果就是用原始遥感数据集的平均值替代。
本实施例中，极高缺失率数据的剔除方法为：在时间上，若任一时刻的静止轨道海洋水色卫星遥感影像的缺失率大于缺失率上限，则将该时刻对应的影像进行剔除；在空间上，若所有静止轨道海洋水色卫星遥感影像中任一位置点的时间序列数据缺失率大于缺失率上限，则将该位置点在所有影像中的像元值进行剔除，被剔除的像元不作为数据缺失点。缺失率上限一般可以设为95％。
应用DINEOF方法重构的原始数据在时间上应保持连续，且数据间隔时间相同。对于静止海洋水色卫星数据，其数据本身在当日内是连续且时间间隔一致的，但在高频日间数据之后，存在长时间没有数据的夜间间隙，夜间间隙通常超过12小时，且在冬季更长。另一方面，由于云层的覆盖时间通常长于静止海洋水色卫星的监测周期，数据出现区域数据连续高缺失率甚至完全缺失的情况较多，这进一步降低了数据的连续性。在连续性较差的数据中 会存在一些数据的非正常变化特征，但这些特征并不能代表数据的整个时空变异性，导致DINEOF方法对少量原始数据和非正常变化过度拟合，结果重建图像不一定反映该给定图像的应有的时空变化，并造成重构数据有明显的断层，出现非常不合理的连续时次图像。本发明考虑一种滤波来对原始遥感数据进行平滑，与一般的平滑方式(如中值滤波、均值滤波等)不同，该滤波在消除突变、噪声的同时需要保留数据原有的变化特征。对此，本发明引入拉普拉斯平滑的方法，构建一个适用于时间序列数据平滑的滤波，并应用于DINEOF的重构过程中。下面首先介绍该滤波过程的理论。
对于数据矩阵X(m×n)的m个时间序列应用拉普拉斯滤波，可理解为将作用于一维时间序列的拉普拉斯平滑滤波(n×1)扩展到二维，构建的拉普拉斯滤波L是n×n的矩阵，对矩阵X在时间域中进行过滤：
Figure PCTCN2019097374-appb-000031
其中，
Figure PCTCN2019097374-appb-000032
是经过拉普拉斯平滑滤波后的原始矩阵。理论上该滤波可直接对X进行过滤，但对原始矩阵的处理计算量过大，需要计算m次。考虑到
Figure PCTCN2019097374-appb-000033
在奇异值分解(SVD)计算中，只需要最重要得前几个时间特征模态和空间特征模态，通过特征向量分解来计算P个最大奇异值和奇异向量，方程为：
Figure PCTCN2019097374-appb-000034
其中，v i是时间特征模态V的第i列，ρ i是相应的奇异值，i＝1，...，P。因为
Figure PCTCN2019097374-appb-000035
是n×n的矩阵，将其定义为时间协方差矩阵A。可将该滤波L与
Figure PCTCN2019097374-appb-000036
的SVD分解相结合，改对原始矩阵过滤为对时间协方差矩阵过滤。对经拉普拉斯平滑滤波过滤的矩阵
Figure PCTCN2019097374-appb-000037
进行特征向量分解表达，可知经过拉普拉斯平滑滤波的时间协方差矩阵：
Figure PCTCN2019097374-appb-000038
根据上式可知，拉普拉斯滤波L先对时间协方差矩阵A的列进行过滤，然后再过滤时间协方差矩阵A的行。由此在大幅减少计算量的同时(从m次降为2n次)，也降低了对缺失数据的敏感性(时间协方差矩阵中没有缺失值)。通过DINEOF分解的均方根误差可以确定最优的α和P的值。一般α值保持恒定，优化P值即可。
基于上述理论，本实施例中的经验正交函数分解重构过程主要通过步骤S2来实现，下面对其具体方法进行详细描述：
S21：将经过预处理的海洋水色遥感数据集的时空场设为m×n的二维数据矩阵X 0，m为静止轨道海洋水色卫星遥感影像中的空间位置点个数，n为海洋水色遥感数据集中影像代表的时刻总数，且m>n，其中矩阵X 0的每一行分别代表由一个空间位置点在不同时刻的像素值构成的时间序列，每一列分别代表一个时刻所有空间位置点的值，也就是某一时刻的一幅GOCI影像上的像素值。因此，矩阵X 0的m行分别代表m个不同空间位置点的时间序列值，n列分别代表n个不同时刻所有空间位置点的值。记I为数据缺失点集合(包括原始缺失的以及在S1中被识别的异常像元)，矩阵中的缺失值用NaN表示，记P为保留模态数，P≤min(m,n)。
S22：计算矩阵X 0中所有有效数据的平均值
Figure PCTCN2019097374-appb-000039
令
Figure PCTCN2019097374-appb-000040
表示将X 0中所有的有值元素的值减去平均值
Figure PCTCN2019097374-appb-000041
从矩阵X中随机抽取部分有效数据作为交叉验证集X c，一般为有效数据总量的1％，作为判断最优的保留模态数的标准。视X中位于交叉验证集同一位置的点也为缺失点，赋值为NaN。将矩阵X中所有数据缺失点以及被抽取到交叉验证集X c的元素均赋值为0，即用0替换X中所有值为NaN的点。
S23：对保留模态数P进行迭代，计算不同P值下的空间特征模态U P、时间特征模态V P，以及交叉验证集的重构值与原始值的均方根误差R。P为正自然数，在[1,k MAX]的取值区间内，开始迭代之前初始化P＝1，然后没迭代一轮P＝P+1，最大迭代次数k MAX≤min(m,n)。
每一轮迭代均执行步骤S231～S234：
S231：对矩阵X进行奇异值分解，通过特征向量分解来计算P个最大奇异值和奇异向量，进而获得m×P的空间特征模态U P、P×P的奇异值矩阵S P和n×P的时间特征模态V P，计算方程为：
Figure PCTCN2019097374-appb-000042
Figure PCTCN2019097374-appb-000043
式中，u i是空间特征模态U P的第i列，v i是时间特征模态V P的第i列，ρ i是与u i和v i对应 的奇异值，i＝1,...,P；X TX为n×n的时间协方差矩阵；上标T代表转置(下同)。
由此，可以获得当前保留模态数P下的空间特征模态U P、时间特征模态V P以及奇异值矩阵S P，供后续调用。
S232：使用拉普拉斯平滑滤波先对所述时间协方差矩阵的列进行过滤，然后再对所述时间协方差矩阵的行进行过滤。在对时间协方差矩阵的行或列进行过滤的时候，均将其视为包含n个时刻数据的时间序列数据集合C＝{c 1,c 2,···,c n}，C对应的时间集合为T＝{t 1,t 2,···,t n}，c l为第l个时刻t l的时间序列数据，l＝1,2,…,n，且满足t 1<t 2<····<t n。拉普拉斯平滑滤波公式为：
Figure PCTCN2019097374-appb-000044
Figure PCTCN2019097374-appb-000045
t′ l＝t l-1+t l
式中，
Figure PCTCN2019097374-appb-000046
表示拉普拉斯平滑滤波前的数据c l，
Figure PCTCN2019097374-appb-000047
表示拉普拉斯平滑滤波后的数据c l；α为平滑强度，0<α<1。上标s表示第s次拉普拉斯平滑滤波，拉普拉斯平滑滤波次数s取值为1,2,…,s max。s初始值为1，每进行一次滤波s＝s+1。
滤波的平滑效果取决于α和最大迭代次数s max的值。该滤波在保留时间序列变化特征的同时实现了对时间序列数据的平滑，并考虑了时间相对性的影响。
S233：每完成一次拉普拉斯平滑滤波后，计算矩阵X中所有数据缺失点的重构值。如果矩阵X中位于矩阵X第i行第j列的元素为缺失点，则该缺失点的重构值
Figure PCTCN2019097374-appb-000048
为：
Figure PCTCN2019097374-appb-000049
式中：i，j为X的空间和时间下标，(u p) i是空间特征模态U P的第p列中第i个元素，
Figure PCTCN2019097374-appb-000050
是时间特征模态V P的第p列中第j个元素，ρ p是相对应的奇异值。
由此，可以获得矩阵X中的所有数据缺失点的重构值，由于在前述的步骤中从矩阵X中取出了部分数据作为交叉验证集，因此可以在滤波后基于交叉验证集的真值与预测的重构值来计算误差。
S234：计算所述交叉验证集中所有元素的重构值与真值的均方根误差R：
Figure PCTCN2019097374-appb-000051
式中：N为所述交叉验证集中的数据总个数；
Figure PCTCN2019097374-appb-000052
分别为交叉验证集中第p个元素的重构值、真值；
若均方根误差R大于误差阈值，则将S233中得到的矩阵X中所有数据缺失点的重构值更新至矩阵X中，在该更新过程中，仅对数据缺失点进行更新，而其他位置保持不变。因此，更新后的矩阵X re可以表示为：
Figure PCTCN2019097374-appb-000053
式中，
Figure PCTCN2019097374-appb-000054
为由缺失点的重构值构成的矩阵。
然后以更新后的X re重新作为矩阵X，重复步骤S231～S233。每次重复都需要计算均方根误差R，直至均方根误差R收敛至阈值以下时，停止重复，并将收敛的均方根误差R记为当前保留模态数P下的均方根误差R P。但考虑到计算效率，需要设置一个重复的次数上限，即前述的s max。若重复步骤S231～S233进行奇异值分解和拉普拉斯平滑滤波的次数达到最大次数s max时，不再继续重复，直接将最终的均方根误差R记为当前保留模态数P下的均方根误差R P。
S24：在上述步骤S23执行中，保留模态数P的迭代范围为[1,k MAX]，可以得到不同保留模态数P＝1,2,…,k MAX下的均方根误差R P，对R P进行比较，以R P最小对应的P值作为最优的保留模态数
Figure PCTCN2019097374-appb-000055
获得最优的时间特征模态
Figure PCTCN2019097374-appb-000056
和最优的空间特征模态
Figure PCTCN2019097374-appb-000057
若当P＝a时，均方根误差R最小，因此最优的保留模态数
Figure PCTCN2019097374-appb-000058
经过上述S2步骤，最终可以确定最优的保留模态数
Figure PCTCN2019097374-appb-000059
获得最优的奇异值矩阵
Figure PCTCN2019097374-appb-000060
时间特征模态
Figure PCTCN2019097374-appb-000061
和最优的空间特征模态
Figure PCTCN2019097374-appb-000062
如前所述，DINEOF依赖原始数据的统计信息进行重构，对于缺失率大于95％的极高缺失率数据无法准确重构，因此需要在重构前进行剔除，这部分剔除数据没有参与S2步骤的拉普拉斯滤波过程，需要依赖于步骤S3进行插值重构。
利用DINEOF方法重构获得的时间特征模态本身反映的就是数据集的时间序列变化特 征，且相邻时次的数据在相同模态对应的时间特征模态值上具有的相关性和连续性，因此对时间特征模态进行插值获得极高缺失率时次对应的时间特征模态向量，是一种可行的方法。但对于静止海洋水色卫星数据，该方法存在一定的缺陷，无法直接用于重构极高缺失率时次的数据。主要是因为，对于静止海洋水色卫星数据重构获得的时间特征模态，在其时间分辨率尺度上的连续性仅限于日内。也就是说，利用非同日的邻近数据插值，需要将时间特征模态与时间分辨率对应，进行非等距离插值，对于存在周期性的特征模态数据，使用传统方法插值效果较差。另一方面，如果当日受云雾覆盖，导致日内所有时次数据都将是连续完全缺失数据。这种情况下使用邻近数值按时间间隔插值，其插值结果与实际情况难以相符。针对这些问题，本发明采用一种降低时间分辨率的做法，利用高时间分辨遥感数据在低时间分辨率上仍具有良好周期性的特点，提出先对重构获得的时间特征模态按较低时间分辨率分解，构建多个连续性较高且时间间隔一致的子时间特征模态，然后采用经验模态分解方法再分解子时间特征模态，对分解结果中的缺失信号利用三次样条法插值填补，最后合并成时间连续、间隔一致的子时间特征模态，进而获得新的包含所有时次的时间特征模态，该方法对时间特征模态进行了两次分解后再插值，故本发明称其为“时间模态二次分解插值”。该方法的核心，一方面在于利用静止海洋水色卫星数据的数据获取时间在低时间分辨率下是等时间间隔的，即连续多日相同时刻获取的数据具有相同的时间间隔(24小时)，数据存在周期性。另一方面，选择对经过经验模态分解获得的各固有模态函数分量进行插值，而不是直接对子时间特征模态向量进行插值，是因为经过分解的信号波动性相比原始信号降低，利于插值填补，同时固有模态函数分量包含了更多原始向量没有的细节，可以达到更高的精确度。
本发明中，为了方便表述，将不同日期中相同时刻(例如不同日期的早上8点)获得的数据称为同一时次的数据。相邻两日的同一时次的数据间隔是固定的24h。
本实施例中，步骤S3进行“时间模态二次分解插值”的具体方法如下：
S31：将最优的时间特征模态
Figure PCTCN2019097374-appb-000063
中的每一列特征模态向量进行分解，分解时以日为周期将不同日中属于同一时次的数据分解到一个子时间特征模态中，第i列特征模态向量v i(t)被分解为：
v i(t)＝{v i1(d1),v i2(dk),···,v iK(dK)}
式中：
Figure PCTCN2019097374-appb-000064
k＝1,2,…,K。K为静止轨道海洋水色卫星遥感影像的日内数据频次总数，在本发明中采用的GOCI数据每天从北京时间8:30到15:30扫描该区域8次，获得时间间隔为1小时的8幅影像，因此K取值为8。v ik(dk)为v i(t)中第k个时次对应的子时间特征模态，其中包含v i(t)中所有日期的第k个时次的数据；dk表示v ik(dk)中的数据包括的日期，且由于数据中的极高缺失率数据被剔除，因此分解后的子时间特征模态中这些数据对应的日期是缺失的，因此用dk表示，即dk中不包括所述极高缺失率数据对应的缺失日期，而v i(t)中的t也同样用于表示数据缺失。后续步骤将对这些缺失日期位置的数据进行插值补全。
S32：对于每一列的特征模态向量中的每一个子时间特征模态v ik(dk)进行经验分解：
Figure PCTCN2019097374-appb-000065
式中：n为海洋水色遥感数据集中影像代表的时刻总数，c j(dk)是第j个固有模态函数分量，r n(dk)为分解的余量；
S33：分别针对每一个子时间特征模态v ik(dk)，利用三次样条插值法对所有固有模态函数分量c j(dk)以及余量信号r n(dk)中数据缺失日期位置进行插值填补，然后再逆向采用之前经验分解公式，重新叠加所有插值填补后的固有模态函数分量和余量信号，获得时间序列完整的子时间特征模态v ik(Dk)，此时的子时间特征模态中数据缺失日期位置已经被插值填补了相应数据，因此用Dk表示，即Dk包含有第k个时次的所有连续日期。
S34：将同一列特征模态向量的K个子时间特征模态v ik(Dk)重新叠加合并，得到包含有所有时刻的完整特征模态向量v i(T)，T表示包含了极高缺失率数据对应的时刻；
S35：将所有的
Figure PCTCN2019097374-appb-000066
列完整特征模态向量v i(T)重新按序叠加合并成矩阵形式，获得时间序列连续的完整时间特征模态V T。
由此，经过前述的S1～S3，已经获得了完整的空间特征模态以及完整的时间特征模态，既可以按照DINEOF方法重构获得时空间全覆盖的静止海洋水色卫星数据集。下面简述步骤S4的具体重构方法：
S41：利用S2中得到的空间特征模态
Figure PCTCN2019097374-appb-000067
以及S3中得到的完整时间特征模态V T，重新计算矩阵X中所有数据缺失点的重构值，其中位于矩阵X第i行第j列的缺失点的重构值
Figure PCTCN2019097374-appb-000068
为：
Figure PCTCN2019097374-appb-000069
式中：上标T均表示转置，与“包含了极高缺失率数据对应的时刻T”有所区别；
Figure PCTCN2019097374-appb-000070
是空间特征模态
Figure PCTCN2019097374-appb-000071
的第p列中第i个元素，
Figure PCTCN2019097374-appb-000072
是时间特征模态V T的第p列中第j个元素，
Figure PCTCN2019097374-appb-000073
是与
Figure PCTCN2019097374-appb-000074
和
Figure PCTCN2019097374-appb-000075
相对应的奇异值；
S42：将S41中得到的矩阵X中所有数据缺失点的重构值更新至矩阵X中，在该更新过程中，仅对数据缺失点进行更新，而矩阵X其他位置保持不变，由此得到的更新后矩阵即为补全缺失值的海洋水色遥感数据集的时空场。通过该更新后矩阵，即可重新分解还原，获得时空间全覆盖的静止海洋水色卫星数据集，该数据集中的所有静止轨道海洋水色卫星遥感影像在时空上均是连续完整的。
下面基于上述实施例方法，将其应用至具体的实例中对其效果进行展示。具体的过程如前所述，不再赘述，下面主要展示其具体参数设置和实现效果。
应用例
下面以杭州湾2017年的GOCI总悬浮物浓度数据为例，对本发明进行具体描述，其具体步骤如下：
1)采用韩国海洋卫星中心提供的GOCI Level-1B数据，利用紫外大气校正方法(UV-AC)进行大气校正，使用经验模型反演总悬浮物浓度数据。对GOCI数据通过TSM反演模型反演，获得杭州湾2017年的总悬浮物(Total Suspended Matter，TSM)反演结果产品(GOCI TSM)2920幅。每幅GOCI TSM数据产品的覆盖范围是120.55°E～122.35°E，29.85°N～30.95°N，像元尺寸是0.005°×0.005°。2920幅GOCI TSM数据的总体缺失率为63.5％，在空间上，研究区域有效位置点的时间序列数据的缺失率(缺失时次)最高99.93％(2918幅)，最低53.94％(1575幅)。在时间上，2920幅GOCI TSM数据中的空间缺失率最高100％，最低0.89％，缺失率大于95％的数据有1124幅，占比38.49％。
2)按照前述的步骤S1，对2017年的2920幅GOCI总悬浮物浓度数据，设L＝9，即采用9×9的空间窗口进行集中度检测9×9的空间窗口进行集中度检测，集中度检测权重w conc和边缘检测权重w prox均为
Figure PCTCN2019097374-appb-000076
设置检测门槛O threshold为3，将异常像元作为数据缺失点去除其中的原始数据。因为DINEOF对原始数据缺失率的限制，进一步剔除图像缺失率大于95％对应时次的原始数据和时间序列缺失率大于95％对应位置点的所有像素。数据剔除操作使得可利用的GOCI TSM数据数量从2920减少到1794，减少了38.56％。海洋像素从39274减少到37152，减少了5.4％，主要剔除的像素主要集中的在近岸水域。剔除后重构数据的总体缺失率为38.18％。在DINEOF内计算基础是基于原始数据的均值和协方差，TSM数据不呈现高斯分布，需要在DINEOF之前执行原始数据的转换。因此，本发明对原始数据进行以10为底的对数化处理。
3)按照前述的S2步骤对预处理后的海洋水色遥感数据集的时空场构建二维数据矩阵，进行经验正交函数分解重构，且在每次迭代分解前使用拉普拉斯平滑滤波对时间协方差矩阵进行过滤；重构后，得到最优的保留模态数，以及对应的时间特征模态和空间特征模态。本实施例中，α设置为0.01，而s max设置为30。最终，当P＝169时，均方根误差R最小，因此最优的保留模态数
Figure PCTCN2019097374-appb-000077
获得该模态数下的奇异值矩阵
Figure PCTCN2019097374-appb-000078
时间特征模态
Figure PCTCN2019097374-appb-000079
和最优的空间特征模态
Figure PCTCN2019097374-appb-000080
4)按照前述的S3步骤将时间特征模态
Figure PCTCN2019097374-appb-000081
按列分解为特征模态向量，每一列特征模态向量再按照数据所属的时次分解成8个子时间特征模态，再将每个子时间特征模态进行经验模态分解，利用三次样条插值法补全分解后固有模态函数分量以及余量中的极高缺失率数据所对应位置的数据，重新叠加合并获得时间序列连续的完整时间特征模态；
5)按照前述的S4步骤，利用S2中得到的空间特征模态以及S3中得到的完整时间特征模态，重构获得时空间全覆盖的静止海洋水色卫星数据集。
将本发明的重构方法命名为DINEOF-G，其相比原有的经典DINEOF方法(
Figure PCTCN2019097374-appb-000082
BARTH A,RIXEN M,et al.Reconstruction of incomplete oceanographic data sets using empirical orthogonal functions:application to the Adriatic Sea  surface temperature[J].Ocean Modelling.2005,9(4):325-346.)在保留模态数(Mode)、重构精度(MSE、RMSE、MAE)、信号解释率和重构率方面均有提升，如下表1所示。在重构精度上提高了8％(按MSE计算)，数据重构率从58％提高到94％，提高了36％。
表1
Figure PCTCN2019097374-appb-000083
为了进一步分析本发明提出的DINEOF-G方法中各步骤的效果及对重构结果的影响，采用调整步骤及实验参数的方式，又设计了不同做法的对比实验，具体方案和试验结果如表2所示，表中试验参数表示所执行的步骤(包含缺失率参数)。
表2
Figure PCTCN2019097374-appb-000084
其中试验二的参数与经典DINEOF一致，试验五的参数与本发明提出的DINEOF-G一致。试验参数顺序与试验中的操作流程顺序一致。从具体试验参数和试验结果分析可知：对比试验二和试验三，异常像元检测可以提高重构精度，虽然损失了部分数据，但对重构率影响较少，即通过异常像元检测对原始数据在时空缺失率上的影响有限。对比试验三和试验四，拉普拉斯平滑滤波再进一步提高重构精度的同时，明显提升了保留模态数，也相应提升了信号解释率。反映拉普拉斯平滑滤波能较好的平滑由夜间间隙和数据连续缺失导致的数据断层，恢复了原始数据的部分特征。对比试验四和试验五，试验五在试验四重构结果的基础上，利用时间模态二次分解插值重构了空间缺失率大于95％时次的数据。试验五与试验四相比重构精度稍有下降，主要是引入了额外的样本观测值(缺失率高于95％，低于100％的数据)，但仍优于使用同样样本的试验一，这也反映出时间模态二次分解插值结果的精度较高。
为了直观展示重构效果，在本发明的DINEOF-G重构方法下，获得了杭州湾2017年365天每天8时次的GOCI TSM数据重构结果。图2是2017年6月2日的TSM原始数据与重构数据(上八张为原始图像，下八张为重构图像)，该日8个原始数据的总体缺失率为76.68％，其中8:30和10:30数据缺失率超过95％。可以看出，利用时间模态二次分解插值，对非连续缺失数据进行重构时，虽然没有使用9:30数据对应的时间特征模态进行插值，但数据仍具有良好的连续性，动态特征明显，这是经典DINEOF方法使用平均值重构高缺失率数据无法达到的效果。对于15:30数据中出现的云迹尾流，重构数据也能较好的过滤处理。对原始数据有值像元值与重构数据相同区域的重构像元值进行精度评价，评价结果如图3所示，原始值与重构预测值集中分布在y＝x直线附近，RSM为0.0098，RMSE为0.0989，MAE为0.0645，表明本发明的方法的各精度指标较为理想，对于静止海洋水色卫星数据的重构具有十分重要的实际应用价值。
以上所述的实施例只是本发明的一种较佳的方案，然其并非用以限制本发明。有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型。因此凡采取等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。

Claims (10)

  1. 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法，其特征在于它的步骤如下：
    S1：对待重构的原始海洋水色遥感数据集中的每幅静止轨道海洋水色卫星遥感影像进行预处理，将异常像元设为数据缺失点，同时剔除数据集中时间或空间上缺失率高于上限的极高缺失率数据；
    S2：以预处理后的海洋水色遥感数据集的时空场构建二维数据矩阵，进行经验正交函数分解重构，且在每次迭代分解前使用拉普拉斯平滑滤波对时间协方差矩阵进行过滤；重构后，得到最优的保留模态数，以及对应的时间特征模态和空间特征模态；
    S3：将所述的时间特征模态按列分解为特征模态向量，每一列特征模态向量再按照数据所属的时次分解成多个子时间特征模态，再将每个子时间特征模态进行经验模态分解，利用三次样条插值法补全分解后固有模态函数分量以及余量中的极高缺失率数据所对应位置的数据，重新叠加合并获得时间序列连续的完整时间特征模态；
    S4：利用S2中得到的空间特征模态以及S3中得到的完整时间特征模态，重构获得时空间全覆盖的静止海洋水色卫星数据集。
  2. 根据权利要求1所述的重构方法，其特征在于：所述步骤S1中，针对每幅静止轨道海洋水色卫星遥感影像，按照S11～S15进行异常像元识别处理：
    S11：针对待重构的原始海洋水色遥感数据集中的每幅静止轨道海洋水色卫星遥感影像，分别以影像中每个有值像元为中心设定一个边长为L的空间窗口；
    S12：分别针对每个有值像元，以该像元对应的空间窗口内所有有值像元作为样本数据集O，计算该像元的集中度O conc：
    Figure PCTCN2019097374-appb-100001
    式中，O i为待计算集中度的当前像元值，O m为样本数据集O的中位数；标准差δ的计算公式为：δ＝k*MAD(O)，k为比例因子常量，MAD(O)为样本数据集O的中位数绝对偏差；
    S13：分别针对每个有值像元，计算该像元的边缘检测值O prox，若该像元周围存在无值像元则O prox取正数值，若该像元周围不存在无值像元则O prox取零值；
    S14：分别针对每个有值像元，计算该像元的异常检测结果O final：
    O final＝w proxO prox+w concO conc
    式中，w prox为边缘检测权重，w conc为集中度检测权重，且满足w prox+w conc＝1；
    S15：根据预设的异常像元检测门槛值O threshold，将O final大于门槛值O threshold的像元作为异常像元，并将异常像元设为数据缺失点。
  3. 根据权利要求1所述的重构方法，其特征在于：所述步骤S1中，所述的极高缺失率数据的剔除方法为：
    在时间上，若任一时刻的静止轨道海洋水色卫星遥感影像的缺失率大于缺失率上限，则将该时刻对应的影像进行剔除；在空间上，若所有静止轨道海洋水色卫星遥感影像中任一位置点的时间序列数据缺失率大于缺失率上限，则将该位置点在所有影像中的像元值进行剔除。
  4. 根据权利要求3所述的重构方法，其特征在于：所述的缺失率上限为95％。
  5. 根据权利要求1所述的重构方法，其特征在于：所述步骤S2的具体方法如下：
    S21：将经过预处理的海洋水色遥感数据集的时空场设为m×n的二维数据矩阵X 0，m为静止轨道海洋水色卫星遥感影像中的空间位置点个数，n为海洋水色遥感数据集中影像代表的时刻总数，且m>n，其中矩阵X 0的m行分别代表m个不同空间位置点的时间序列值，n列分别代表n个不同时刻所有空间位置点的值；
    S22：计算矩阵X 0中所有有效数据的平均值
    Figure PCTCN2019097374-appb-100002
    令
    Figure PCTCN2019097374-appb-100003
    从矩阵X中随机抽取部分有效数据作为交叉验证集X c，将矩阵X中所有数据缺失点以及被抽取到交叉验证集X c的元素均赋值为0；
    S23：对保留模态数P进行迭代，计算不同P值下的空间特征模态U P、时间特征模态V P，以及交叉验证集的重构值与原始值的均方根误差R；其中每一轮迭代均执行步骤S231～S234：
    S231：对矩阵X进行奇异值分解，通过特征向量分解来计算P个最大奇异值和奇异向量，进而获得m×P的空间特征模态U P、P×P的奇异值矩阵S P和n×P的时间特征模态V P，计算 方程为：
    Figure PCTCN2019097374-appb-100004
    Figure PCTCN2019097374-appb-100005
    式中，u i是空间特征模态U P的第i列，v i是时间特征模态V P的第i列，ρ i是与u i和v i对应的奇异值，i＝1,...,P；X TX为n×n的时间协方差矩阵；
    S232：使用拉普拉斯平滑滤波先对所述时间协方差矩阵的列进行过滤，然后再对所述时间协方差矩阵的行进行过滤；对于时间协方差矩阵的行或列，均将其视为包含n个时刻数据的时间序列数据集合C＝{c 1,c 2,…,c n}，C对应的时间集合为T＝{t 1,t 2,…,t n}，c l为第l个时刻t l的时间序列数据，l＝1,2,…,n，且满足t 1<t 2<····<t n；拉普拉斯平滑滤波公式为：
    Figure PCTCN2019097374-appb-100006
    Figure PCTCN2019097374-appb-100007
    t′ l＝t l-1+t l
    式中，
    Figure PCTCN2019097374-appb-100008
    表示拉普拉斯平滑滤波前的数据c l，
    Figure PCTCN2019097374-appb-100009
    表示拉普拉斯平滑滤波后的数据c l；α为平滑强度，0<α<1；
    S233：完成一次拉普拉斯平滑滤波后，计算矩阵X中所有数据缺失点的重构值，其中位于矩阵X第i行第j列的缺失点的重构值
    Figure PCTCN2019097374-appb-100010
    为：
    Figure PCTCN2019097374-appb-100011
    式中：(u p) i是空间特征模态U P的第p列中第i个元素，
    Figure PCTCN2019097374-appb-100012
    是时间特征模态V P的第p列中第j个元素，ρ p是相对应的奇异值；
    S234：计算所述交叉验证集中所有元素的重构值与真值的均方根误差R：
    Figure PCTCN2019097374-appb-100013
    式中：N为所述交叉验证集中的数据总个数；
    Figure PCTCN2019097374-appb-100014
    分别为交叉验证集中第p个元素的 重构值、真值；
    若均方根误差R大于误差阈值，则将S233中得到的矩阵X中所有数据缺失点的重构值更新至矩阵X中，然后以更新后的矩阵X重复步骤S231～S233，直至均方根误差R收敛至阈值以下时，停止重复，并将收敛的均方根误差R记为当前保留模态数P下的均方根误差R P；
    S24：对保留模态数P＝1,2,…,k MAX下的均方根误差R P进行比较，以R P最小对应的P值作为最优的保留模态数
    Figure PCTCN2019097374-appb-100015
    获得最优的时间特征模态
    Figure PCTCN2019097374-appb-100016
    和最优的空间特征模态
    Figure PCTCN2019097374-appb-100017
  6. 根据权利要求5所述的重构方法，其特征在于：所述步骤S23中，保留模态数P的迭代范围为[1,k MAX]，最大迭代次数k MAX≤min(m,n)。
  7. 根据权利要求5所述的重构方法，其特征在于：所述步骤S23中，若重复步骤S231～S233进行奇异值分解和拉普拉斯平滑滤波的次数达到最大次数时，停止重复，并将最终的均方根误差R记为当前保留模态数P下的均方根误差R P。
  8. 根据权利要求1所述的重构方法，其特征在于：所述步骤S3的具体方法如下：
    S31：将最优的时间特征模态
    Figure PCTCN2019097374-appb-100018
    中的每一列特征模态向量进行分解，分解时以日为周期将不同日中属于同一时次的数据分解到一个子时间特征模态中，第i列特征模态向量v i(t)被分解为：
    v i(t)＝{v i1(d1),v i2(dk),···,v iK(dK)}
    式中：K为静止轨道海洋水色卫星遥感影像的日内数据频次总数；v ik(dk)为v i(t)中第k个时次对应的子时间特征模态，其中包含v i(t)中所有日期的第k个时次的数据；
    Figure PCTCN2019097374-appb-100019
    k＝1,2,…,K；dk表示v ik(dk)中的数据包括的日期，dk中不包括所述极高缺失率数据对应的缺失日期；
    S32：对于每一列的特征模态向量中的每一个子时间特征模态v ik(dk)进行经验分解：
    Figure PCTCN2019097374-appb-100020
    式中：n为海洋水色遥感数据集中影像代表的时刻总数，c j(dk)是第j个固有模态函数分量，r n(dk)为分解的余量；
    S33：分别针对每一个子时间特征模态v ik(dk)，利用三次样条插值法对所有固有模态函数分量c j(dk)以及余量信号r n(dk)中数据缺失日期位置进行插值填补，然后再重新叠加所有插值填补后的固有模态函数分量和余量信号，获得时间序列完整的子时间特征模态v ik(Dk)，Dk包含有第k个时次的所有连续日期；
    S34：将同一列特征模态向量的K个子时间特征模态v ik(Dk)重新叠加合并，得到包含有所有时刻的完整特征模态向量v i(T)；
    S35：将所有的
    Figure PCTCN2019097374-appb-100021
    列完整特征模态向量v i(T)重新叠加合并，获得时间序列连续的完整时间特征模态V T。
  9. 根据权利要求8所述的重构方法，其特征在于：所述的静止轨道海洋水色卫星遥感影像的日内数据频次总数K为8。
  10. 根据权利要求1所述的重构方法，其特征在于：所述步骤S4的具体方法如下：
    S41：利用S2中得到的空间特征模态
    Figure PCTCN2019097374-appb-100022
    以及S3中得到的完整时间特征模态V T，重新计算矩阵X中所有数据缺失点的重构值，其中位于矩阵X第i行第j列的缺失点的重构值
    Figure PCTCN2019097374-appb-100023
    为：
    Figure PCTCN2019097374-appb-100024
    式中：
    Figure PCTCN2019097374-appb-100025
    是空间特征模态
    Figure PCTCN2019097374-appb-100026
    的第p列中第i个元素，
    Figure PCTCN2019097374-appb-100027
    是时间特征模态V T的第p列中第j个元素，
    Figure PCTCN2019097374-appb-100028
    是相对应的奇异值；
    S42：将S41中得到的矩阵X中所有数据缺失点的重构值更新至矩阵X中，矩阵X的其余元素值不变，得到补全缺失值的海洋水色遥感数据集的时空场，进而获得时空间全覆盖的静止海洋水色卫星数据集。
PCT/CN2019/097374 2019-07-04 2019-07-23 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法 WO2021000361A1 (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2020556955A JP7004844B2 (ja) 2019-07-04 2019-07-23 経験的直交関数分解法による静止海色衛星データの再構築方法
US17/042,214 US11790580B2 (en) 2019-07-04 2019-07-23 Method for reconstructing geostationary ocean color satellite data based on data interpolating empirical orthogonal functions

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201910600999.6 2019-07-04
CN201910600999 2019-07-04

Publications (1)

Publication Number Publication Date
WO2021000361A1 true WO2021000361A1 (zh) 2021-01-07

Family

ID=68255160

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2019/097374 WO2021000361A1 (zh) 2019-07-04 2019-07-23 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法

Country Status (4)

Country Link
US (1) US11790580B2 (zh)
JP (1) JP7004844B2 (zh)
CN (1) CN110378858B (zh)
WO (1) WO2021000361A1 (zh)

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113408644A (zh) * 2021-07-02 2021-09-17 南京信息工程大学 卫星数据重构方法、探测上层海洋对台风响应的方法
CN113554105A (zh) * 2021-07-28 2021-10-26 桂林电子科技大学 一种基于时空融合的物联网缺失数据补全方法
CN113627296A (zh) * 2021-08-02 2021-11-09 齐鲁空天信息研究院 叶绿素a浓度插值方法、装置、电子设备及存储介质
CN113835397A (zh) * 2021-06-30 2021-12-24 浙江大学 基于b样条曲线和路径积分的线性数控加工路径平滑方法
CN113935956A (zh) * 2021-09-23 2022-01-14 中国矿业大学(北京) 一种二向混合建模矿区土壤含水量数据缺失修复方法
CN114124361A (zh) * 2022-01-27 2022-03-01 广东工业大学 一种海洋感知数据的融合通信方法及系统
CN117112514A (zh) * 2023-10-23 2023-11-24 山东同利新材料有限公司 基于对氯甲基苯乙烯生产数据的记录存储方法
CN117274503A (zh) * 2023-11-17 2023-12-22 中国科学院国家空间科学中心 一种面向空间环境数据的多粒度动态可视化方法与系统
CN117745777A (zh) * 2024-01-05 2024-03-22 哈尔滨工业大学 一种基于时空配准的遥感检测地表密集异常元去除方法
CN117807380A (zh) * 2024-01-02 2024-04-02 中国科学院西北生态环境资源研究院 一种时间序列数据补全方法、装置、存储介质及电子设备
CN118396842A (zh) * 2024-06-26 2024-07-26 中国科学院空天信息创新研究院 时间序列遥感影像缺失区域重建方法、装置和电子设备

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110378858B (zh) * 2019-07-04 2021-04-09 浙江大学 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法
CN111275307B (zh) * 2020-01-16 2023-09-05 生态环境部华南环境科学研究所 一种水质自动在线站高频连续观测数据质量控制方法
CN111664937B (zh) * 2020-06-08 2021-04-30 西安电子科技大学 一种部分相干高斯光束辐照度波动时间特性的确定方法
CN113344805B (zh) * 2021-05-18 2023-06-13 南京信息工程大学 基于eof的自动气象站观测资料修复方法
CN114238847B (zh) * 2021-10-29 2023-02-10 中国人民解放军61540部队 一种基于海洋实测数据的表层准地转重构方法及系统
CN117272812B (zh) * 2023-09-26 2024-07-09 昆明理工大学 一种低纬小区域电离层模型构建方法
CN117093806B (zh) * 2023-10-16 2024-02-09 国家卫星气象中心(国家空间天气监测预警中心) 基于遥感的全空间覆盖海上大气co2柱浓度计算方法
CN117095134B (zh) * 2023-10-18 2023-12-22 中科星图深海科技有限公司 一种三维海洋环境数据插值处理方法
CN117668477B (zh) * 2024-01-31 2024-04-26 山东科技大学 一种海洋大数据智能轻量化处理方法及系统
CN117975296B (zh) * 2024-04-01 2024-06-11 自然资源部第一海洋研究所 一种异常海洋锋面的卫星遥感探测提取识别方法及系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103020478A (zh) * 2012-12-28 2013-04-03 杭州师范大学 一种海洋水色遥感产品真实性检验的方法
CN104217098A (zh) * 2014-07-14 2014-12-17 浙江大学 一种走航断面监测与卫星遥感结合的海域碳收支计算方法
US9697614B2 (en) * 2014-12-08 2017-07-04 Mitsubishi Electric Research Laboratories, Inc. Method for segmenting and tracking content in videos using low-dimensional subspaces and sparse vectors
CN110378858A (zh) * 2019-07-04 2019-10-25 浙江大学 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法

Family Cites Families (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4003869B2 (ja) * 2002-04-05 2007-11-07 独立行政法人科学技術振興機構 地球観測衛星データのノイズ除去処理方法、ノイズ除去処理プログラム、ノイズ除去処理プログラムを記録した記録媒体
US7379592B2 (en) * 2003-01-21 2008-05-27 United States Of America As Represented By The Secretary Of The Navy System and method for significant dust detection and enhancement of dust images over land and ocean
US6868361B2 (en) * 2003-07-29 2005-03-15 Council Of Scientific And Industrial Research Method of determining the volume scattering function of ocean waters in the backward direction using a satellite ocean color sensor
KR100654368B1 (ko) * 2005-12-26 2006-12-05 한국항공우주연구원 정지궤도 해양관측 카메라의 해색정보 최적화 검출방법
US7908103B2 (en) * 2007-05-21 2011-03-15 Nilanjan Senroy System and methods for determining masking signals for applying empirical mode decomposition (EMD) and for demodulating intrinsic mode functions obtained from application of EMD
KR101291675B1 (ko) * 2012-02-02 2013-08-01 주식회사 해양기술이앤지 스마트폰을 이용한 위성자료 분석 시스템
CN102798384B (zh) * 2012-07-03 2014-09-17 天津大学 一种基于压缩采样的海洋遥感图像水色水温监测方法
KR101381292B1 (ko) * 2012-12-28 2014-04-04 한국해양과학기술원 인공위성 시스템의 제어 방법 및 장치.
KR101381293B1 (ko) * 2012-12-31 2014-04-04 한국해양과학기술원 관측영상 기반의 인공위성 성능 평가 방법 및 장치.
CN103400022A (zh) * 2013-06-08 2013-11-20 杭州师范大学 一种海表面温度遥感数据集等纬度重构方法
KR101576167B1 (ko) * 2014-05-09 2015-12-09 한국해양과학기술원 해수면 흐름 모니터링 장치
CN104361602B (zh) * 2014-11-26 2017-12-08 中国科学院遥感与数字地球研究所 一种基于modis图像的水体颜色检测方法及装置
KR101731234B1 (ko) * 2016-02-04 2017-05-02 한국해양과학기술원 위성 또는 항공기에서 촬영된 영상의 슬롯간 편차 보정 방법
KR101717167B1 (ko) * 2016-05-24 2017-03-16 한국해양과학기술원 해양 프론트 기반의 해양 정보 분석 장치 및 방법
KR101880616B1 (ko) * 2017-08-03 2018-07-23 한국해양과학기술원 해상풍과 해무 위성정보를 이용한 해무 예측 방법
CN108776719A (zh) * 2018-04-27 2018-11-09 南京信息工程大学 基于经验模态分解的高能电子通量小时预报模型建立方法
CN108981957B (zh) 2018-05-31 2021-01-05 西北工业大学 基于自组织神经网络及经验正交函数的水下温度场重构方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103020478A (zh) * 2012-12-28 2013-04-03 杭州师范大学 一种海洋水色遥感产品真实性检验的方法
CN104217098A (zh) * 2014-07-14 2014-12-17 浙江大学 一种走航断面监测与卫星遥感结合的海域碳收支计算方法
US9697614B2 (en) * 2014-12-08 2017-07-04 Mitsubishi Electric Research Laboratories, Inc. Method for segmenting and tracking content in videos using low-dimensional subspaces and sparse vectors
CN110378858A (zh) * 2019-07-04 2019-10-25 浙江大学 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法

Cited By (19)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113835397A (zh) * 2021-06-30 2021-12-24 浙江大学 基于b样条曲线和路径积分的线性数控加工路径平滑方法
CN113835397B (zh) * 2021-06-30 2024-02-09 浙江大学 基于b样条曲线和路径积分的线性数控加工路径平滑方法
CN113408644A (zh) * 2021-07-02 2021-09-17 南京信息工程大学 卫星数据重构方法、探测上层海洋对台风响应的方法
CN113408644B (zh) * 2021-07-02 2023-07-14 南京信息工程大学 卫星数据重构方法、探测上层海洋对台风响应的方法
CN113554105B (zh) * 2021-07-28 2023-04-18 桂林电子科技大学 一种基于时空融合的物联网缺失数据补全方法
CN113554105A (zh) * 2021-07-28 2021-10-26 桂林电子科技大学 一种基于时空融合的物联网缺失数据补全方法
CN113627296A (zh) * 2021-08-02 2021-11-09 齐鲁空天信息研究院 叶绿素a浓度插值方法、装置、电子设备及存储介质
CN113627296B (zh) * 2021-08-02 2023-10-27 齐鲁空天信息研究院 叶绿素a浓度插值方法、装置、电子设备及存储介质
CN113935956A (zh) * 2021-09-23 2022-01-14 中国矿业大学(北京) 一种二向混合建模矿区土壤含水量数据缺失修复方法
CN113935956B (zh) * 2021-09-23 2022-03-25 中国矿业大学(北京) 一种二向混合建模矿区土壤含水量数据缺失修复方法
CN114124361B (zh) * 2022-01-27 2022-04-26 广东工业大学 一种海洋感知数据的融合通信方法及系统
CN114124361A (zh) * 2022-01-27 2022-03-01 广东工业大学 一种海洋感知数据的融合通信方法及系统
CN117112514A (zh) * 2023-10-23 2023-11-24 山东同利新材料有限公司 基于对氯甲基苯乙烯生产数据的记录存储方法
CN117112514B (zh) * 2023-10-23 2024-01-09 山东同利新材料有限公司 基于对氯甲基苯乙烯生产数据的记录存储方法
CN117274503A (zh) * 2023-11-17 2023-12-22 中国科学院国家空间科学中心 一种面向空间环境数据的多粒度动态可视化方法与系统
CN117274503B (zh) * 2023-11-17 2024-02-09 中国科学院国家空间科学中心 一种面向空间环境数据的多粒度动态可视化方法与系统
CN117807380A (zh) * 2024-01-02 2024-04-02 中国科学院西北生态环境资源研究院 一种时间序列数据补全方法、装置、存储介质及电子设备
CN117745777A (zh) * 2024-01-05 2024-03-22 哈尔滨工业大学 一种基于时空配准的遥感检测地表密集异常元去除方法
CN118396842A (zh) * 2024-06-26 2024-07-26 中国科学院空天信息创新研究院 时间序列遥感影像缺失区域重建方法、装置和电子设备

Also Published As

Publication number Publication date
JP2021533427A (ja) 2021-12-02
CN110378858A (zh) 2019-10-25
CN110378858B (zh) 2021-04-09
US20230154081A1 (en) 2023-05-18
JP7004844B2 (ja) 2022-01-21
US11790580B2 (en) 2023-10-17

Similar Documents

Publication Publication Date Title
WO2021000361A1 (zh) 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法
CN106897707B (zh) 基于多源中分的特征影像时间序列合成方法及装置
CN111932457A (zh) 遥感影像高时空融合处理算法及装置
CN108830796A (zh) 基于谱空结合和梯度域损失的高光谱图像超分辨重构方法
Wang et al. MCT-Net: Multi-hierarchical cross transformer for hyperspectral and multispectral image fusion
CN114565843B (zh) 一种时间序列遥感图像融合方法
CN109829872B (zh) 一种用于内陆水体遥感的多时相多源遥感影像融合方法
WO2024036739A1 (zh) 水库水储量反演方法和装置
CN112052589B (zh) 基于modis逐日积雪产品的积雪覆盖比例的估算方法和装置
CN109033543A (zh) 一种地表异质区植被覆盖度估算方法、装置及设备
Sun et al. Unsupervised 3D tensor subspace decomposition network for spatial-temporal-spectral fusion of hyperspectral and multispectral images
CN114820741B (zh) 一种高光谱影像全波段超分重建方法
CN113705340B (zh) 一种基于雷达遥感数据的深度学习变化检测方法
Chen et al. Fusion of Hyperspectral-Multispectral images joining Spatial-Spectral Dual-Dictionary and structured sparse Low-rank representation
CN109359264B (zh) 一种基于modis的叶绿素产品降尺度方法及装置
Fang et al. MIMO-SST: Multi-Input Multi-Output Spatial-Spectral Transformer for Hyperspectral and Multispectral Image Fusion
Sun et al. Coupled temporal variation information estimation and resolution enhancement for remote sensing spatial–temporal–spectral fusion
CN111273376B (zh) 降尺度的海表净辐射确定方法、系统、设备及存储介质
CN116704362A (zh) 结合时空特征的grace卫星空窗期数据重构方法、系统及设备
Fayad et al. Vision Transformers, a new approach for high-resolution and large-scale mapping of canopy heights
Liu et al. SAR image super-resolution based on TV-regularization using gradient profile prior
CN115222837A (zh) 真彩云图生成方法、装置、电子设备及存储介质
Wu et al. Super-resolution Reconstruction of Night-light Images Based on Improved SRCNN
Kremezi et al. Data fusion for increasing monitoring capabilities of Sentinel optical data in marine environment
Mu et al. Prediction of North Atlantic Oscillation Index Associated with the Sea Level Pressure Using DWT‐LSTM and DWT‐ConvLSTM Networks

Legal Events

Date Code Title Description
ENP Entry into the national phase

Ref document number: 2020556955

Country of ref document: JP

Kind code of ref document: A

121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19935865

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 19935865

Country of ref document: EP

Kind code of ref document: A1