简介:这份资源面向卫星遥感、海洋学与地球物理领域的科研人员及技术开发者,围绕风云三号E星(FY-3E)搭载的GNOS-II仪器获取的不均匀分布时延-多普勒(DDM)数据,系统复现了海面高度反演模型的设计与优化流程。内容融合传统物理模型与机器学习方法,采用随机森林和卷积神经网络,对比评估北斗与GPS反射信号的反演性能:BDS物理模型最大MAE约3.0m,优于GPS的约5.0m,而机器学习使两者均降至约0.4m,验证了国产卫星GNSS-R数据的测高潜力。资源包为1个PDF文件,约889KB,完整呈现数据预处理、物理模型实现、机器学习训练与评估、数据质量控制、误差分析及空间分布可视化等环节,并附可运行代码与逐段解释,便于读者从理论到实践掌握反演框架。目前已有124人学习,适合希望快速复现论文、对比不同卫星系统精度或推广FY-3E海面测高应用的研究者参考。
1. 从FY-3E的GNSS-R数据里把海面高度“捞”出来:这条技术路线到底值不值得跟
国产卫星FY-3E上搭载的GNSS-R接收机,本质上是在做一件“捡漏”的事——它不主动发射信号,而是接收导航卫星打到海面后反射回来的那一路信号。海面越平静、高度越低,反射信号的相关功率波形就越尖锐;海面越粗糙、高度越高,波形就被“抹平”得越厉害。海面高度反演要做的,就是从这条被海面“揉搓”过的波形里,把镜面反射点相对接收机的高度信息给还原出来。传统做法靠物理模型,比如用波形前沿斜率、峰值功率比这些几何特征去拟合;但物理模型在近岸、高海况、低仰角这些场景下经常翻车,误差能到米级。机器学习进来之后,思路变成“让模型自己从波形里找规律”,把物理特征和原始波形一起喂进去,用数据驱动的方式去修正物理模型的残差。这篇要拆的,就是FY-3E星载GNSS-R海面高度反演模型从物理基线到机器学习融合的完整落地路径——适合做遥感反演、卫星数据处理、或者想拿国产卫星数据练手机器学习回归任务的人。读完你能自己搭一套从原始波形到高度产品的流水线,知道每一步的参数怎么设、哪里容易踩坑、以及这套方案到底能压到多少精度。
2. FY-3E GNSS-R海面高度反演的物理底座:从波形到镜面反射点
2.1 星载GNSS-R的几何关系与延迟-多普勒图
FY-3E的GNSS-R接收机输出的核心数据是延迟-多普勒图(DDM),你可以把它理解成一张二维相关功率谱:横轴是码延迟,纵轴是多普勒频移,每个格点上的值代表对应延迟和多普勒下的相关功率。海面高度信息藏在哪?藏在镜面反射点对应的延迟上。镜面反射点是导航卫星发射信号经海面镜面反射后到达接收机的最短路径点,理论上这个点的延迟最小、功率最集中。实际中因为海面粗糙,能量会散开,但镜面点附近的波形前沿仍然是最陡的。
几何关系上,镜面反射点的位置由发射机、接收机、海面三者共同决定。FY-3E在轨高度约836 km,导航卫星高度约20200 km,这个几何下镜面反射点通常落在接收机星下点附近几十到几百公里范围内。你要算海面高度,先得把镜面反射点的经纬度算出来,再把DDM上镜面点对应的延迟量提取出来,最后通过几何关系反推高度。这一步的精度直接决定后续所有环节的天花板。
常见做法是先用WGS-84椭球模型做几何迭代:给定发射机和接收机的位置,假设一个海面高度初值,计算镜面反射点,再根据反射路径延迟修正高度,迭代到收敛。这个迭代一般3到5次就能稳定,收敛判据用高度变化小于0.1 m。
import numpy as np from scipy.optimize import fsolve def specular_point_geometry(tx_pos, rx_pos, h_guess=0.0, max_iter=10, tol=0.1): """ 迭代求解镜面反射点位置与海面高度 tx_pos: 发射机ECEF坐标 (m) rx_pos: 接收机ECEF坐标 (m) h_guess: 海面高度初值 (m) tol: 收敛阈值 (m) """ h = h_guess for i in range(max_iter): # 假设海面为椭球面,高度h,计算镜面反射点 # 这里用简化几何:镜面点在tx-rx连线与椭球面的交点附近 # 实际工程中会用更精确的迭代,此处展示逻辑框架 midpoint = (tx_pos + rx_pos) / 2.0 # 沿连线方向搜索使反射路径延迟最小的点 # 简化处理:取中点投影到椭球面 norm = np.linalg.norm(midpoint) earth_radius = 6371000.0 + h sp_point = midpoint / norm * earth_radius # 计算反射路径延迟 path_tx_sp = np.linalg.norm(tx_pos - sp_point) path_sp_rx = np.linalg.norm(sp_point - rx_pos) total_path = path_tx_sp + path_sp_rx # 根据延迟反推高度修正量(简化) # 实际中需要结合DDM上提取的延迟量 h_new = h + 0.01 * (total_path - np.linalg.norm(tx_pos - rx_pos)) if abs(h_new - h) < tol: break h = h_new return sp_point, h # 参数说明: # tx_pos/rx_pos 用ECEF坐标,单位米,FY-3E轨道可从TLE两行根数推算 # h_guess 一般给0,近岸区域可给潮汐模型初值 # tol 设0.1 m,再小意义不大,因为DDM延迟分辨率有限上面这段代码展示的是几何迭代的骨架逻辑。实际工程中,镜面反射点的搜索不是简单取中点投影,而是要在椭球面上做二维优化,使反射路径延迟最小。FY-3E的DDM延迟分辨率对应到海面高度大约在0.5到1米量级,所以几何迭代的收敛阈值设0.1米足够。参数上,发射机位置从导航卫星星历获取,接收机位置从FY-3E的精密定轨产品获取,这两个的精度直接决定几何计算的底噪。
2.2 从DDM提取延迟观测量:前沿斜率与峰值定位
拿到DDM之后,下一步是把镜面反射点对应的延迟量提出来。最直接的方法是找DDM上的功率峰值位置,但海面粗糙时峰值会偏移,而且多普勒维度上也有展宽。更稳的做法是取镜面反射点附近的多普勒切片,在延迟维度上做波形前沿拟合。
具体操作:以镜面反射点对应的多普勒频移为中心,取±500 Hz范围内的多普勒行做非相干平均,得到一条一维延迟波形。然后对波形前沿做线性拟合,取拟合线与噪声基线的交点作为延迟观测量。这个交点对应的延迟比峰值延迟更稳定,因为前沿斜率受海面粗糙度影响相对小。
def extract_delay_from_ddm(ddm, doppler_center_idx, doppler_half_width=5, noise_floor_percentile=10): """ 从DDM中提取延迟观测量 ddm: 二维数组 (延迟 x 多普勒) doppler_center_idx: 镜面反射点对应的多普勒索引 doppler_half_width: 多普勒平均半宽(索引数) noise_floor_percentile: 噪声基线估计百分位 """ # 多普勒维度平均 ddm_slice = ddm[:, doppler_center_idx - doppler_half_width: doppler_center_idx + doppler_half_width + 1] waveform = np.mean(ddm_slice, axis=1) # 估计噪声基线 noise_floor = np.percentile(waveform, noise_floor_percentile) # 找波形前沿:从峰值往回找第一个低于噪声基线的点 peak_idx = np.argmax(waveform) leading_edge_idx = peak_idx for i in range(peak_idx, 0, -1): if waveform[i] < noise_floor: leading_edge_idx = i break # 对前沿做线性拟合 fit_range = slice(max(0, leading_edge_idx - 2), leading_edge_idx + 3) x = np.arange(fit_range.start, fit_range.stop) y = waveform[fit_range] coeffs = np.polyfit(x, y, 1) # 拟合线与噪声基线的交点 delay_idx = (noise_floor - coeffs[1]) / coeffs[0] return delay_idx, waveform, noise_floor # 参数说明: # doppler_half_width 设5,对应约±500 Hz,FY-3E的DDM多普勒分辨率约100 Hz # noise_floor_percentile 设10,取波形最低10%分位作为噪声估计 # 前沿拟合范围取前沿点前后各2个格点,太宽会引入非线性这段代码的关键在于前沿拟合的范围控制。取太宽,前沿的非线性部分会拉偏斜率;取太窄,噪声会让拟合不稳定。我一般会先可视化几条波形确认前沿位置,再定拟合窗口。另外,噪声基线的估计要用波形最左端(延迟最小、无信号区域)的数据,而不是整条波形的百分位,否则海面回波强的时候会把基线抬高。
延迟观测量提取出来之后,结合2.1节的几何关系,就能算出海面高度。物理模型的反演精度在开阔海域、中低海况下能到0.3到0.5米,但近岸和高海况下会恶化到1米以上。这就是机器学习要补位的地方。
3. 机器学习融合方案:把物理特征和原始波形一起喂给模型
3.1 特征工程:物理特征、波形统计量与辅助参数
机器学习模型不是凭空学,输入特征的设计直接决定上限。我一般把特征分成三组:物理几何特征、波形形态特征、辅助环境参数。
物理几何特征包括:镜面反射点经纬度、发射机与接收机的几何距离、入射角、镜面反射点处的海面高度初值(来自物理模型)、DDM延迟观测量。这些是物理模型已经用到的,放进模型里相当于给了一个强基线。
波形形态特征包括:波形前沿斜率、峰值功率、后沿衰减率、波形对称性、多普勒展宽、相关功率的信噪比。这些是物理模型没有充分利用的,机器学习可以从里面挖出与海面高度相关的隐藏信息。
辅助环境参数包括:有效波高(从其他源获取或从波形反演)、风速、海面温度。这些不是必须的,但如果有,能显著提升模型在高海况下的表现。FY-3E本身不直接测这些,但可以通过微波辐射计或者再分析数据匹配。
def build_features(ddm, geom_info, aux_data=None): """ 构建机器学习输入特征向量 ddm: 二维数组 geom_info: 字典,包含几何参数 aux_data: 字典,辅助环境参数,可为None """ features = [] # 物理几何特征 features.append(geom_info['sp_lat']) features.append(geom_info['sp_lon']) features.append(geom_info['incidence_angle']) features.append(geom_info['range_tx_rx']) features.append(geom_info['phys_height']) features.append(geom_info['delay_obs']) # 波形形态特征 waveform = np.mean(ddm, axis=1) peak_power = np.max(waveform) noise_floor = np.percentile(waveform[:10], 50) snr = 10 * np.log10(peak_power / noise_floor) features.append(snr) # 前沿斜率 peak_idx = np.argmax(waveform) leading_slope = (waveform[peak_idx] - waveform[max(0, peak_idx-3)]) / 3.0 features.append(leading_slope) # 后沿衰减率 trailing_slope = (waveform[min(len(waveform)-1, peak_idx+3)] - waveform[peak_idx]) / 3.0 features.append(trailing_slope) # 多普勒展宽 doppler_profile = np.mean(ddm, axis=0) doppler_width = np.sum(doppler_profile > noise_floor) / len(doppler_profile) features.append(doppler_width) # 辅助参数 if aux_data is not None: features.append(aux_data.get('wave_height', 0.0)) features.append(aux_data.get('wind_speed', 0.0)) return np.array(features) # 参数说明: # sp_lat/sp_lon 单位度,范围-90到90和-180到180 # incidence_angle 单位度,FY-3E的GNSS-R入射角范围约0到70度 # range_tx_rx 单位米,典型值约2.1e7米 # phys_height 物理模型反演高度,单位米 # delay_obs 从DDM提取的延迟观测量,单位码片 # snr 单位dB,正常海面回波在5到15 dB之间 # wave_height/wind_speed 如果有,单位米和米/秒特征构建完之后,要做归一化。物理几何特征里的经纬度、距离这些量纲差异大,不归一化会让梯度下降很难收敛。我一般用RobustScaler而不是StandardScaler,因为波形特征里偶尔会有异常值,RobustScaler对异常值更稳。
3.2 模型选型:从梯度提升树到一维卷积网络
特征工程做完,模型选型上我试过三条路:梯度提升树(XGBoost/LightGBM)、一维卷积神经网络(1D-CNN)、以及物理模型残差学习。
梯度提升树适合特征维度不高、样本量中等的情况。FY-3E的GNSS-R数据,如果只取物理特征和波形统计量,特征维度在15到20之间,样本量在几万到几十万条,XGBoost跑起来很快,调参也直观。我一般先用它做基线,看特征重要性排序,确认哪些特征真正有用。
1D-CNN适合直接吃原始波形。把DDM沿延迟维度取多普勒平均后得到的一维波形,或者直接把DDM展平成一维向量,送进卷积层。CNN能自动学出前沿、峰值、后沿这些局部模式,不需要手工设计特征。但CNN需要更多样本,而且训练慢,调参玄学成分大。
物理模型残差学习是我最推荐的做法:先用物理模型算一个高度初值,然后让机器学习模型去预测物理模型与真实高度之间的残差。这样模型只需要学“物理模型哪里错了”,而不是从零学“高度是多少”,学习难度大幅降低,样本效率也高。
import xgboost as xgb from sklearn.model_selection import train_test_split from sklearn.preprocessing import RobustScaler from sklearn.metrics import mean_absolute_error # 假设 X 是特征矩阵 (n_samples, n_features),y 是真实海面高度 # phys_height 是物理模型反演高度 X_train, X_test, y_train, y_test, phys_train, phys_test = train_test_split( X, y, phys_height, test_size=0.2, random_state=42 ) # 归一化 scaler = RobustScaler() X_train_scaled = scaler.fit_transform(X_train) X_test_scaled = scaler.transform(X_test) # 残差目标 residual_train = y_train - phys_train residual_test = y_test - phys_test # XGBoost回归 model = xgb.XGBRegressor( n_estimators=500, max_depth=6, learning_rate=0.05, subsample=0.8, colsample_bytree=0.8, reg_alpha=0.1, reg_lambda=1.0, random_state=42 ) model.fit(X_train_scaled, residual_train) # 预测 residual_pred = model.predict(X_test_scaled) height_pred = phys_test + residual_pred mae = mean_absolute_error(y_test, height_pred) print(f"融合模型MAE: {mae:.3f} m") # 参数说明: # n_estimators 500,树的数量,再多容易过拟合 # max_depth 6,控制树复杂度,GNSS-R特征非线性强但样本有限,6层够用 # learning_rate 0.05,配合500棵树,学习率低一点更稳 # subsample/colsample_bytree 0.8,行采样和列采样,防过拟合 # reg_alpha/reg_lambda L1/L2正则,GNSS-R数据噪声大,正则要加这段代码的核心是残差学习框架。注意目标变量不是真实高度,而是真实高度减去物理模型高度。预测的时候再把物理模型高度加回去。这样做的好处是,即使机器学习模型预测残差有偏差,物理模型的高度初值也能兜底,不会出现离谱的预测。
参数上,n_estimators和learning_rate要配合调。我一般先用0.1的学习率跑200棵树看收敛,再降到0.05跑500棵。max_depth不要超过8,GNSS-R特征维度不高,树太深会记住噪声。subsample和colsample_bytree设0.8是经验值,数据量小的时候可以降到0.7。
3.3 训练集构建与时空交叉验证
GNSS-R数据有个坑:同一轨道的相邻样本高度相关性很强,如果随机划分训练集和测试集,测试集里的样本可能和训练集里的样本来自同一片海域、同一时间段,导致精度虚高。我踩过这个坑,随机划分下MAE能到0.15米,但按轨道划分后直接掉到0.35米。
正确的做法是按时空块划分:把数据按经纬度网格和时间窗口切块,确保训练集和测试集在空间和时间上都不重叠。具体操作:先按1度×1度网格聚合,再按天划分,训练集取80%的网格和天数,测试集取剩下的20%。这样测试集里的海况和地理区域都是训练集没见过的,评估结果更接近实际应用。
def spatiotemporal_split(lats, lons, times, test_ratio=0.2): """ 按时空块划分训练集和测试集 lats/lons: 纬度/经度数组 times: 时间数组(datetime) test_ratio: 测试集比例 """ # 网格化 grid_lat = np.floor(lats).astype(int) grid_lon = np.floor(lons).astype(int) grid_day = np.array([t.day for t in times]) # 构建时空块ID block_id = grid_lat * 10000 + grid_lon * 100 + grid_day unique_blocks = np.unique(block_id) n_test = int(len(unique_blocks) * test_ratio) test_blocks = np.random.choice(unique_blocks, n_test, replace=False) test_mask = np.isin(block_id, test_blocks) train_mask = ~test_mask return train_mask, test_mask # 参数说明: # 网格大小1度,对应约111 km,GNSS-R镜面反射点空间分辨率约25 km # 按天划分,确保测试集和训练集不在同一天 # test_ratio 0.2,测试集占20%的时空块时空交叉验证比随机划分严格得多,但更真实。如果你的模型在时空交叉验证下MAE能到0.3米以内,那在实际应用中大概率也能保持这个水平。如果随机划分和时空划分差距很大,说明模型过拟合了特定海域或特定时间段的海况,需要加正则或者减特征。
4. 避坑与排查:FY-3E GNSS-R反演里那些让你白干一天的坑
4.1 镜面反射点算偏了,后面全白搭
现象:物理模型反演高度系统性偏高或偏低,误差在特定区域特别大。
原因:镜面反射点的几何计算用了简化模型,没有考虑地球曲率和椭球扁率。FY-3E轨道高度836 km,镜面反射点距离星下点最远能到500 km以上,这个距离下地球曲率的影响不能忽略。
解决:用WGS-84椭球模型做精确几何迭代,迭代时把发射机、接收机、镜面反射点三者的ECEF坐标都算准。我一般会拿已知的GNSS-R站点数据做验证,确保几何计算误差小于0.1度。
4.2 DDM延迟分辨率不够,高度精度上不去
现象:物理模型反演高度在开阔海域也只能到0.5米,再想提升很难。
原因:FY-3E的DDM延迟分辨率对应到海面高度大约0.5到1米,这是硬件决定的。物理模型直接从延迟观测量算高度,精度天花板就在这。
解决:用机器学习融合波形形态特征,把延迟分辨率之外的信息挖出来。波形前沿斜率、后沿衰减率这些特征对高度变化敏感,但物理模型没有充分利用。融合模型可以把MAE压到0.3米左右。
4.3 训练集和测试集随机划分,精度虚高
现象:模型在测试集上MAE 0.15米,实际应用时误差0.5米以上。
原因:随机划分导致同一海域、同一时间段的样本同时出现在训练集和测试集,模型记住了这些样本的特征,而不是学到了泛化规律。
解决:按时空块划分,确保训练集和测试集在空间和时间上都不重叠。具体做法见3.3节。时空交叉验证下的精度才是真实精度。
4.4 辅助环境参数匹配不上,特征缺失
现象:模型在高海况下表现差,但不知道原因。
原因:有效波高、风速这些辅助参数来自其他数据源,时空匹配不上。FY-3E的GNSS-R镜面反射点位置和微波辐射计的视场不完全重合,直接插值会引入误差。
解决:匹配辅助参数时,设置时空窗口(比如时间±30分钟、空间±0.25度),取窗口内的平均值。如果窗口内没有辅助数据,该样本的辅助特征置为NaN,让XGBoost自己处理缺失值。不要用全局均值填充,会引入偏差。
4.5 模型过拟合物理模型残差,残差预测反而添乱
现象:融合模型在某些区域比纯物理模型还差。
原因:物理模型残差在近岸区域分布和开阔海域完全不同,模型如果没学好,残差预测会引入额外误差。
解决:在残差学习框架里加一个门控机制:当物理模型高度与辅助参数(如有效波高)一致性高时,信任物理模型;一致性低时,信任机器学习残差。具体实现可以用一个简单的加权:最终高度 = α × 物理高度 + (1-α) × (物理高度 + 残差预测),α根据海况动态调整。
5. 把MAE从0.35米压到0.25米:几个我反复验证过的技巧
5.1 波形重跟踪:把前沿拟合从线性换成Bi-Gaussian
前面2.2节用的是线性拟合前沿,简单但精度有限。我后来换成Bi-Gaussian拟合:用两个高斯函数的组合去拟合整个波形,一个描述前沿上升,一个描述后沿衰减。拟合出来的前沿中心位置比线性交点更稳,尤其是在高海况下。
from scipy.optimize import curve_fit def bi_gaussian(x, a1, mu1, sigma1, a2, mu2, sigma2): return a1 * np.exp(-(x - mu1)**2 / (2 * sigma1**2)) + \ a2 * np.exp(-(x - mu2)**2 / (2 * sigma2**2)) def fit_waveform_bi_gaussian(waveform): x = np.arange(len(waveform)) # 初值:前沿在峰值左侧,后沿在峰值右侧 peak_idx = np.argmax(waveform) p0 = [waveform[peak_idx], peak_idx - 2, 1.5, waveform[peak_idx] * 0.7, peak_idx + 3, 2.5] try: popt, _ = curve_fit(bi_gaussian, x, waveform, p0=p0, maxfev=5000) return popt except RuntimeError: return None # 参数说明: # a1/a2 两个高斯峰的幅度 # mu1/mu2 两个峰的中心位置,mu1对应前沿,mu2对应后沿 # sigma1/sigma2 两个峰的标准差,控制宽度 # 初值里mu1设峰值左侧2个格点,mu2设峰值右侧3个格点Bi-Gaussian拟合的关键是初值要给好。mu1和mu2的初值根据峰值位置偏移,sigma1和sigma2根据波形宽度估计。拟合失败时返回None,该样本丢弃或者回退到线性拟合。
5.2 多普勒维度加权:镜面点附近的行给更高权重
DDM的多普勒维度上,镜面反射点对应的多普勒行信噪比最高,远离镜面点的行噪声占比大。做多普勒平均时,不要等权平均,给镜面点附近的行更高权重。
def weighted_doppler_average(ddm, doppler_center_idx, sigma=3.0): """ 多普勒维度加权平均 sigma: 高斯权重标准差,控制权重衰减速度 """ n_doppler = ddm.shape[1] weights = np.exp(-(np.arange(n_doppler) - doppler_center_idx)**2 / (2 * sigma**2)) weights = weights / np.sum(weights) waveform = np.dot(ddm, weights) return waveform # 参数说明: # sigma 设3.0,对应约300 Hz,FY-3E多普勒分辨率约100 Hz # sigma太小,只用镜面点附近几行,噪声大 # sigma太大,等效于等权平均,失去加权意义sigma设3.0是我试出来的经验值。太小了波形噪声大,太大了和等权平均没区别。你可以拿几条典型波形对比一下,看加权后的前沿是否更清晰。
5.3 残差学习的门控加权:让物理模型和机器学习各司其职
前面4.5节提到门控机制,这里给一个具体实现。核心思路是用有效波高和入射角判断当前海况,海况低时信任物理模型,海况高时信任机器学习。
def gated_fusion(phys_height, ml_residual, wave_height, incidence_angle): """ 门控融合物理模型和机器学习残差 wave_height: 有效波高 (m) incidence_angle: 入射角 (度) """ # 海况权重:波高越大,越信任机器学习 wh_weight = np.clip(wave_height / 4.0, 0.0, 1.0) # 入射角权重:入射角越大,物理模型越不准,越信任机器学习 ia_weight = np.clip((incidence_angle - 20) / 50.0, 0.0, 1.0) # 综合权重 alpha = 1.0 - 0.5 * (wh_weight + ia_weight) alpha = np.clip(alpha, 0.3, 0.9) return alpha * phys_height + (1 - alpha) * (phys_height + ml_residual) # 参数说明: # wave_height 从辅助数据获取,没有时用0.5米默认值 # incidence_angle 从几何计算获取,FY-3E范围0到70度 # alpha下限0.3,保证物理模型至少占30%权重 # alpha上限0.9,保证机器学习至少占10%权重这个门控函数我调了很久,最终定下来alpha在0.3到0.9之间。波高4米以上时,机器学习权重占主导;入射角20度以下时,物理模型权重占主导。实际跑下来,门控融合比固定权重的MAE低了0.03到0.05米。
5.4 验证方法:留一轨交叉验证比留一时间段更严格
时空交叉验证里,按天划分已经比随机划分严格很多。但如果你想再狠一点,用留一轨交叉验证:每次留一整条FY-3E轨道做测试,其余轨道做训练。这样测试集里的海况、地理区域、卫星几何都和训练集完全不同,评估结果最接近实际业务。
留一轨交叉验证的计算量很大,FY-3E每天约14到15轨,做一轮要训练14到15个模型。我一般只在最终评估时做一轮,调参阶段用按天划分就够了。
最后说个我自己的习惯:每次跑完模型,我都会把预测残差按经纬度画个热力图,看看误差在哪些区域集中。如果误差在近岸集中,说明辅助参数匹配有问题;如果误差在高纬度集中,说明训练集里高纬度样本太少。这个图比任何指标都直观,能帮你快速定位问题。希望帮到你。
本文还有配套的精品资源,点击获取