198 lines
7.3 KiB
Plaintext
198 lines
7.3 KiB
Plaintext
================================================================================
|
||
地形校正模块 - 已知问题与改进方向
|
||
================================================================================
|
||
|
||
【问题描述】
|
||
-----------
|
||
当前版本的 slope_aspectV2.py 在处理推扫式无人机高光谱图像时，存在一个关键假设：
|
||
整幅图像使用单一的太阳天顶角和方位角来计算 cosine_i（入射角余弦）。
|
||
|
||
然而，实际数据采集情况是：
|
||
- 高光谱图像由3个架次的无人机飞行拼接而成
|
||
- 每个架次的飞行时间不同（通常间隔数分钟至数十分钟）
|
||
- 因此，不同扫描线（或不同图像区域）对应的太阳角度是不同的
|
||
- 但目前代码仅基于DEM中心点坐标计算一组太阳角度，并应用于整幅图像
|
||
|
||
【技术影响】
|
||
-----------
|
||
这个简化会导致以下问题：
|
||
|
||
1. 时间跨度误差
|
||
- 假设3个架次间隔30分钟，太阳方位角可能变化7-8度（取决于纬度）
|
||
- 在日出/日落时段，太阳高度角变化可达10-15度
|
||
- 这会导致不同架次交界处的 cosine_i 出现系统性偏差
|
||
|
||
2. 地形校正不准确
|
||
- 东向坡面：过早架次会低估照明，过晚架次会高估照明
|
||
- 西向坡面：与东向坡面相反
|
||
- 交界区域会出现明显的"拼接缝"现象
|
||
|
||
3. BRDF校正累积误差
|
||
- cosine_i 是许多BRDF模型的关键输入参数
|
||
- 角度误差会传播到后续的BRDF地形校正结果
|
||
|
||
【当前代码实现】
|
||
---------------
|
||
文件: slope/slope_aspectV2.py
|
||
|
||
当前流程：
|
||
1. 读取DEM -> 2. 计算坡度/坡向 -> 3. 获取DEM中心点坐标(lat, lon)
|
||
4. 根据单一时间计算一组太阳角度(solar_zn, solar_az)
|
||
5. 将这组角度应用于整幅图像的所有像素
|
||
|
||
关键代码段（约110-130行）：
|
||
```python
|
||
# 获取DEM中心点坐标
|
||
latitude, longitude = get_dem_center_coords(dem) # 仅一个中心点！
|
||
|
||
# 计算太阳角度（仅一组）
|
||
solar_zn_deg, solar_az_deg = calc_solar_angles(latitude, longitude, observation_time)
|
||
|
||
# 应用于整幅图像
|
||
cosine_i = calc_cosine_i(solar_zn_rad, solar_az_rad, aspect_rad, slope_rad)
|
||
```
|
||
|
||
【需要的改进方向】
|
||
-----------------
|
||
接手人员需要进行以下修改：
|
||
|
||
方案A: 分架次处理（推荐）
|
||
----------------------------
|
||
如果已知3个架次的地理分界线或时间范围：
|
||
|
||
1. 输入格式修改
|
||
- 接收3组时间参数（或从飞行日志自动提取）
|
||
- 例如: -t1 "2024-06-15 08:30:00" -t2 "2024-06-15 09:00:00" -t3 "2024-06-15 09:30:00"
|
||
|
||
2. 空间分区处理
|
||
- 根据每架次的覆盖范围（可从高光谱图像的GLT或地理坐标确定）
|
||
- 将DEM分割为3个区域
|
||
- 每个区域使用对应架次的时间计算太阳角度
|
||
|
||
3. 修改后的伪代码示例：
|
||
```python
|
||
# 定义3个架次的时间和地理范围
|
||
flight_segments = [
|
||
{'time': '2024-06-15 08:30:00', 'row_range': (0, 1200)}, # 架次1: 行0-1199
|
||
{'time': '2024-06-15 09:00:00', 'row_range': (1200, 2800)}, # 架次2: 行1200-2799
|
||
{'time': '2024-06-15 09:30:00', 'row_range': (2800, 4000)}, # 架次3: 行2800-3999
|
||
]
|
||
|
||
# 初始化cosine_i数组
|
||
cosine_i = np.zeros_like(slope_data)
|
||
|
||
# 分段计算
|
||
for segment in flight_segments:
|
||
dt = parse_datetime(segment['time'])
|
||
solar_zn, solar_az = calc_solar_angles(lat, lon, dt)
|
||
|
||
y1, y2 = segment['row_range']
|
||
cosine_i[y1:y2, :] = calc_cosine_i(
|
||
np.deg2rad(solar_zn), np.deg2rad(solar_az),
|
||
aspect_rad[y1:y2, :], slope_rad[y1:y2, :]
|
||
)
|
||
```
|
||
|
||
方案B: 逐行/逐像素太阳角度（更精确）
|
||
-------------------------------------
|
||
如果需要更高精度（考虑太阳角度在单架次内的缓慢变化）：
|
||
|
||
1. 输入要求
|
||
- 高光谱图像的每一行对应的时间戳（可从IMU/GPS数据提取）
|
||
- 或者：每个像素的精确UTC时间
|
||
|
||
2. 逐行计算
|
||
```python
|
||
# 假设times是一个与图像行数相同长度的一维数组
|
||
# times[i] = 第i行的观测时间（datetime对象）
|
||
|
||
for row in range(n_rows):
|
||
dt = observation_times[row] # 该行的具体时间
|
||
solar_zn, solar_az = calc_solar_angles(lat, lon, dt)
|
||
|
||
cosine_i[row, :] = calc_cosine_i(
|
||
np.deg2rad(solar_zn), np.deg2rad(solar_az),
|
||
aspect_rad[row, :], slope_rad[row, :]
|
||
)
|
||
```
|
||
|
||
3. 性能优化
|
||
- 由于pvlib计算较慢，可对时间相近的行批量计算
|
||
- 例如：每10行计算一次太阳角度，中间行插值
|
||
|
||
方案C: 从飞行日志自动提取（最自动化）
|
||
--------------------------------------
|
||
如果飞行平台记录了GPS/IMU数据：
|
||
|
||
1. 输入数据
|
||
- 读取POS数据（位置、姿态、时间）
|
||
- 或从高光谱图像的GLT（地理位置查找表）提取
|
||
|
||
2. 建立时间-位置映射
|
||
- 根据无人机位置和时间建立每行图像的观测时间
|
||
|
||
3. 分区域计算cosine_i
|
||
|
||
【数据结构参考】
|
||
---------------
|
||
典型的推扫式高光谱数据结构：
|
||
|
||
高光谱图像 (lines, samples, bands)
|
||
- lines: 扫描线数（飞行方向的像素数）
|
||
- samples: 每行的像素数（垂直飞行方向）
|
||
- 第i行对应无人机的第i个位置
|
||
|
||
关键信息来源：
|
||
1. 图像头文件 (.hdr) - 通常包含起始时间和积分时间
|
||
2. GLT文件 (Geographic Lookup Table) - 每个像素的地理坐标
|
||
3. 飞行日志/POS数据 - 每行对应的时间戳
|
||
4. 帧辅助数据 (Frame Auxiliary Data) - 部分传感器自带每帧时间
|
||
|
||
【建议的输入接口修改】
|
||
----------------------
|
||
当前命令行接口：
|
||
python slope_aspectV2.py -i dem.tif -o output.dat -t "2024-06-15 10:00:00"
|
||
|
||
建议扩展为：
|
||
# 方式1: 手动指定3个架次的时间和行范围
|
||
python slope_aspectV2.py -i dem.tif -o output.dat \
|
||
--segments "0,1200,2024-06-15T08:30:00" \
|
||
--segments "1200,2800,2024-06-15T09:00:00" \
|
||
--segments "2800,4000,2024-06-15T09:30:00"
|
||
|
||
# 方式2: 从飞行日志文件自动读取
|
||
python slope_aspectV2.py -i dem.tif -o output.dat --time-log flight_log.csv
|
||
|
||
# 方式3: 从高光谱图像头文件提取时间信息
|
||
python slope_aspectV2.py -i dem.tif -o output.dat --hyperspectral-ref image.hdr
|
||
|
||
【验证方法】
|
||
-----------
|
||
修改后应验证以下指标：
|
||
|
||
1. 连续性检查
|
||
- 目视检查3个架次交界处的cosine_i是否连续
|
||
- 统计交界行前后10行的cosine_i差异
|
||
|
||
2. 物理合理性
|
||
- 东向坡的cosine_i应随时间增加而增加（上午飞行）
|
||
- 西向坡的cosine_i应随时间增加而减小（上午飞行）
|
||
|
||
3. 与GLT对比
|
||
- 如果高光谱有GLT，可根据每像素的太阳角度公式验证
|
||
|
||
【相关文件】
|
||
-----------
|
||
- slope_aspectV2.py - 当前主程序（需要修改）
|
||
- angle_compute.py - 角度计算工具（可能包含相关逻辑）
|
||
- 高光谱图像的.hdr文件 - 可能包含时间元数据
|
||
- 飞行日志/POS数据 - 外部时间参考
|
||
|
||
【优先级建议】
|
||
--------------
|
||
1. 【高】首先实现方案A（分架次处理）- 工作量适中，效果显著
|
||
2. 【中】评估是否需要方案B（逐行）- 取决于太阳角度变化幅度
|
||
3. 【低】长期可考虑方案C（全自动）- 需要整合飞行数据解析
|
||
|
||
简单来说在计算的时候，需要加入角度文件，进行逐像素的cosi计算
|