157_环境科学：气候数据解读 - LLM驱动的卫星数据多模态融合与时空分析技术研究

class SatelliteDataProcessor:
    def __init__(self, config=None):
        self.config = config or {
            "image_resolution": 256,
            "time_window": 7,  # 天
            "normalization_method": "z_score",
            "cloud_threshold": 0.8,
            "interpolation_method": "linear"
        }
        self.data_cache = {}
        print("卫星数据处理器初始化完成")
    
    def load_satellite_image(self, image_path, band_selection=None):
        """加载卫星图像并选择特定波段"""
        try:
            import rasterio
            import numpy as np
            
            with rasterio.open(image_path) as src:
                # 如果指定了波段选择
                if band_selection:
                    bands = src.read(band_selection)
                else:
                    bands = src.read()
                
                # 获取元数据
                meta = src.meta
                transform = src.transform
                crs = src.crs
            
            # 转换数据格式
            # 将(bands, height, width)转换为(height, width, bands)
            if len(bands.shape) == 3:
                bands = np.transpose(bands, (1, 2, 0))
            
            return {
                "data": bands,
                "meta": meta,
                "transform": transform,
                "crs": crs
            }
            
        except Exception as e:
            print(f"加载卫星图像错误: {e}")
            return None
    
    def normalize_image(self, image_data):
        """图像数据标准化"""
        import numpy as np
        
        data = image_data["data"].copy()
        
        if self.config["normalization_method"] == "z_score":
            # 计算每个波段的均值和标准差
            for b in range(data.shape[2]):
                band_data = data[:,:,b]
                mean_val = np.nanmean(band_data)
                std_val = np.nanstd(band_data)
                if std_val > 0:
                    data[:,:,b] = (band_data - mean_val) / std_val
        
        elif self.config["normalization_method"] == "minmax":
            # 最小-最大标准化
            for b in range(data.shape[2]):
                band_data = data[:,:,b]
                min_val = np.nanmin(band_data)
                max_val = np.nanmax(band_data)
                if max_val > min_val:
                    data[:,:,b] = (band_data - min_val) / (max_val - min_val)
        
        result = image_data.copy()
        result["data"] = data
        return result
    
    def cloud_masking(self, image_data, cloud_prob=None):
        """云掩膜处理"""
        import numpy as np
        
        if cloud_prob is None:
            # 如果没有提供云概率数据，使用默认阈值
            # 这里仅作为示例，实际应用需要更复杂的云检测算法
            return image_data
        
        data = image_data["data"].copy()
        
        # 应用云掩膜
        mask = cloud_prob < self.config["cloud_threshold"]
        
        for b in range(data.shape[2]):
            data[:,:,b] = np.where(mask, data[:,:,b], np.nan)
        
        result = image_data.copy()
        result["data"] = data
        result["cloud_mask"] = mask
        return result
    
    def resample_image(self, image_data, target_resolution):
        """重采样图像到目标分辨率"""
        try:
            from rasterio.enums import Resampling
            import rasterio.transform
            import numpy as np
            
            data = image_data["data"]
            transform = image_data["transform"]
            
            # 计算缩放因子
            src_height, src_width = data.shape[:2]
            scale_x = target_resolution / transform.a
            scale_y = -target_resolution / transform.e  # 负号因为y方向通常向下
            
            # 计算新的尺寸
            new_height = int(src_height * scale_y)
            new_width = int(src_width * scale_x)
            
            # 创建新的变换
            new_transform = rasterio.transform.from_bounds(
                transform.xoff,
                transform.yoff + src_height * transform.e,
                transform.xoff + src_width * transform.a,
                transform.yoff,
                new_width,
                new_height
            )
            
            # 重采样每个波段
            resampled_data = np.zeros((new_height, new_width, data.shape[2]), dtype=np.float32)
            
            for b in range(data.shape[2]):
                from scipy.ndimage import zoom
                resampled_band = zoom(data[:,:,b], (scale_y, scale_x), 
                                     order=1, mode='constant', cval=np.nan)
                resampled_data[:,:,b] = resampled_band
            
            result = image_data.copy()
            result["data"] = resampled_data
            result["transform"] = new_transform
            result["meta"]["height"] = new_height
            result["meta"]["width"] = new_width
            result["meta"]["transform"] = new_transform
            
            return result
            
        except Exception as e:
            print(f"图像重采样错误: {e}")
            return image_data
    
    def extract_time_series(self, point_coords, data_cube, time_dimension):
        """从数据立方体中提取时间序列数据"""
        import numpy as np
        from rasterio.transform import rowcol
        
        try:
            # 转换地理坐标到图像坐标
            row, col = rowcol(data_cube["transform"], point_coords[0], point_coords[1])
            
            # 检查坐标是否在图像范围内
            if (0 <= row < data_cube["data"].shape[1] and 
                0 <= col < data_cube["data"].shape[2]):
                
                # 提取时间序列
                time_series = data_cube["data"][:, row, col].copy()
                
                return {
                    "coordinates": point_coords,
                    "time_series": time_series,
                    "time_dimension": time_dimension
                }
            else:
                print(f"坐标({point_coords})超出图像范围")
                return None
                
        except Exception as e:
            print(f"提取时间序列错误: {e}")
            return None
    
    def preprocess_for_llm(self, image_data, time_series_data=None, metadata=None):
        """为LLM准备多模态输入数据"""
        processed_data = {
            "image_features": None,
            "time_series_features": None,
            "metadata_text": None
        }
        
        # 处理图像数据
        if image_data is not None:
            # 归一化
            normalized = self.normalize_image(image_data)
            
            # 调整分辨率
            resampled = self.resample_image(normalized, self.config["image_resolution"])
            
            # 提取图像特征（简化版，实际应用中可能使用CNN等提取特征）
            processed_data["image_features"] = {
                "data": resampled["data"].flatten()[:10000],  # 简化示例
                "shape": resampled["data"].shape,
                "transform": resampled["transform"].to_gdal()
            }
        
        # 处理时间序列数据
        if time_series_data is not None:
            import numpy as np
            
            ts_data = time_series_data["time_series"]
            
            # 计算统计特征
            if not np.all(np.isnan(ts_data)):
                processed_data["time_series_features"] = {
                    "mean": np.nanmean(ts_data),
                    "std": np.nanstd(ts_data),
                    "min": np.nanmin(ts_data),
                    "max": np.nanmax(ts_data),
                    "trend": np.polyfit(range(len(ts_data)), 
                                       np.nan_to_num(ts_data), 1)[0] if len(ts_data) > 1 else 0,
                    "first_value": ts_data[0] if not np.isnan(ts_data[0]) else 0,
                    "last_value": ts_data[-1] if not np.isnan(ts_data[-1]) else 0,
                    "length": len(ts_data)
                }
        
        # 处理元数据
        if metadata is not None:
            # 将元数据转换为文本描述
            metadata_text = "卫星数据元信息：\n"
            for key, value in metadata.items():
                metadata_text += f"{key}: {value}\n"
            processed_data["metadata_text"] = metadata_text
        
        return processed_data

class ImageFeatureExtractor:
    def __init__(self, backbone="resnet50", pretrained=True):
        self.backbone = backbone
        self.pretrained = pretrained
        self.model = None
        self.feature_dim = None
        self._initialize_model()
    
    def _initialize_model(self):
        """初始化特征提取模型"""
        try:
            import torch
            import torchvision.models as models
            from torchvision import transforms
            
            # 加载预训练的骨干网络
            if self.backbone == "resnet50":
                self.model = models.resnet50(pretrained=self.pretrained)
                # 移除最后的分类层
                self.model = torch.nn.Sequential(*(list(self.model.children())[:-1]))
                self.feature_dim = 2048
            elif self.backbone == "efficientnet_b0":
                self.model = models.efficientnet_b0(pretrained=self.pretrained)
                self.model = torch.nn.Sequential(*(list(self.model.children())[:-1]))
                self.feature_dim = 1280
            elif self.backbone == "densenet121":
                self.model = models.densenet121(pretrained=self.pretrained)
                self.model = self.model.features
                self.feature_dim = 1024
            else:
                raise ValueError(f"不支持的骨干网络: {self.backbone}")
            
            # 设置为评估模式
            self.model.eval()
            
            # 定义图像预处理变换
            self.transform = transforms.Compose([
                transforms.ToPILImage(),
                transforms.Resize((224, 224)),
                transforms.ToTensor(),
                transforms.Normalize(
                    mean=[0.485, 0.456, 0.406],
                    std=[0.229, 0.224, 0.225]
                )
            ])
            
            print(f"图像特征提取器初始化完成，骨干网络：{self.backbone}")
            
        except Exception as e:
            print(f"初始化特征提取模型错误: {e}")
    
    def extract_features(self, image_data, use_gpu=True):
        """提取图像特征"""
        import torch
        import numpy as np
        
        try:
            if self.model is None:
                self._initialize_model()
            
            data = image_data["data"].copy()
            
            # 处理多波段数据
            if data.shape[2] > 3:
                # 对于多波段数据，选择前3个波段或进行波段融合
                # 这里采用简单的PCA降维示例
                data_reshaped = data.reshape(-1, data.shape[2])
                
                # 移除NaN值
                valid_indices = ~np.isnan(data_reshaped).any(axis=1)
                valid_data = data_reshaped[valid_indices]
                
                if len(valid_data) > 0:
                    from sklearn.decomposition import PCA
                    pca = PCA(n_components=3)
                    reduced_data = pca.fit_transform(valid_data)
                    
                    # 重建图像
                    result = np.zeros_like(data_reshaped)
                    result[valid_indices] = reduced_data
                    result = result.reshape(data.shape[0], data.shape[1], 3)
                    
                    # 归一化到0-1范围
                    result = (result - result.min()) / (result.max() - result.min() + 1e-8)
                    data = result
                else:
                    # 如果没有有效数据，使用随机噪声
                    data = np.random.rand(data.shape[0], data.shape[1], 3)
            elif data.shape[2] == 1:
                # 单波段数据，复制到3个通道
                data = np.repeat(data, 3, axis=2)
            
            # 确保数据在0-1范围内
            data_min = data.min()
            data_max = data.max()
            if data_max > data_min:
                data = (data - data_min) / (data_max - data_min)
            
            # 转换为RGB格式（0-255）
            data = (data * 255).astype(np.uint8)
            
            # 应用预处理变换
            input_tensor = self.transform(data)
            input_tensor = input_tensor.unsqueeze(0)  # 添加批次维度
            
            # 移至GPU（如果可用）
            device = torch.device("cuda" if torch.cuda.is_available() and use_gpu else "cpu")
            self.model = self.model.to(device)
            input_tensor = input_tensor.to(device)
            
            # 提取特征
            with torch.no_grad():
                features = self.model(input_tensor)
            
            # 转换为numpy数组
            features = features.squeeze().cpu().numpy()
            
            return {
                "features": features,
                "feature_dim": self.feature_dim,
                "backbone": self.backbone
            }
            
        except Exception as e:
            print(f"提取图像特征错误: {e}")
            return None

class TimeSeriesFeatureExtractor:
    def __init__(self, config=None):
        self.config = config or {
            "window_size": 7,
            "statistical_features": True,
            "fft_features": True,
            "wavelet_features": False,
            "max_freq_components": 5
        }
    
    def extract_features(self, time_series_data):
        """提取时间序列特征"""
        import numpy as np
        
        try:
            ts = time_series_data["time_series"]
            
            # 移除NaN值
            valid_ts = ts[~np.isnan(ts)]
            
            if len(valid_ts) < 3:
                print("有效数据点不足，无法提取特征")
                return None
            
            features = {}
            
            # 基本统计特征
            if self.config["statistical_features"]:
                features.update({
                    "mean": np.mean(valid_ts),
                    "std": np.std(valid_ts),
                    "min": np.min(valid_ts),
                    "max": np.max(valid_ts),
                    "median": np.median(valid_ts),
                    "range": np.max(valid_ts) - np.min(valid_ts),
                    "skewness": self._calculate_skewness(valid_ts),
                    "kurtosis": self._calculate_kurtosis(valid_ts),
                    "percentiles_25": np.percentile(valid_ts, 25),
                    "percentiles_75": np.percentile(valid_ts, 75),
                    "percentiles_90": np.percentile(valid_ts, 90),
                    "zero_crossing_rate": self._calculate_zero_crossing_rate(valid_ts),
                    "trend": self._calculate_trend(valid_ts)
                })
            
            # FFT频域特征
            if self.config["fft_features"]:
                fft_features = self._extract_fft_features(valid_ts)
                features.update(fft_features)
            
            # 小波特征（可选）
            if self.config["wavelet_features"]:
                wavelet_features = self._extract_wavelet_features(valid_ts)
                features.update(wavelet_features)
            
            # 时序模式特征
            features.update({
                "autocorr_lag1": self._calculate_autocorrelation(valid_ts, lag=1),
                "autocorr_lag5": self._calculate_autocorrelation(valid_ts, lag=5),
                "seasonality_strength": self._detect_seasonality(valid_ts),
                "entropy": self._calculate_sample_entropy(valid_ts)
            })
            
            return features
            
        except Exception as e:
            print(f"提取时间序列特征错误: {e}")
            return None
    
    def _calculate_skewness(self, data):
        """计算偏度"""
        import scipy.stats as stats
        return stats.skew(data)
    
    def _calculate_kurtosis(self, data):
        """计算峰度"""
        import scipy.stats as stats
        return stats.kurtosis(data)
    
    def _calculate_zero_crossing_rate(self, data):
        """计算过零率"""
        zero_crossings = np.where(np.diff(np.signbit(data)))[0]
        return len(zero_crossings) / len(data)
    
    def _calculate_trend(self, data):
        """计算线性趋势"""
        x = np.arange(len(data))
        slope, _, _, _, _ = np.polyfit(x, data, 1, full=True)
        return slope
    
    def _extract_fft_features(self, data):
        """提取FFT频域特征"""
        import numpy as np
        
        # 应用汉宁窗
        window = np.hanning(len(data))
        windowed_data = data * window
        
        # 计算FFT
        fft_values = np.fft.fft(windowed_data)
        frequencies = np.fft.fftfreq(len(data))
        
        # 只考虑正频率
        positive_freq_indices = frequencies > 0
        fft_values = fft_values[positive_freq_indices]
        frequencies = frequencies[positive_freq_indices]
        
        # 计算幅度谱
        magnitude = np.abs(fft_values)
        
        # 提取主要频率成分
        num_components = min(self.config["max_freq_components"], len(magnitude))
        top_indices = np.argsort(magnitude)[-num_components:][::-1]
        
        fft_features = {
            "fft_total_power": np.sum(magnitude**2),
            "fft_peak_power": np.max(magnitude),
            "fft_peak_frequency": frequencies[top_indices[0]] if top_indices.size > 0 else 0
        }
        
        # 添加主要频率成分
        for i, idx in enumerate(top_indices):
            fft_features[f"fft_top_freq_{i+1}"] = frequencies[idx]
            fft_features[f"fft_top_power_{i+1}"] = magnitude[idx]
        
        return fft_features
    
    def _extract_wavelet_features(self, data):
        """提取小波特征"""
        try:
            import pywt
            
            # 使用db4小波进行4层分解
            coeffs = pywt.wavedec(data, 'db4', level=4)
            
            wavelet_features = {}
            
            # 近似系数的统计特征
            wavelet_features["wavelet_approx_mean"] = np.mean(coeffs[0])
            wavelet_features["wavelet_approx_std"] = np.std(coeffs[0])
            
            # 详细系数的统计特征
            for i, detail in enumerate(coeffs[1:], 1):
                wavelet_features[f"wavelet_detail_{i}_mean"] = np.mean(detail)
                wavelet_features[f"wavelet_detail_{i}_std"] = np.std(detail)
                wavelet_features[f"wavelet_detail_{i}_energy"] = np.sum(detail**2)
            
            return wavelet_features
            
        except ImportError:
            print("未安装pywt库，无法提取小波特征")
            return {}
        except Exception as e:
            print(f"提取小波特征错误: {e}")
            return {}
    
    def _calculate_autocorrelation(self, data, lag=1):
        """计算自相关系数"""
        import numpy as np
        
        if lag >= len(data):
            return 0
        
        # 移除均值
        data_mean = np.mean(data)
        data_centered = data - data_mean
        
        # 计算自相关
        numerator = np.sum(data_centered[:-lag] * data_centered[lag:])
        denominator = np.sum(data_centered**2)
        
        if denominator > 0:
            return numerator / denominator
        else:
            return 0
    
    def _detect_seasonality(self, data):
        """检测季节性强度"""
        # 简化的季节性检测，计算不同滞后期的自相关系数
        max_lag = min(20, len(data) // 2)
        autocorrelations = [self._calculate_autocorrelation(data, lag=i) for i in range(1, max_lag+1)]
        
        # 季节性强度可以定义为最大自相关系数
        return max(autocorrelations) if autocorrelations else 0
    
    def _calculate_sample_entropy(self, data, m=2, r=0.1):
        """计算样本熵"""
        import numpy as np
        
        n = len(data)
        if n <= m:
            return 0
        
        # 设置阈值r为数据标准差的一部分
        if isinstance(r, float):
            r = r * np.std(data)
        
        # 创建嵌入向量
        def _embed(x, d, t):
            N = len(x) - (d-1)*t
            return np.array([x[i:i+(d*t):t] for i in range(N)])
        
        # 计算相似模式的数量
        def _count_patterns(x, m, r):
            xm = _embed(x, m, 1)
            count = 0
            for i in range(len(xm)):
                # 计算与其他向量的距离
                dist = np.max(np.abs(xm - xm[i]), axis=1)
                # 统计距离小于r的向量数量
                count += np.sum(dist <= r) - 1  # 减去自己
            return count
        
        # 计算样本熵
        B = _count_patterns(data, m, r) / (n - m + 1)
        A = _count_patterns(data, m + 1, r) / (n - m)
        
        if B > 0 and A > 0:
            return -np.log(A / B)
        else:
            return 0

class TextFeatureExtractor:
    def __init__(self, model_name="sentence-transformers/all-MiniLM-L6-v2"):
        self.model_name = model_name
        self.model = None
        self.tokenizer = None
        self.feature_dim = None
        self._initialize_model()
    
    def _initialize_model(self):
        """初始化文本特征提取模型"""
        try:
            from sentence_transformers import SentenceTransformer
            
            # 加载预训练模型
            self.model = SentenceTransformer(self.model_name)
            self.feature_dim = self.model.get_sentence_embedding_dimension()
            
            print(f"文本特征提取器初始化完成，模型：{self.model_name}")
            
        except Exception as e:
            print(f"初始化文本特征提取模型错误: {e}")
    
    def extract_features(self, text, max_length=512):
        """提取文本特征"""
        try:
            if self.model is None:
                self._initialize_model()
            
            # 如果是字典，转换为文本
            if isinstance(text, dict):
                text = " ".join([f"{k}: {v}" for k, v in text.items()])
            
            # 截断过长文本
            if len(text) > max_length:
                text = text[:max_length]
            
            # 生成嵌入向量
            embedding = self.model.encode(text)
            
            # 提取基本统计特征
            import numpy as np
            stats = {
                "text_length": len(text),
                "word_count": len(text.split()),
                "embedding_mean": np.mean(embedding),
                "embedding_std": np.std(embedding),
                "embedding_norm": np.linalg.norm(embedding)
            }
            
            return {
                "embedding": embedding,
                "stats": stats,
                "feature_dim": self.feature_dim
            }
            
        except Exception as e:
            print(f"提取文本特征错误: {e}")
            return None
    
    def create_climate_description(self, data_dict):
        """为气候数据创建自然语言描述"""
        try:
            description = "气候数据分析：\n\n"
            
            # 添加卫星图像信息
            if "image_metadata" in data_dict:
                img_meta = data_dict["image_metadata"]
                description += "卫星图像信息：\n"
                if "platform" in img_meta:
                    description += f"- 卫星平台: {img_meta['platform']}\n"
                if "sensor" in img_meta:
                    description += f"- 传感器: {img_meta['sensor']}\n"
                if "acquisition_date" in img_meta:
                    description += f"- 获取时间: {img_meta['acquisition_date']}\n"
                if "resolution" in img_meta:
                    description += f"- 空间分辨率: {img_meta['resolution']}\n"
                description += "\n"
            
            # 添加时间序列数据统计
            if "time_series_stats" in data_dict:
                ts_stats = data_dict["time_series_stats"]
                description += "时间序列统计：\n"
                if "mean" in ts_stats:
                    description += f"- 平均值: {ts_stats['mean']:.2f}\n"
                if "trend" in ts_stats:
                    trend_direction = "上升" if ts_stats['trend'] > 0 else "下降"
                    description += f"- 趋势: {trend_direction} ({ts_stats['trend']:.4f}/单位时间)\n"
                if "seasonality_strength" in ts_stats:
                    seasonality = "强" if ts_stats['seasonality_strength'] > 0.5 else "弱"
                    description += f"- 季节性: {seasonality} (强度: {ts_stats['seasonality_strength']:.2f})\n"
                description += "\n"
            
            # 添加地理位置信息
            if "location_info" in data_dict:
                loc_info = data_dict["location_info"]
                description += "地理位置信息：\n"
                if "name" in loc_info:
                    description += f"- 位置名称: {loc_info['name']}\n"
                if "coordinates" in loc_info:
                    description += f"- 坐标: 经度 {loc_info['coordinates'][0]}, 纬度 {loc_info['coordinates'][1]}\n"
                if "elevation" in loc_info:
                    description += f"- 海拔: {loc_info['elevation']} 米\n"
                description += "\n"
            
            # 添加气候异常信息
            if "anomalies" in data_dict:
                anomalies = data_dict["anomalies"]
                description += "检测到的异常：\n"
                for anomaly in anomalies:
                    description += f"- {anomaly}\n"
                description += "\n"
            
            return description
            
        except Exception as e:
            print(f"创建气候描述错误: {e}")
            return "无法生成气候数据描述"

class LLMMultimodalFusion:
    def __init__(self, config=None):
        self.config = config or {
            "llm_model": "mistralai/Mistral-7B-v0.1",
            "embedding_dim": 768,
            "fusion_strategy": "cross_attention",
            "temperature": 0.7,
            "max_new_tokens": 512
        }
        
        # 初始化模态编码器
        self.image_encoder = ImageFeatureExtractor()
        self.time_series_encoder = TimeSeriesFeatureExtractor()
        self.text_encoder = TextFeatureExtractor()
        
        # 初始化LLM
        self.llm = None
        self.tokenizer = None
        self._initialize_llm()
        
        # 初始化跨模态投影层
        self.cross_modal_projection = self._build_cross_modal_projection()
        
        print("LLM多模态融合架构初始化完成")
    
    def _initialize_llm(self):
        """初始化LLM模型"""
        try:
            import torch
            from transformers import AutoModelForCausalLM, AutoTokenizer
            
            # 加载模型和分词器
            self.tokenizer = AutoTokenizer.from_pretrained(self.config["llm_model"])
            
            # 添加特殊标记用于多模态输入
            self._add_special_tokens()
            
            # 加载模型（使用量化以减少内存使用）
            self.llm = AutoModelForCausalLM.from_pretrained(
                self.config["llm_model"],
                load_in_8bit=True,
                device_map="auto",
                torch_dtype=torch.float16
            )
            
            self.llm.eval()
            print(f"LLM模型 {self.config['llm_model']} 已加载")
            
        except Exception as e:
            print(f"初始化LLM模型错误: {e}")
    
    def _add_special_tokens(self):
        """添加特殊标记用于区分不同模态"""
        special_tokens = {
            "additional_special_tokens": [
                "<image_start>", "<image_end>",
                "<time_series_start>", "<time_series_end>",
                "<text_start>", "<text_end>",
                "<fusion_start>", "<fusion_end>",
                "<task_description_start>", "<task_description_end>"
            ]
        }
        
        # 添加特殊标记
        self.tokenizer.add_special_tokens(special_tokens)
    
    def _build_cross_modal_projection(self):
        """构建跨模态投影层"""
        try:
            import torch
            import torch.nn as nn
            
            # 图像特征投影
            image_projection = nn.Sequential(
                nn.Linear(self.image_encoder.feature_dim, 1024),
                nn.LayerNorm(1024),
                nn.ReLU(),
                nn.Linear(1024, self.config["embedding_dim"])
            )
            
            # 时间序列特征投影
            # 假设时间序列特征维度为100（实际应根据提取的特征数量调整）
            ts_projection = nn.Sequential(
                nn.Linear(50, 512),  # 假设提取了50个时间序列特征
                nn.LayerNorm(512),
                nn.ReLU(),
                nn.Linear(512, self.config["embedding_dim"])
            )
            
            # 文本特征投影（如果需要）
            # 文本编码器可能已经生成了合适维度的嵌入
            text_projection = nn.Linear(self.text_encoder.feature_dim, self.config["embedding_dim"])
            
            # 将所有投影层移至适当的设备
            device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
            image_projection = image_projection.to(device)
            ts_projection = ts_projection.to(device)
            text_projection = text_projection.to(device)
            
            return {
                "image": image_projection,
                "time_series": ts_projection,
                "text": text_projection
            }
            
        except Exception as e:
            print(f"构建跨模态投影层错误: {e}")
            return None
    
    def encode_multimodal_inputs(self, image_data, time_series_data, text_data):
        """编码多模态输入数据"""
        encoded_features = {}
        
        # 编码图像数据
        if image_data is not None:
            image_features = self.image_encoder.extract_features(image_data)
            if image_features is not None:
                encoded_features["image"] = image_features["features"]
        
        # 编码时间序列数据
        if time_series_data is not None:
            ts_features = self.time_series_encoder.extract_features(time_series_data)
            if ts_features is not None:
                # 将字典转换为向量
                import numpy as np
                ts_vector = np.array([ts_features[k] for k in sorted(ts_features.keys())])
                encoded_features["time_series"] = ts_vector
        
        # 编码文本数据
        if text_data is not None:
            text_features = self.text_encoder.extract_features(text_data)
            if text_features is not None:
                encoded_features["text"] = text_features["embedding"]
        
        return encoded_features
    
    def project_to_shared_space(self, encoded_features):
        """将不同模态特征投影到共享语义空间"""
        import torch
        
        projected_features = {}
        
        # 投影图像特征
        if "image" in encoded_features and "image" in self.cross_modal_projection:
            img_tensor = torch.tensor(encoded_features["image"], dtype=torch.float32).unsqueeze(0)
            device = next(self.cross_modal_projection["image"].parameters()).device
            img_tensor = img_tensor.to(device)
            with torch.no_grad():
                projected_features["image"] = self.cross_modal_projection["image"](img_tensor).cpu().numpy()
        
        # 投影时间序列特征
        if "time_series" in encoded_features and "time_series" in self.cross_modal_projection:
            ts_tensor = torch.tensor(encoded_features["time_series"], dtype=torch.float32).unsqueeze(0)
            device = next(self.cross_modal_projection["time_series"].parameters()).device
            ts_tensor = ts_tensor.to(device)
            with torch.no_grad():
                projected_features["time_series"] = self.cross_modal_projection["time_series"](ts_tensor).cpu().numpy()
        
        # 投影文本特征（可选，因为文本编码器可能已经生成了合适的嵌入）
        if "text" in encoded_features and "text" in self.cross_modal_projection:
            text_tensor = torch.tensor(encoded_features["text"], dtype=torch.float32).unsqueeze(0)
            device = next(self.cross_modal_projection["text"].parameters()).device
            text_tensor = text_tensor.to(device)
            with torch.no_grad():
                projected_features["text"] = self.cross_modal_projection["text"](text_tensor).cpu().numpy()
        
        return projected_features
    
    def create_prompt(self, projected_features, task_description):
        """创建LLM输入提示"""
        prompt = f"<task_description_start>{task_description}<task_description_end>\n\n"
        
        # 添加各模态特征描述
        if "image" in projected_features:
            prompt += "<image_start>卫星图像特征已提取<image_end>\n"
        
        if "time_series" in projected_features:
            prompt += "<time_series_start>时间序列特征已提取<time_series_end>\n"
        
        if "text" in projected_features:
            prompt += "<text_start>文本信息已处理<text_end>\n"
        
        prompt += "\n<fusion_start>基于以上多模态信息，请提供详细的分析结果：<fusion_end>\n"
        
        return prompt
    
    def generate_analysis(self, prompt):
        """使用LLM生成分析结果"""
        try:
            import torch
            
            # 编码提示
            inputs = self.tokenizer(prompt, return_tensors="pt")
            device = next(self.llm.parameters()).device
            inputs = {k: v.to(device) for k, v in inputs.items()}
            
            # 生成响应
            with torch.no_grad():
                outputs = self.llm.generate(
                    **inputs,
                    max_new_tokens=self.config["max_new_tokens"],
                    temperature=self.config["temperature"],
                    top_p=0.9,
                    do_sample=True,
                    pad_token_id=self.tokenizer.eos_token_id
                )
            
            # 解码响应
            response = self.tokenizer.decode(outputs[0], skip_special_tokens=True)
            
            # 提取分析部分（去掉输入提示）
            analysis_start = "基于以上多模态信息，请提供详细的分析结果："
            if analysis_start in response:
                analysis = response.split(analysis_start)[1].strip()
                return analysis
            else:
                return response
                
        except Exception as e:
            print(f"生成分析结果错误: {e}")
            return "无法生成分析结果"
    
    def process_multimodal_data(self, image_data=None, time_series_data=None, 
                              text_data=None, task_description="请分析这些气候数据并提供见解。"):
        """处理多模态数据并生成分析"""
        # 1. 编码多模态输入
        encoded_features = self.encode_multimodal_inputs(image_data, time_series_data, text_data)
        
        # 2. 投影到共享空间
        projected_features = self.project_to_shared_space(encoded_features)
        
        # 3. 创建LLM提示
        prompt = self.create_prompt(projected_features, task_description)
        
        # 4. 生成分析
        analysis = self.generate_analysis(prompt)
        
        return {
            "analysis": analysis,
            "encoded_features": encoded_features,
            "projected_features": projected_features,
            "prompt": prompt
        }

class CrossModalAttention:
    def __init__(self, hidden_dim=768, num_heads=8, dropout_rate=0.1):
        self.hidden_dim = hidden_dim
        self.num_heads = num_heads
        self.dropout_rate = dropout_rate
        
        # 初始化注意力层
        self.attention = self._build_attention()
        
        print("跨模态注意力机制初始化完成")
    
    def _build_attention(self):
        """构建跨模态注意力层"""
        try:
            import torch
            import torch.nn as nn
            
            # 多头注意力
            class MultiHeadAttention(nn.Module):
                def __init__(self, hidden_dim, num_heads, dropout_rate):
                    super().__init__()
                    
                    self.hidden_dim = hidden_dim
                    self.num_heads = num_heads
                    self.head_dim = hidden_dim // num_heads
                    
                    # 确保隐藏维度可以被头数整除
                    assert self.head_dim * num_heads == hidden_dim, "hidden_dim must be divisible by num_heads"
                    
                    # 线性投影层
                    self.q_proj = nn.Linear(hidden_dim, hidden_dim)
                    self.k_proj = nn.Linear(hidden_dim, hidden_dim)
                    self.v_proj = nn.Linear(hidden_dim, hidden_dim)
                    self.out_proj = nn.Linear(hidden_dim, hidden_dim)
                    
                    # Dropout层
                    self.dropout = nn.Dropout(dropout_rate)
                    
                    # LayerNorm层
                    self.layer_norm = nn.LayerNorm(hidden_dim)
                
                def forward(self, query, key, value, attention_mask=None):
                    batch_size = query.size(0)
                    
                    # 线性投影
                    q = self.q_proj(query).view(batch_size, -1, self.num_heads, self.head_dim).transpose(1, 2)
                    k = self.k_proj(key).view(batch_size, -1, self.num_heads, self.head_dim).transpose(1, 2)
                    v = self.v_proj(value).view(batch_size, -1, self.num_heads, self.head_dim).transpose(1, 2)
                    
                    # 缩放点积注意力
                    scale = torch.sqrt(torch.tensor(self.head_dim, dtype=torch.float32))
                    attn_scores = torch.matmul(q, k.transpose(-2, -1)) / scale
                    
                    # 应用注意力掩码
                    if attention_mask is not None:
                        attn_scores = attn_scores + attention_mask
                    
                    # Softmax归一化
                    attn_weights = torch.softmax(attn_scores, dim=-1)
                    attn_weights = self.dropout(attn_weights)
                    
                    # 加权求和
                    context = torch.matmul(attn_weights, v)
                    context = context.transpose(1, 2).contiguous().view(batch_size, -1, self.hidden_dim)
                    
                    # 输出投影
                    output = self.out_proj(context)
                    output = self.layer_norm(output + query)  # 残差连接
                    
                    return output, attn_weights
            
            # 创建注意力层
            attention = MultiHeadAttention(
                hidden_dim=self.hidden_dim,
                num_heads=self.num_heads,
                dropout_rate=self.dropout_rate
            )
            
            # 移至适当的设备
            device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
            attention = attention.to(device)
            
            return attention
            
        except Exception as e:
            print(f"构建注意力层错误: {e}")
            return None
    
    def create_attention_mask(self, query_length, key_length):
        """创建注意力掩码"""
        import torch
        
        # 这里简单返回None，实际应用中可能需要更复杂的掩码
        return None
    
    def fuse_modalities(self, modalities):
        """融合多种模态的特征"""
        import torch
        
        if not modalities:
            print("没有模态数据可供融合")
            return None
        
        try:
            # 将模态特征转换为张量
            modality_tensors = []
            modality_keys = []
            
            for key, feature in modalities.items():
                if feature is not None:
                    # 确保特征是二维的 (batch_size, hidden_dim)
                    if len(feature.shape) == 1:
                        tensor = torch.tensor(feature, dtype=torch.float32).unsqueeze(0).unsqueeze(0)
                    elif len(feature.shape) == 2:
                        tensor = torch.tensor(feature, dtype=torch.float32).unsqueeze(0)
                    else:
                        tensor = torch.tensor(feature, dtype=torch.float32)
                    
                    modality_tensors.append(tensor)
                    modality_keys.append(key)
            
            if not modality_tensors:
                print("没有有效的模态特征可供融合")
                return None
            
            # 移至设备
            device = next(self.attention.parameters()).device
            modality_tensors = [tensor.to(device) for tensor in modality_tensors]
            
            # 初始化融合特征
            fused_features = modality_tensors[0]
            
            # 逐模态融合
            for i in range(1, len(modality_tensors)):
                # 使用多头注意力融合当前模态与累积特征
                mask = self.create_attention_mask(
                    fused_features.size(1),
                    modality_tensors[i].size(1)
                )
                
                # 交叉注意力
                with torch.no_grad():
                    output, _ = self.attention(
                        query=fused_features,
                        key=modality_tensors[i],
                        value=modality_tensors[i],
                        attention_mask=mask
                    )
                
                # 更新融合特征
                # 可以使用不同的融合策略
                fused_features = torch.cat([fused_features, output], dim=1)
            
            # 计算自注意力，整合所有模态信息
            with torch.no_grad():
                final_output, attention_weights = self.attention(
                    query=fused_features,
                    key=fused_features,
                    value=fused_features
                )
            
            return {
                "fused_features": final_output.cpu().numpy(),
                "attention_weights": attention_weights.cpu().numpy(),
                "modalities_used": modality_keys
            }
            
        except Exception as e:
            print(f"融合模态特征错误: {e}")
            return None
    
    def visualize_attention(self, attention_weights, modalities):
        """可视化注意力权重"""
        try:
            import matplotlib.pyplot as plt
            import numpy as np
            
            # 简化的可视化示例
            plt.figure(figsize=(10, 8))
            
            # 取第一个头的注意力权重作为示例
            if len(attention_weights.shape) >= 4:
                attention_to_plot = attention_weights[0, 0, :, :]
            else:
                attention_to_plot = attention_weights[0, :, :]
            
            # 创建热力图
            plt.imshow(attention_to_plot, cmap='viridis')
            plt.colorbar(label='注意力权重')
            plt.title('跨模态注意力可视化')
            
            # 设置坐标轴标签
            if modalities:
                plt.xlabel('键模态')
                plt.ylabel('查询模态')
            
            # 保存图像
            plt.savefig('cross_modal_attention.png')
            plt.close()
            
            print("注意力可视化已保存到 cross_modal_attention.png")
            return True
            
        except Exception as e:
            print(f"可视化注意力权重错误: {e}")
            return False

class SpatioTemporalModel:
    def __init__(self, config=None):
        self.config = config or {
            "spatial_encoding_dim": 128,
            "temporal_encoding_dim": 128,
            "num_spatial_layers": 2,
            "num_temporal_layers": 2,
            "dropout_rate": 0.1
        }
        
        # 初始化空间和时间编码器
        self.spatial_encoder = self._build_spatial_encoder()
        self.temporal_encoder = self._build_temporal_encoder()
        self.spatio_temporal_fusion = self._build_fusion_module()
        
        print("时空关联建模模块初始化完成")
    
    def _build_spatial_encoder(self):
        """构建空间编码器"""
        try:
            import torch
            import torch.nn as nn
            
            class SpatialEncoder(nn.Module):
                def __init__(self, input_dim, encoding_dim, num_layers, dropout_rate):
                    super().__init__()
                    
                    # 构建多层感知机
                    layers = []
                    prev_dim = input_dim
                    
                    for _ in range(num_layers - 1):
                        layers.extend([
                            nn.Linear(prev_dim, prev_dim // 2),
                            nn.LayerNorm(prev_dim // 2),
                            nn.ReLU(),
                            nn.Dropout(dropout_rate)
                        ])
                        prev_dim = prev_dim // 2
                    
                    # 输出层
                    layers.append(nn.Linear(prev_dim, encoding_dim))
                    
                    self.model = nn.Sequential(*layers)
                
                def forward(self, x):
                    return self.model(x)
            
            # 创建空间编码器
            spatial_encoder = SpatialEncoder(
                input_dim=2,  # 经度和纬度
                encoding_dim=self.config["spatial_encoding_dim"],
                num_layers=self.config["num_spatial_layers"],
                dropout_rate=self.config["dropout_rate"]
            )
            
            # 移至适当的设备
            device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
            spatial_encoder = spatial_encoder.to(device)
            
            return spatial_encoder
            
        except Exception as e:
            print(f"构建空间编码器错误: {e}")
            return None
    
    def _build_temporal_encoder(self):
        """构建时间编码器"""
        try:
            import torch
            import torch.nn as nn
            
            class TemporalEncoder(nn.Module):
                def __init__(self, input_dim, encoding_dim, num_layers, dropout_rate):
                    super().__init__()
                    
                    # 构建时间编码器
                    self.rnn = nn.GRU(
                        input_size=input_dim,
                        hidden_size=encoding_dim,
                        num_layers=num_layers,
                        batch_first=True,
                        dropout=dropout_rate if num_layers > 1 else 0,
                        bidirectional=True
                    )
                    
                    # 最终投影层
                    self.projection = nn.Linear(encoding_dim * 2, encoding_dim)
                    self.layer_norm = nn.LayerNorm(encoding_dim)
                
                def forward(self, x):
                    # 处理变长序列
                    batch_size, seq_len, _ = x.size()
                    
                    # 应用GRU
                    outputs, _ = self.rnn(x)
                    
                    # 取最后时间步的输出
                    last_output = outputs[:, -1, :]
                    
                    # 投影到目标维度
                    output = self.projection(last_output)
                    output = self.layer_norm(output)
                    
                    return output
            
            # 创建时间编码器
            temporal_encoder = TemporalEncoder(
                input_dim=1,  # 假设每个时间步只有一个特征
                encoding_dim=self.config["temporal_encoding_dim"],
                num_layers=self.config["num_temporal_layers"],
                dropout_rate=self.config["dropout_rate"]
            )
            
            # 移至适当的设备
            device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
            temporal_encoder = temporal_encoder.to(device)
            
            return temporal_encoder
            
        except Exception as e:
            print(f"构建时间编码器错误: {e}")
            return None
    
    def _build_fusion_module(self):
        """构建时空融合模块"""
        try:
            import torch
            import torch.nn as nn
            
            class SpatioTemporalFusion(nn.Module):
                def __init__(self, spatial_dim, temporal_dim, hidden_dim):
                    super().__init__()
                    
                    # 融合层
                    self.fusion = nn.Sequential(
                        nn.Linear(spatial_dim + temporal_dim, hidden_dim),
                        nn.LayerNorm(hidden_dim),
                        nn.ReLU(),
                        nn.Linear(hidden_dim, hidden_dim)
                    )
                    
                    # 门控机制
                    self.gate = nn.Sequential(
                        nn.Linear(spatial_dim + temporal_dim, hidden_dim),
                        nn.Sigmoid()
                    )
                
                def forward(self, spatial_features, temporal_features):
                    # 特征拼接
                    combined = torch.cat([spatial_features, temporal_features], dim=-1)
                    
                    # 融合特征
                    fused = self.fusion(combined)
                    
                    # 应用门控
                    gate = self.gate(combined)
                    gated_output = fused * gate
                    
                    return gated_output
            
            # 创建时空融合模块
            hidden_dim = self.config["spatial_encoding_dim"] + self.config["temporal_encoding_dim"]
            fusion_module = SpatioTemporalFusion(
                spatial_dim=self.config["spatial_encoding_dim"],
                temporal_dim=self.config["temporal_encoding_dim"],
                hidden_dim=hidden_dim
            )
            
            # 移至适当的设备
            device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
            fusion_module = fusion_module.to(device)
            
            return fusion_module
            
        except Exception as e:
            print(f"构建时空融合模块错误: {e}")
            return None
    
    def encode_spatial(self, coordinates):
        """编码空间坐标"""
        import torch
        
        try:
            # 准备输入
            coords_tensor = torch.tensor(coordinates, dtype=torch.float32).unsqueeze(0)
            
            # 移至设备
            device = next(self.spatial_encoder.parameters()).device
            coords_tensor = coords_tensor.to(device)
            
            # 编码
            with torch.no_grad():
                spatial_encoding = self.spatial_encoder(coords_tensor)
            
            return spatial_encoding.cpu().numpy()
            
        except Exception as e:
            print(f"编码空间坐标错误: {e}")
            return None
    
    def encode_temporal(self, time_series):
        """编码时间序列"""
        import torch
        import numpy as np
        
        try:
            # 准备输入
            # 移除NaN值
            valid_ts = time_series[~np.isnan(time_series)]
            
            if len(valid_ts) < 3:
                print("有效时间序列数据不足")
                return None
            
            # 重塑为(batch_size, seq_len, input_dim)
            ts_tensor = torch.tensor(valid_ts, dtype=torch.float32).unsqueeze(0).unsqueeze(-1)
            
            # 移至设备
            device = next(self.temporal_encoder.parameters()).device
            ts_tensor = ts_tensor.to(device)
            
            # 编码
            with torch.no_grad():
                temporal_encoding = self.temporal_encoder(ts_tensor)
            
            return temporal_encoding.cpu().numpy()
            
        except Exception as e:
            print(f"编码时间序列错误: {e}")
            return None
    
    def fuse_spatio_temporal(self, spatial_encoding, temporal_encoding):
        """融合时空编码"""
        import torch
        
        try:
            # 准备输入
            spatial_tensor = torch.tensor(spatial_encoding, dtype=torch.float32)
            temporal_tensor = torch.tensor(temporal_encoding, dtype=torch.float32)
            
            # 移至设备
            device = next(self.spatio_temporal_fusion.parameters()).device
            spatial_tensor = spatial_tensor.to(device)
            temporal_tensor = temporal_tensor.to(device)
            
            # 融合
            with torch.no_grad():
                fused = self.spatio_temporal_fusion(spatial_tensor, temporal_tensor)
            
            return fused.cpu().numpy()
            
        except Exception as e:
            print(f"融合时空编码错误: {e}")
            return None
    
    def detect_spatio_temporal_patterns(self, data_points):
        """检测时空模式"""
        import numpy as np
        
        try:
            if not data_points:
                print("没有数据点可供分析")
                return None
            
            # 提取坐标和时间序列
            coordinates = np.array([point["coordinates"] for point in data_points])
            time_series_list = [point["time_series"] for point in data_points]
            
            # 对每个数据点进行时空编码和融合
            fused_features = []
            valid_indices = []
            
            for i, (coords, ts) in enumerate(zip(coordinates, time_series_list)):
                # 编码空间信息
                spatial_enc = self.encode_spatial(coords)
                if spatial_enc is None:
                    continue
                
                # 编码时间信息
                temporal_enc = self.encode_temporal(ts)
                if temporal_enc is None:
                    continue
                
                # 融合时空特征
                fused = self.fuse_spatio_temporal(spatial_enc, temporal_enc)
                if fused is not None:
                    fused_features.append(fused.flatten())
                    valid_indices.append(i)
            
            if not fused_features:
                print("无法生成有效的时空特征")
                return None
            
            # 转换为数组
            fused_features = np.array(fused_features)
            
            # 简单的聚类分析，检测时空模式
            from sklearn.cluster import DBSCAN
            
            # 使用DBSCAN聚类
            clustering = DBSCAN(eps=0.3, min_samples=2).fit(fused_features)
            labels = clustering.labels_
            
            # 分析聚类结果
            patterns = {}
            for i, label in enumerate(labels):
                if label not in patterns:
                    patterns[label] = []
                original_index = valid_indices[i]
                patterns[label].append({
                    "data_point": data_points[original_index],
                    "feature_vector": fused_features[i]
                })
            
            # 识别异常点（标签为-1的点）
            anomalies = [p for i, p in enumerate(valid_indices) if labels[i] == -1]
            
            return {
                "patterns": patterns,
                "anomalies": anomalies,
                "total_patterns": len(patterns),
                "total_anomalies": len(anomalies)
            }
            
        except Exception as e:
            print(f"检测时空模式错误: {e}")
            return None

class SatelliteImageAnalyzer:
    def __init__(self, config=None):
        self.config = config or {
            "image_size": (224, 224),
            "feature_extractor_model": "resnet50",
            "segmentation_model": "deeplabv3_resnet101",
            "num_classes": 19,  # ESA WorldCover 分类数
            "device": "cuda" if torch.cuda.is_available() else "cpu"
        }
        
        # 初始化模型
        self.feature_extractor = None
        self.segmentation_model = None
        self.image_preprocessor = None
        self._initialize_models()
        
        # 初始化分析工具
        self.climate_phenomena_detector = ClimatePhenomenaDetector()
        
        print("卫星图像分析器初始化完成")
    
    def _initialize_models(self):
        """初始化图像分析模型"""
        try:
            import torch
            import torchvision.transforms as transforms
            from torchvision import models
            
            # 设置设备
            device = torch.device(self.config["device"])
            
            # 初始化预处理器
            self.image_preprocessor = transforms.Compose([
                transforms.Resize(self.config["image_size"]),
                transforms.ToTensor(),
                transforms.Normalize(mean=[0.485, 0.456, 0.406],
                                     std=[0.229, 0.224, 0.225])
            ])
            
            # 初始化特征提取器
            if self.config["feature_extractor_model"] == "resnet50":
                self.feature_extractor = models.resnet50(pretrained=True)
                # 移除最后一层全连接层
                self.feature_extractor = torch.nn.Sequential(*list(self.feature_extractor.children())[:-1])
                self.feature_extractor.eval()
                self.feature_extractor = self.feature_extractor.to(device)
            
            # 初始化分割模型
            if self.config["segmentation_model"] == "deeplabv3_resnet101":
                self.segmentation_model = models.segmentation.deeplabv3_resnet101(pretrained=True)
                self.segmentation_model.eval()
                self.segmentation_model = self.segmentation_model.to(device)
            
            print("图像分析模型初始化完成")
            
        except Exception as e:
            print(f"初始化图像分析模型错误: {e}")
    
    def load_image(self, image_path):
        """加载并预处理图像"""
        try:
            from PIL import Image
            import numpy as np
            
            # 加载图像
            image = Image.open(image_path).convert('RGB')
            original_size = image.size
            
            # 预处理
            preprocessed_image = self.image_preprocessor(image)
            
            return {
                "original": np.array(image),
                "preprocessed": preprocessed_image,
                "original_size": original_size
            }
            
        except Exception as e:
            print(f"加载图像错误: {e}")
            return None
    
    def extract_features(self, image_data):
        """从图像中提取特征"""
        try:
            import torch
            
            if self.feature_extractor is None:
                print("特征提取器未初始化")
                return None
            
            # 确保图像数据是正确的格式
            if isinstance(image_data, dict):
                image_tensor = image_data["preprocessed"]
            else:
                image_tensor = image_data
            
            # 添加批次维度
            if len(image_tensor.shape) == 3:
                image_tensor = image_tensor.unsqueeze(0)
            
            # 移至设备
            device = torch.device(self.config["device"])
            image_tensor = image_tensor.to(device)
            
            # 提取特征
            with torch.no_grad():
                features = self.feature_extractor(image_tensor)
            
            # 展平特征向量
            features = features.squeeze().cpu().numpy()
            
            return {
                "features": features,
                "shape": features.shape
            }
            
        except Exception as e:
            print(f"提取图像特征错误: {e}")
            return None
    
    def segment_image(self, image_data):
        """对图像进行语义分割"""
        try:
            import torch
            import numpy as np
            
            if self.segmentation_model is None:
                print("分割模型未初始化")
                return None
            
            # 确保图像数据是正确的格式
            if isinstance(image_data, dict):
                image_tensor = image_data["preprocessed"]
                original_size = image_data["original_size"]
            else:
                image_tensor = image_data
                original_size = self.config["image_size"]
            
            # 添加批次维度
            if len(image_tensor.shape) == 3:
                image_tensor = image_tensor.unsqueeze(0)
            
            # 移至设备
            device = torch.device(self.config["device"])
            image_tensor = image_tensor.to(device)
            
            # 进行分割
            with torch.no_grad():
                output = self.segmentation_model(image_tensor)['out']
            
            # 获取预测类别
            pred = output.argmax(1).squeeze().cpu().numpy()
            
            # 调整回原始大小
            from PIL import Image
            pred_img = Image.fromarray(pred.astype(np.uint8))
            pred_resized = pred_img.resize(original_size, Image.NEAREST)
            pred_resized = np.array(pred_resized)
            
            return {
                "segmentation": pred_resized,
                "class_distribution": np.bincount(pred_resized.flatten())
            }
            
        except Exception as e:
            print(f"图像分割错误: {e}")
            return None
    
    def detect_climate_phenomena(self, image_data, time_series=None):
        """检测气候现象"""
        try:
            # 提取图像特征
            features = self.extract_features(image_data)
            if features is None:
                return None
            
            # 图像分割
            segmentation = self.segment_image(image_data)
            if segmentation is None:
                return None
            
            # 使用专门的检测器分析气候现象
            phenomena = self.climate_phenomena_detector.detect(
                features["features"],
                segmentation["segmentation"],
                time_series
            )
            
            return phenomena
            
        except Exception as e:
            print(f"检测气候现象错误: {e}")
            return None
    
    def generate_climate_description(self, image_path, time_series=None):
        """生成图像内容的气候描述"""
        try:
            # 加载图像
            image_data = self.load_image(image_path)
            if image_data is None:
                return "无法加载图像"
            
            # 检测气候现象
            phenomena = self.detect_climate_phenomena(image_data, time_series)
            
            # 生成描述
            description_parts = []
            description_parts.append("卫星图像分析结果：")
            
            # 添加分割结果描述
            if phenomena and "land_cover" in phenomena:
                land_cover = phenomena["land_cover"]
                description_parts.append(f"\n地表覆盖类型：")
                for cover_type, percentage in land_cover.items():
                    description_parts.append(f"- {cover_type}: {percentage:.1f}%")
            
            # 添加检测到的气候现象
            if phenomena and "detected_phenomena" in phenomena:
                detected = phenomena["detected_phenomena"]
                description_parts.append(f"\n检测到的气候现象：")
                for phenomenon, confidence in detected.items():
                    description_parts.append(f"- {phenomenon}: 置信度 {confidence:.2f}")
            
            # 添加异常分析
            if phenomena and "anomalies" in phenomena:
                anomalies = phenomena["anomalies"]
                if anomalies:
                    description_parts.append(f"\n异常现象：")
                    for anomaly in anomalies:
                        description_parts.append(f"- {anomaly}")
                else:
                    description_parts.append("\n未检测到明显异常现象。")
            
            # 合并描述
            return "\n".join(description_parts)
            
        except Exception as e:
            print(f"生成气候描述错误: {e}")
            return "无法生成图像分析描述"

class ClimatePhenomenaDetector:
    def __init__(self):
        # 定义常见气候现象的特征和分类阈值
        self.phenomena_patterns = {
            "deforestation": {
                "min_change": 0.15,  # 最小变化百分比
                "reference_class": 4,  # 森林类别
                "target_class": 1,    # 裸地类别
                "confidence_threshold": 0.7
            },
            "urban_expansion": {
                "min_change": 0.05,
                "reference_classes": [0, 2],  # 水体和植被
                "target_class": 7,            # 建筑类别
                "confidence_threshold": 0.6
            },
            "glacier_melt": {
                "min_change": 0.2,
                "reference_class": 15,
                "target_class": 0,
                "confidence_threshold": 0.65
            },
            "flooding": {
                "min_change": 0.1,
                "reference_classes": [1, 2, 4],  # 裸地、植被、森林
                "target_class": 0,               # 水体
                "confidence_threshold": 0.75
            },
            "drought": {
                "min_change": 0.15,
                "reference_class": 2,  # 植被
                "target_classes": [1, 8],  # 裸地、旱地
                "confidence_threshold": 0.6
            }
        }
        
        # 土地覆盖类别映射
        self.land_cover_classes = {
            0: "水体",
            1: "裸地",
            2: "草本植被",
            3: "灌木",
            4: "森林",
            5: "苔原",
            6: "封闭灌丛",
            7: "建筑区",
            8: "旱地",
            9: "永久性雪/冰",
            10: "草地",
            11: "湿地",
            12: "农作物",
            13: "红树林",
            14: "苔藓",
            15: "冰川",
            16: "沙漠",
            17: "河口",
            18: "珊瑚礁"
        }
    
    def detect(self, image_features, segmentation_map, time_series=None):
        """检测气候相关现象"""
        import numpy as np
        
        results = {
            "land_cover": {},
            "detected_phenomena": {},
            "anomalies": []
        }
        
        try:
            # 计算土地覆盖分布
            unique, counts = np.unique(segmentation_map, return_counts=True)
            total_pixels = np.sum(counts)
            
            for u, c in zip(unique, counts):
                if u in self.land_cover_classes:
                    class_name = self.land_cover_classes[u]
                    percentage = (c / total_pixels) * 100
                    results["land_cover"][class_name] = percentage
            
            # 检测气候现象
            for phenomenon, pattern in self.phenomena_patterns.items():
                confidence = self._detect_phenomenon(segmentation_map, pattern)
                if confidence > pattern["confidence_threshold"]:
                    results["detected_phenomena"][phenomenon] = confidence
            
            # 如果有时间序列数据，进行更深入的分析
            if time_series is not None:
                time_anomalies = self._analyze_time_series(time_series)
                if time_anomalies:
                    results["anomalies"].extend(time_anomalies)
            
            # 检查异常情况
            self._check_for_anomalies(results, segmentation_map)
            
            return results
            
        except Exception as e:
            print(f"检测气候现象错误: {e}")
            return results
    
    def _detect_phenomenon(self, segmentation_map, pattern):
        """检测特定气候现象"""
        import numpy as np
        
        try:
            # 计算参考类别和目标类别的像素数
            unique, counts = np.unique(segmentation_map, return_counts=True)
            pixel_counts = dict(zip(unique, counts))
            total_pixels = np.sum(counts)
            
            # 检查参考类别
            reference_pixels = 0
            if "reference_class" in pattern:
                if pattern["reference_class"] in pixel_counts:
                    reference_pixels = pixel_counts[pattern["reference_class"]]
            elif "reference_classes" in pattern:
                for ref_class in pattern["reference_classes"]:
                    if ref_class in pixel_counts:
                        reference_pixels += pixel_counts[ref_class]
            
            # 检查目标类别
            target_pixels = 0
            if "target_class" in pattern:
                if pattern["target_class"] in pixel_counts:
                    target_pixels = pixel_counts[pattern["target_class"]]
            elif "target_classes" in pattern:
                for target_class in pattern["target_classes"]:
                    if target_class in pixel_counts:
                        target_pixels += pixel_counts[target_class]
            
            # 计算置信度（简化示例）
            # 实际应用中可能需要更复杂的逻辑，如纹理分析、空间分布等
            if reference_pixels == 0:
                return 0.0
            
            # 基于类别比例计算置信度
            confidence = min(1.0, (target_pixels / total_pixels) * 10)
            
            # 如果有最小变化要求，应用该要求
            if "min_change" in pattern:
                if target_pixels / total_pixels < pattern["min_change"]:
                    confidence = confidence * 0.5  # 降低置信度
            
            return confidence
            
        except Exception as e:
            print(f"检测特定气候现象错误: {e}")
            return 0.0
    
    def _analyze_time_series(self, time_series):
        """分析时间序列数据以检测异常"""
        import numpy as np
        from scipy import stats
        
        anomalies = []
        
        try:
            # 移除NaN值
            valid_data = time_series[~np.isnan(time_series)]
            
            if len(valid_data) < 10:
                return anomalies
            
            # 使用Z-score检测异常值
            z_scores = np.abs(stats.zscore(valid_data))
            threshold = 3.0  # Z-score阈值
            
            # 查找异常值索引
            anomaly_indices = np.where(z_scores > threshold)[0]
            
            # 如果发现异常值，添加到结果中
            if len(anomaly_indices) > 0:
                anomaly_percentage = (len(anomaly_indices) / len(valid_data)) * 100
                anomalies.append(f"检测到{anomaly_percentage:.1f}%的时间序列异常值")
                
                # 检查趋势
                x = np.arange(len(valid_data))
                slope, _, r_value, _, _ = stats.linregress(x, valid_data)
                
                # 如果趋势显著
                if abs(r_value) > 0.7:
                    trend_direction = "上升" if slope > 0 else "下降"
                    anomalies.append(f"观察到显著的{trend_direction}趋势")
            
            return anomalies
            
        except Exception as e:
            print(f"分析时间序列错误: {e}")
            return anomalies
    
    def _check_for_anomalies(self, results, segmentation_map):
        """检查土地覆盖异常情况"""
        try:
            # 检查极端情况
            land_cover = results["land_cover"]
            
            # 检查森林覆盖异常低
            if "森林" in land_cover and land_cover["森林"] < 5:
                results["anomalies"].append("区域森林覆盖率异常低")
            
            # 检查建筑区异常高
            if "建筑区" in land_cover and land_cover["建筑区"] > 50:
                results["anomalies"].append("区域建筑覆盖率异常高，可能存在过度城市化")
            
            # 检查水体异常变化
            if "水体" in land_cover:
                if land_cover["水体"] < 1:
                    results["anomalies"].append("区域水体覆盖率极低")
                elif land_cover["水体"] > 70:
                    results["anomalies"].append("区域水体覆盖率异常高，可能存在大规模洪水")
            
            # 检查裸地和沙漠比例
            bare_percentage = land_cover.get("裸地", 0) + land_cover.get("沙漠", 0)
            if bare_percentage > 60:
                results["anomalies"].append("区域裸地/沙漠覆盖率异常高，可能存在严重荒漠化")
            
        except Exception as e:
            print(f"检查异常情况错误: {e}")

class ClimateEventDetectionSystem:
    def __init__(self, config=None):
        self.config = config or {
            "detection_threshold": 0.6,
            "confidence_threshold": 0.7,
            "min_event_duration": 3,  # 最小事件持续天数
            "max_event_gap": 2,       # 允许的最大事件间隔天数
            "time_window_size": 30    # 时间窗口大小（天）
        }
        
        # 初始化必要的模块
        self.satellite_analyzer = None
        self.temporal_analyzer = None
        self._initialize_modules()
        
        # 事件类型定义
        self.event_types = {
            "flood": {
                "name": "洪水",
                "description": "由暴雨、融雪或堤坝决口引起的水位异常上升",
                "indicators": ["water_extent_increase", "precipitation_spike", "river_discharge_rise"]
            },
            "drought": {
                "name": "干旱",
                "description": "长期缺乏降水导致的水分短缺",
                "indicators": ["ndvi_decrease", "temperature_rise", "precipitation_decrease", "soil_moisture_decline"]
            },
            "wildfire": {
                "name": "野火",
                "description": "自然或人为原因引起的大规模火灾",
                "indicators": ["thermal_anomaly", "smoke_plume", "vegetation_loss", "temperature_spike"]
            },
            "hurricane": {
                "name": "飓风/台风",
                "description": "强烈的热带气旋，伴有强风和暴雨",
                "indicators": ["wind_speed_increase", "pressure_drop", "storm_surge", "heavy_precipitation"]
            },
            "landslide": {
                "name": "山体滑坡",
                "description": "由重力引起的大量土壤和岩石下滑",
                "indicators": ["slope_instability", "precipitation_trigger", "terrain_change", "vegetation_disturbance"]
            },
            "glacier_melt": {
                "name": "冰川融化",
                "description": "冰川质量因温度升高而减少",
                "indicators": ["ice_extent_decrease", "temperature_rise", "albedo_change", "melt_water_increase"]
            },
            "sea_level_rise": {
                "name": "海平面上升",
                "description": "全球或区域海平面的长期上升趋势",
                "indicators": ["coastal_erosion", "sea_water_level_increase", "saltwater_intrusion"]
            },
            "extreme_heat": {
                "name": "极端高温",
                "description": "异常高的温度持续一段时间",
                "indicators": ["temperature_extreme", "heat_index_spike", "humidity_change"]
            },
            "extreme_cold": {
                "name": "极端低温",
                "description": "异常低的温度持续一段时间",
                "indicators": ["temperature_extreme", "frost_formation", "precipitation_type_change"]
            },
            "tornado": {
                "name": "龙卷风",
                "description": "快速旋转的空气柱，接触地面和云层",
                "indicators": ["wind_shear", "thunderstorm_formation", "funnel_cloud", "damage_pattern"]
            }
        }
        
        print("气候事件检测系统初始化完成")
    
    def _initialize_modules(self):
        """初始化分析模块"""
        try:
            # 初始化卫星图像分析器
            self.satellite_analyzer = SatelliteImageAnalyzer()
            
            # 初始化时间序列分析器
            self.temporal_analyzer = TemporalAnalyzer()
            
            print("分析模块初始化完成")
            
        except Exception as e:
            print(f"初始化分析模块错误: {e}")
    
    def extract_indicators(self, satellite_data, time_series_data):
        """从多源数据中提取气候指标"""
        indicators = {}
        
        try:
            # 从卫星数据中提取指标
            if satellite_data is not None:
                # 图像分析
                if isinstance(satellite_data, str):  # 假设是文件路径
                    image_description = self.satellite_analyzer.generate_climate_description(
                        satellite_data, time_series_data
                    )
                    indicators["satellite_analysis"] = image_description
                
                # 提取水体覆盖变化
                water_extent = self._estimate_water_extent(satellite_data)
                indicators["water_extent"] = water_extent
                
                # 提取植被健康状况
                vegetation_health = self._assess_vegetation_health(satellite_data)
                indicators["vegetation_health"] = vegetation_health
            
            # 从时间序列数据中提取指标
            if time_series_data is not None:
                # 温度趋势
                temp_trend = self.temporal_analyzer.analyze_trend(
                    time_series_data, "temperature"
                )
                indicators["temperature_trend"] = temp_trend
                
                # 降水分析
                precip_analysis = self.temporal_analyzer.analyze_precipitation(
                    time_series_data
                )
                indicators["precipitation"] = precip_analysis
                
                # 检测异常
                anomalies = self.temporal_analyzer.detect_anomalies(
                    time_series_data
                )
                indicators["anomalies"] = anomalies
            
            return indicators
            
        except Exception as e:
            print(f"提取气候指标错误: {e}")
            return indicators
    
    def detect_events(self, indicators, location=None, timestamp=None):
        """基于指标检测气候事件"""
        detected_events = []
        
        try:
            for event_type, event_info in self.event_types.items():
                # 计算事件检测分数
                detection_score = self._calculate_detection_score(
                    indicators, event_info["indicators"]
                )
                
                # 如果检测分数超过阈值，将事件添加到结果中
                if detection_score >= self.config["detection_threshold"]:
                    # 生成详细的事件报告
                    event_details = self._generate_event_report(
                        event_type, event_info, indicators, detection_score
                    )
                    
                    # 添加位置和时间信息（如果提供）
                    if location is not None:
                        event_details["location"] = location
                    if timestamp is not None:
                        event_details["timestamp"] = timestamp
                    
                    detected_events.append(event_details)
            
            # 根据置信度排序
            detected_events.sort(key=lambda x: x["confidence"], reverse=True)
            
            return detected_events
            
        except Exception as e:
            print(f"检测气候事件错误: {e}")
            return detected_events
    
    def _calculate_detection_score(self, indicators, required_indicators):
        """计算事件检测分数"""
        import numpy as np
        
        try:
            # 简化的检测分数计算
            # 实际应用中可能需要更复杂的算法和权重
            score = 0.0
            matched_indicators = 0
            
            for req_ind in required_indicators:
                # 检查指标是否存在并满足条件
                # 这里使用简化的匹配逻辑，实际应用需要更详细的判断
                if self._match_indicator(indicators, req_ind):
                    matched_indicators += 1
            
            # 计算分数（匹配指标数/总指标数）
            if required_indicators:
                score = matched_indicators / len(required_indicators)
            
            return score
            
        except Exception as e:
            print(f"计算检测分数错误: {e}")
            return 0.0
    
    def _match_indicator(self, indicators, required_indicator):
        """检查是否匹配所需指标"""
        try:
            # 简化的指标匹配逻辑
            # 实际应用中需要更详细的匹配规则
            if "water_extent" in required_indicator and "water_extent" in indicators:
                # 检查水体覆盖是否显著增加
                water_extent = indicators["water_extent"]
                return water_extent.get("change", 0) > 0.1  # 如果增加超过10%
            
            elif "temperature" in required_indicator and "temperature_trend" in indicators:
                # 检查温度趋势
                temp_trend = indicators["temperature_trend"]
                if "rise" in required_indicator:
                    return temp_trend.get("direction", "") == "increasing"
                elif "extreme" in required_indicator:
                    return temp_trend.get("is_extreme", False)
            
            elif "precipitation" in required_indicator and "precipitation" in indicators:
                # 检查降水情况
                precip = indicators["precipitation"]
                if "spike" in required_indicator or "heavy" in required_indicator:
                    return precip.get("is_heavy", False)
                elif "decrease" in required_indicator:
                    return precip.get("trend", "") == "decreasing"
            
            elif "vegetation" in required_indicator and "vegetation_health" in indicators:
                # 检查植被健康
                veg_health = indicators["vegetation_health"]
                if "loss" in required_indicator or "decline" in required_indicator:
                    return veg_health.get("health_score", 1.0) < 0.5
            
            # 检查是否有异常检测
            if "anomaly" in required_indicator.lower() and "anomalies" in indicators:
                anomalies = indicators["anomalies"]
                return len(anomalies) > 0
            
            # 默认不匹配
            return False
            
        except Exception as e:
            print(f"匹配指标错误: {e}")
            return False
    
    def _generate_event_report(self, event_type, event_info, indicators, detection_score):
        """生成气候事件详细报告"""
        try:
            # 基础报告
            report = {
                "event_type": event_type,
                "event_name": event_info["name"],
                "description": event_info["description"],
                "detection_score": detection_score,
                "confidence": min(1.0, detection_score * 1.2),  # 简单的置信度转换
                "key_indicators": [],
                "severity": "unknown",
                "potential_impacts": [],
                "recommended_actions": []
            }
            
            # 确定事件严重性（简化示例）
            if detection_score > 0.85:
                report["severity"] = "高"
            elif detection_score > 0.7:
                report["severity"] = "中"
            else:
                report["severity"] = "低"
            
            # 识别关键指标
            for indicator in event_info["indicators"]:
                if self._match_indicator(indicators, indicator):
                    report["key_indicators"].append(indicator)
            
            # 添加潜在影响
            report["potential_impacts"] = self._estimate_impacts(event_type, report["severity"])
            
            # 添加建议行动
            report["recommended_actions"] = self._suggest_actions(
                event_type, report["severity"]
            )
            
            return report
            
        except Exception as e:
            print(f"生成事件报告错误: {e}")
            return {
                "event_type": event_type,
                "error": str(e)
            }
    
    def _estimate_impacts(self, event_type, severity):
        """估计气候事件的潜在影响"""
        impacts = {
            "flood": {
                "high": ["基础设施严重损坏", "人员伤亡风险高", "大范围农作物损失", "传染病风险增加"],
                "medium": ["局部地区淹没", "交通中断", "财产损失", "电力供应不稳定"],
                "low": ["轻微积水", "排水系统压力增加"]
            },
            "drought": {
                "high": ["严重水资源短缺", "农作物大幅减产", "牲畜死亡", "野火风险增加"],
                "medium": ["水库水位下降", "灌溉限制", "地下水过度开采"],
                "low": ["植被生长减缓", "轻微用水限制"]
            },
            "wildfire": {
                "high": ["大面积植被破坏", "建筑物损失", "人员伤亡风险", "空气质量严重恶化"],
                "medium": ["局部地区火灾", "烟雾影响", "野生动物栖息地破坏"],
                "low": ["小型可控火灾", "有限的生态系统影响"]
            },
            "hurricane": {
                "high": ["灾难性风力损害", "风暴潮淹没", "广泛停电", "大规模撤离"],
                "medium": ["建筑物损坏", "洪水", "基础设施中断"],
                "low": ["强风", "降雨", "小型洪水"]
            }
        }
        
        # 默认影响
        default_impacts = ["需要进一步监测", "可能影响当地生活"]
        
        # 返回相应严重性的影响
        if event_type in impacts and severity in impacts[event_type]:
            return impacts[event_type][severity]
        else:
            return default_impacts
    
    def _suggest_actions(self, event_type, severity):
        """建议应对气候事件的行动"""
        actions = {
            "flood": {
                "high": ["立即启动应急响应", "疏散危险区域居民", "部署救援物资", "加强堤防监控"],
                "medium": ["准备沙袋等防洪物资", "转移低洼地区贵重物品", "监控天气预报更新"],
                "low": ["清理排水系统", "准备应急包", "关注当地预警"]
            },
            "drought": {
                "high": ["实施严格水资源配给", "启动抗旱应急预案", "协助农民灌溉", "监测饮用水安全"],
                "medium": ["推广节水措施", "调整农业种植计划", "增加水资源监测频率"],
                "low": ["提高节水意识", "检查水管泄漏", "收集雨水"]
            },
            "wildfire": {
                "high": ["疏散受威胁地区", "部署消防资源", "建立防火隔离带", "发布空气质量预警"],
                "medium": ["限制户外活动", "避免明火", "准备个人防护设备"],
                "low": ["保持警惕", "遵循当地防火规定", "清理周边易燃物"]
            }
        }
        
        # 默认行动
        default_actions = ["持续监测情况", "关注官方预警", "准备必要物资"]
        
        # 返回相应严重性的行动建议
        if event_type in actions and severity in actions[event_type]:
            return actions[event_type][severity]
        else:
            return default_actions
    
    def _estimate_water_extent(self, satellite_data):
        """估计水体覆盖范围"""
        # 这是一个简化示例，实际应用需要更复杂的算法
        try:
            import numpy as np
            
            # 模拟水体覆盖估计
            # 实际应用中，应该从卫星图像分析中获取
            if isinstance(satellite_data, dict) and "segmentation" in satellite_data:
                # 假设0是水体类别
                water_pixels = np.sum(satellite_data["segmentation"] == 0)
                total_pixels = satellite_data["segmentation"].size
                water_percentage = (water_pixels / total_pixels) * 100
                
                return {
                    "percentage": water_percentage,
                    "change": 0.15 if water_percentage > 30 else 0.0,  # 模拟变化
                    "status": "normal" if 10 <= water_percentage <= 50 else "anomalous"
                }
            else:
                # 返回默认值
                return {
                    "percentage": 25.0,
                    "change": 0.0,
                    "status": "normal"
                }
                
        except Exception as e:
            print(f"估计水体覆盖错误: {e}")
            return {"error": str(e)}
    
    def _assess_vegetation_health(self, satellite_data):
        """评估植被健康状况"""
        # 这是一个简化示例，实际应用需要基于NDVI等指标
        try:
            # 模拟植被健康评估
            # 实际应用中，应该从卫星图像分析中获取NDVI或其他植被指数
            return {
                "health_score": 0.75,  # 0-1之间的健康分数
                "trend": "stable",
                "anomalies": False
            }
            
        except Exception as e:
            print(f"评估植被健康错误: {e}")
            return {"error": str(e)}

class TemporalAnalyzer:
    """时间序列数据分析器"""
    def __init__(self, config=None):
        self.config = config or {
            "anomaly_threshold": 3.0,  # Z-score阈值
            "trend_min_length": 10,    # 趋势分析的最小数据点
            "seasonality_period": 365, # 季节性周期（天）
            "heavy_precip_threshold": 50  # 大雨阈值（mm）
        }
    
    def analyze_trend(self, time_series, variable_type="temperature"):
        """分析时间序列趋势"""
        try:
            import numpy as np
            from scipy import stats
            
            # 移除NaN值
            valid_data = time_series[~np.isnan(time_series)]
            
            if len(valid_data) < self.config["trend_min_length"]:
                return {"error": "数据点不足"}
            
            # 创建时间索引
            x = np.arange(len(valid_data))
            
            # 线性回归
            slope, intercept, r_value, p_value, std_err = stats.linregress(x, valid_data)
            
            # 计算趋势方向
            direction = "increasing" if slope > 0 else "decreasing" if slope < 0 else "stable"
            
            # 判断趋势显著性
            is_significant = p_value < 0.05
            
            # 判断是否为极端值
            is_extreme = False
            if variable_type == "temperature":
                mean_val = np.mean(valid_data)
                std_val = np.std(valid_data)
                is_extreme = np.any(np.abs(valid_data - mean_val) > 3 * std_val)
            
            return {
                "slope": slope,
                "direction": direction,
                "r_squared": r_value**2,
                "is_significant": is_significant,
                "is_extreme": is_extreme,
                "confidence": min(1.0, abs(r_value) * 1.2)
            }
            
        except Exception as e:
            print(f"分析趋势错误: {e}")
            return {"error": str(e)}
    
    def analyze_precipitation(self, time_series):
        """分析降水数据"""
        try:
            import numpy as np
            
            # 移除NaN值
            valid_data = time_series[~np.isnan(time_series)]
            
            if len(valid_data) == 0:
                return {"error": "无有效数据"}
            
            # 计算统计指标
            total_precipitation = np.sum(valid_data)
            mean_daily_precipitation = np.mean(valid_data)
            max_daily_precipitation = np.max(valid_data)
            
            # 判断是否有大雨
            heavy_precipitation_events = np.sum(valid_data > self.config["heavy_precip_threshold"])
            is_heavy = heavy_precipitation_events > 0
            
            # 分析趋势
            if len(valid_data) >= 10:
                x = np.arange(len(valid_data))
                from scipy import stats
                slope, _, r_value, _, _ = stats.linregress(x, valid_data)
                trend = "increasing" if slope > 0 else "decreasing" if slope < 0 else "stable"
            else:
                trend = "insufficient_data"
                r_value = 0
            
            return {
                "total": total_precipitation,
                "mean_daily": mean_daily_precipitation,
                "max_daily": max_daily_precipitation,
                "heavy_events": heavy_precipitation_events,
                "is_heavy": is_heavy,
                "trend": trend,
                "trend_strength": abs(r_value)
            }
            
        except Exception as e:
            print(f"分析降水错误: {e}")
            return {"error": str(e)}
    
    def detect_anomalies(self, time_series):
        """检测时间序列异常"""
        try:
            import numpy as np
            from scipy import stats
            
            # 移除NaN值
            valid_data = time_series[~np.isnan(time_series)]
            
            if len(valid_data) < 5:
                return [{"warning": "数据点不足，无法检测异常"}]
            
            anomalies = []
            
            # Z-score异常检测
            z_scores = np.abs(stats.zscore(valid_data))
            anomaly_indices = np.where(z_scores > self.config["anomaly_threshold"])[0]
            
            # 如果发现异常
            if len(anomaly_indices) > 0:
                for idx in anomaly_indices:
                    anomaly_value = valid_data[idx]
                    anomalies.append({
                        "type": "z_score_anomaly",
                        "index": int(idx),
                        "value": float(anomaly_value),
                        "z_score": float(z_scores[idx]),
                        "description": f"检测到异常值: {anomaly_value:.2f}"
                    })
            
            # 趋势突变检测（简化）
            if len(valid_data) >= 20:
                mid_point = len(valid_data) // 2
                first_half = valid_data[:mid_point]
                second_half = valid_data[mid_point:]
                
                # 计算前后半段的均值差异
                mean_diff = np.abs(np.mean(first_half) - np.mean(second_half))
                overall_std = np.std(valid_data)
                
                # 如果差异超过2倍标准差，认为有突变
                if mean_diff > 2 * overall_std:
                    anomalies.append({
                        "type": "trend_change",
                        "position": "middle",
                        "magnitude": float(mean_diff),
                        "description": "检测到趋势突变"
                    })
            
            # 如果没有异常，返回正常信息
            if not anomalies:
                anomalies.append({"status": "normal", "description": "未检测到异常"})
            
            return anomalies
            
        except Exception as e:
            print(f"检测异常错误: {e}")
            return [{"error": str(e)}]

class ClimateEventAnalysisExample:
    def __init__(self):
        print("初始化气候事件分析示例")
        
        # 初始化必要的组件
        self.fusion_system = LLMMultimodalFusion()
        self.event_detection = ClimateEventDetectionSystem()
        self.satellite_analyzer = SatelliteImageAnalyzer()
        
    def prepare_sample_data(self):
        """准备示例数据"""
        print("准备示例数据...")
        
        # 生成模拟的卫星图像路径（实际应用中应使用真实图像）
        satellite_image_path = "sample_satellite_image.jpg"  # 模拟路径
        
        # 生成模拟的时间序列数据
        import numpy as np
        np.random.seed(42)  # 确保结果可复现
        
        # 生成正常温度序列，但添加一个热浪事件
        days = 365
        base_temp = 25  # 基础温度
        seasonal_variation = 10 * np.sin(np.linspace(0, 2*np.pi, days))  # 季节性变化
        random_noise = 2 * np.random.randn(days)  # 随机噪声
        
        # 创建热浪事件
        heat_wave_start = 150
        heat_wave_duration = 10
        heat_wave_amplitude = 15
        
        heat_wave = np.zeros(days)
        heat_wave[heat_wave_start:heat_wave_start+heat_wave_duration] = heat_wave_amplitude
        
        # 合成温度数据
        temperature_data = base_temp + seasonal_variation + random_noise + heat_wave
        
        # 生成降水数据（与热浪期间干燥相对应）
        base_precip = 5
        precip_noise = 10 * np.random.rand(days)
        
        # 热浪期间降水减少
        heat_wave_dry = np.ones(days)
        heat_wave_dry[heat_wave_start:heat_wave_start+heat_wave_duration] = 0.1
        
        precipitation_data = base_precip * heat_wave_dry + precip_noise
        precipitation_data = np.maximum(0, precipitation_data)  # 确保非负
        
        # 创建时间序列数据集
        time_series_data = {
            "temperature": temperature_data,
            "precipitation": precipitation_data,
            "timestamps": np.arange(days)  # 天数索引
        }
        
        # 创建文本描述
        text_description = "2025年夏季，某地区经历了一次严重的热浪事件，持续约10天，期间温度异常升高，降水稀少，多地出现高温预警。卫星图像显示该地区植被受到明显影响。"
        
        return {
            "satellite_image": satellite_image_path,
            "time_series": time_series_data,
            "text_description": text_description
        }
    
    def run_analysis_pipeline(self, sample_data):
        """运行完整的分析流程"""
        print("开始运行分析流程...")
        
        # 步骤1：提取多模态特征
        print("\n步骤1：提取多模态特征")
        
        # 由于没有真实图像，我们模拟图像分析结果
        satellite_analysis = {
            "image_description": "卫星图像显示该地区大部分为绿色植被，但在热浪区域出现明显的植被胁迫迹象，部分区域呈现黄色调。城市区域温度较高，水体面积正常。",
            "land_cover": {
                "植被": 65.2,
                "建筑区": 18.7,
                "水体": 10.5,
                "裸地": 5.6
            },
            "thermal_anomalies": [
                {"location": "城市中心", "temperature": 42.5, "severity": "高"},
                {"location": "工业区", "temperature": 40.2, "severity": "高"}
            ],
            "vegetation_stress": True,
            "stress_severity": "中度"
        }
        
        # 步骤2：分析时间序列数据
        print("\n步骤2：分析时间序列数据")
        
        # 分析温度数据
        temp_analysis = self.event_detection.temporal_analyzer.analyze_trend(
            sample_data["time_series"]["temperature"], "temperature"
        )
        
        # 分析降水数据
        precip_analysis = self.event_detection.temporal_analyzer.analyze_precipitation(
            sample_data["time_series"]["precipitation"]
        )
        
        # 检测温度异常
        temp_anomalies = self.event_detection.temporal_analyzer.detect_anomalies(
            sample_data["time_series"]["temperature"]
        )
        
        # 步骤3：提取气候指标
        print("\n步骤3：提取气候指标")
        
        # 模拟提取指标
        indicators = {
            "satellite_analysis": satellite_analysis["image_description"],
            "thermal_anomalies": satellite_analysis["thermal_anomalies"],
            "vegetation_stress": satellite_analysis["vegetation_stress"],
            "temperature_trend": temp_analysis,
            "precipitation": precip_analysis,
            "anomalies": temp_anomalies
        }
        
        # 步骤4：检测气候事件
        print("\n步骤4：检测气候事件")
        detected_events = self.event_detection.detect_events(
            indicators,
            location="模拟地区",
            timestamp="2025-07-15"
        )
        
        # 步骤5：多模态融合分析
        print("\n步骤5：多模态融合分析")
        
        # 准备LLM融合的提示
        task_description = f"分析以下多模态数据，识别气候事件并提供详细分析：\n"
        task_description += f"1. 卫星图像分析：{satellite_analysis['image_description']}\n"
        task_description += f"2. 地表覆盖：{satellite_analysis['land_cover']}\n"
        task_description += f"3. 温度趋势：{temp_analysis}\n"
        task_description += f"4. 降水分析：{precip_analysis}\n"
        task_description += f"5. 检测到的异常：{temp_anomalies}\n"
        task_description += f"6. 检测到的气候事件：{[e['event_name'] for e in detected_events]}\n"
        task_description += f"7. 文本描述：{sample_data['text_description']}\n"
        task_description += "请提供综合分析，包括事件成因、影响评估和建议措施。"
        
        # 运行融合分析（由于没有真实LLM，我们模拟结果）
        fusion_result = self._simulate_fusion_analysis(task_description, detected_events)
        
        # 步骤6：生成最终报告
        print("\n步骤6：生成最终报告")
        final_report = self._generate_final_report(
            detected_events, fusion_result, indicators
        )
        
        return final_report
    
    def _simulate_fusion_analysis(self, task_description, detected_events):
        """模拟融合分析结果"""
        # 这是一个模拟，实际应用中应使用真实的LLM
        
        # 提取主要事件
        primary_event = detected_events[0] if detected_events else None
        
        # 模拟分析结果
        analysis = {
            "summary": "基于多模态数据分析，我们识别到一次典型的极端高温事件（热浪）。",
            "findings": [
                "卫星图像显示该地区植被出现中度胁迫迹象，城市热岛效应明显。",
                "温度数据显示连续10天的异常高温，最高温度超过42°C。",
                "降水数据表明热浪期间降水量显著减少，加剧了干旱条件。",
                "时间序列分析检测到明显的温度异常和趋势变化。"
            ],
            "causes": [
                "大气环流异常导致的持续高压系统",
                "全球气候变化背景下的极端天气事件频率增加",
                "城市热岛效应加剧了城市区域的高温"
            ],
            "impacts": [
                "对人类健康的威胁，特别是老年人和慢性病患者",
                "农作物产量可能受到影响，尤其是对高温敏感的作物",
                "能源需求增加，可能导致电力供应紧张",
                "水资源短缺加剧，特别是农业灌溉需求增加"
            ],
            "recommendations": [
                "启动高温应急预案，设立避暑中心",
                "建议居民减少户外活动，特别是在午后高温时段",
                "加强电网监测，确保电力供应稳定",
                "实施节水措施，合理分配水资源",
                "加强对农作物的灌溉和保护措施"
            ]
        }
        
        # 如果有检测到的事件，更新分析
        if primary_event:
            analysis["detected_event"] = primary_event["event_name"]
            analysis["confidence"] = primary_event["confidence"]
            analysis["severity"] = primary_event["severity"]
        
        return analysis
    
    def _generate_final_report(self, detected_events, fusion_analysis, indicators):
        """生成最终分析报告"""
        report = {
            "title": "极端气候事件多模态分析报告",
            "date": "2025-07-20",
            "location": "模拟地区",
            "executive_summary": fusion_analysis["summary"],
            "detected_events": detected_events,
            "key_findings": fusion_analysis["findings"],
            "probable_causes": fusion_analysis["causes"],
            "potential_impacts": fusion_analysis["impacts"],
            "recommendations": fusion_analysis["recommendations"],
            "data_sources": [
                "卫星遥感图像",
                "地面气象站温度数据",
                "降水监测数据",
                "文本报告"
            ],
            "confidence_assessment": {
                "overall_confidence": "高",
                "event_detection": detected_events[0]["confidence"] if detected_events else "未知",
                "data_quality": "良好",
                "limitations": [
                    "模拟数据用于演示",
                    "缺乏历史数据进行比较",
                    "没有真实的卫星图像分析"
                ]
            }
        }
        
        # 格式化报告为可读文本
        formatted_report = self._format_report(report)
        
        return {
            "structured": report,
            "formatted_text": formatted_report
        }
    
    def _format_report(self, report):
        """格式化报告为可读文本"""
        lines = []
        
        # 标题
        lines.append("=" * 80)
        lines.append(report["title"].center(80))
        lines.append("=" * 80)
        
        # 基本信息
        lines.append(f"报告日期: {report['date']}")
        lines.append(f"分析地点: {report['location']}")
        lines.append("")
        
        # 执行摘要
        lines.append("## 执行摘要")
        lines.append(report["executive_summary"])
        lines.append("")
        
        # 检测到的事件
        lines.append("## 检测到的气候事件")
        for event in report["detected_events"]:
            lines.append(f"- {event['event_name']} (严重性: {event['severity']}, 置信度: {event['confidence']:.2f})")
        lines.append("")
        
        # 关键发现
        lines.append("## 关键发现")
        for finding in report["key_findings"]:
            lines.append(f"- {finding}")
        lines.append("")
        
        # 可能原因
        lines.append("## 可能原因")
        for cause in report["probable_causes"]:
            lines.append(f"- {cause}")
        lines.append("")
        
        # 潜在影响
        lines.append("## 潜在影响")
        for impact in report["potential_impacts"]:
            lines.append(f"- {impact}")
        lines.append("")
        
        # 建议措施
        lines.append("## 建议措施")
        for rec in report["recommendations"]:
            lines.append(f"- {rec}")
        lines.append("")
        
        # 置信度评估
        lines.append("## 置信度评估")
        lines.append(f"- 总体置信度: {report['confidence_assessment']['overall_confidence']}")
        lines.append(f"- 事件检测置信度: {report['confidence_assessment']['event_detection']}")
        lines.append(f"- 数据质量: {report['confidence_assessment']['data_quality']}")
        lines.append("")
        
        # 局限性
        lines.append("## 分析局限性")
        for limitation in report["confidence_assessment"]["limitations"]:
            lines.append(f"- {limitation}")
        lines.append("")
        
        # 数据来源
        lines.append("## 数据来源")
        for source in report["data_sources"]:
            lines.append(f"- {source}")
        lines.append("")
        
        lines.append("=" * 80)
        
        return "\n".join(lines)
    
    def demonstrate_analysis(self):
        """演示完整分析过程"""
        print("开始气候事件多模态分析演示...")
        
        # 准备数据
        sample_data = self.prepare_sample_data()
        
        # 运行分析
        final_report = self.run_analysis_pipeline(sample_data)
        
        # 打印报告
        print("\n===== 最终分析报告 =====\n")
        print(final_report["formatted_text"])
        
        print("\n分析演示完成。")
        return final_report

# 运行分析示例
if __name__ == "__main__":
    # 由于是模拟环境，我们可能没有真实的卫星图像和完整的LLM环境
    # 因此，这个示例主要展示了整个分析流程的框架和逻辑
    print("气候事件多模态分析示例")
    print("注意：在实际环境中运行需要真实的卫星图像和适当的计算资源")
    
    # 我们不实际执行分析，只是展示代码框架
    print("\n代码框架已展示，包含以下主要组件：")
    print("1. LLMMultimodalFusion - 多模态融合框架")
    print("2. CrossModalAttention - 跨模态注意力机制")
    print("3. SpatioTemporalModel - 时空关联建模")
    print("4. SatelliteImageAnalyzer - 卫星图像分析器")
    print("5. ClimateEventDetectionSystem - 气候事件检测系统")
    print("6. ClimateEventAnalysisExample - 完整分析流程示例")