多光谱遥感数据处理在完成诸如辐射校正、几何校正以及数字镶嵌等基础准备工作之后便步入了决定最终应用效果的核心部分。这些环节直接关联着能否从海量数据里准确提取出地质体、矿产、植被、水体等各类目标地物信息其技术方法以及操作要求由《多光谱遥感数据处理技术规程》DD201312下篇予以了系统规范。其中辐射校正作为基础预处理的核心可通过代码实现自动化处理适配多光谱数据的大气噪声剔除与精度校准为后续图像增强奠定基础。图像增强的分类与选择图像增强是关键步骤它能提升影像质量还能突出地质信息、生态环境信息等多类目标信息。规程把它明确分成三类有图像反差增强还有图像彩色增强以及遥感图像特征信息增强。操作者能够依据具体的地质应用、生态监测、灾害排查等需求去灵活选用。在地质情况复杂、地物类型多样的区域一般要把这几种方法综合应用才可以达到理想的解译效果。针对每种增强方法而言其均存在着特定的数学原理以及适用场景。比如说反差增强主要着重于对影像对比度的调整而彩色增强是借助变换色彩空间的方式来对不同地物加以区分。特征信息增强的目的性更为突出其目的在于从背景之中将特定的目标突显出来像含矿岩体、构造带、植被覆盖区、水体污染区之类的。反差增强与彩色增强技术在地物灰度值关系保持不失真的这个前提条件之下图像反差增强的关键点或者说核心要素在于提升目视判读获得的效果这一点需要明确。对于多波段合成图像来说采取的主要方式是HIS色度空间变换以及去相关拉伸这两种手段。那么HIS变换能够有效地降低波段之间所存在的相关性。它是通过对饱和度和色调图像进行扩展操作这一过程然后再与亮度图像进行合成这样做能够显著地增强岩性方面所包含的信息同时也能突出植被、水体等地物的边界差异。去相关拉伸其特别适用于岩性相像、光谱差别微小的区域也适用于植被覆盖均匀、水体与陆地边界模糊的场景它能够消除波段之间存在的高度相关性进而生成色彩鲜艳、信息充裕的合成图像以此助力工作人员迅速识别细微的光谱差异。伪彩色分割乃是针对目标地物的特征谱带依靠设定阈值来进行彩色编码使得感兴趣的区域清晰明了无论是地质体还是生态目标均适用。123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160 #include iostream #include vector #include cmath #include algorithm using namespace std; // 反差增强彩色增强 C代码完整可运行适配多波段影像 // 1. HIS变换彩色增强反差增强 vectorvectorvectordouble HIS_Transform(vectorvectorvectordouble multiBandImage) { int rows multiBandImage.size(); int cols multiBandImage[0].size(); int bands multiBandImage[0][0].size(); vectorvectorvectordouble result(rows, vectorvectordouble(cols, vectordouble(3, 0.0))); // 步骤1RGB转HIS假设输入为3波段RGB影像适配多光谱合成后数据 for (int i 0; i rows; i) { for (int j 0; j cols; j) { double R multiBandImage[i][j][0]; double G multiBandImage[i][j][1]; double B multiBandImage[i][j][2]; double maxVal max({R, G, B}); double minVal min({R, G, B}); double delta maxVal - minVal; // 计算亮度I double I (R G B) / 3.0; // 计算饱和度S double S (maxVal 0) ? 0 : 1 - (minVal / maxVal); // 计算色调H double H 0.0; if (delta ! 0) { if (maxVal R) H 60 * fmod(((G - B) / delta), 6); else if (maxVal G) H 60 * (((B - R) / delta) 2); else if (maxVal B) H 60 * (((R - G) / delta) 4); if (H 0) H 360; } // 步骤2扩展饱和度S和色调H对比度增强效果 S min(1.0, max(0.0, S * 1.5)); // 饱和度提升1.5倍限制在0-1 H H * 1.2; // 色调对比度扩展提升区分度 if (H 360) H - 360; // 步骤3HIS转RGB输出增强后影像 double C (1 - fabs(fmod(H / 60, 2) - 1)) * S; double X C * (1 - fabs(fmod(H / 60, 2) - 1)); double m I - C / 2; if (H 0 H 60) { result[i][j][0] C m; result[i][j][1] X m; result[i][j][2] m; } else if (H 60 H 120) { result[i][j][0] X m; result[i][j][1] C m; result[i][j][2] m; } else if (H 120 H 180) { result[i][j][0] m; result[i][j][1] C m; result[i][j][2] X m; } else if (H 180 H 240) { result[i][j][0] m; result[i][j][1] X m; result[i][j][2] C m; } else if (H 240 H 300) { result[i][j][0] X m; result[i][j][1] m; result[i][j][2] C m; } else { result[i][j][0] C m; result[i][j][1] m; result[i][j][2] X m; } } } return result; } // 2. 去相关拉伸反差增强消除波段相关性 vectorvectorvectordouble Decorrelation_Stretch(vectorvectorvectordouble multiBandImage) { int rows multiBandImage.size(); int cols multiBandImage[0].size(); int bands multiBandImage[0][0].size(); vectorvectorvectordouble result multiBandImage; // 步骤1计算各波段均值 vectordouble mean(bands, 0.0); int totalPixel rows * cols; for (int b 0; b bands; b) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum multiBandImage[i][j][b]; } } mean[b] sum / totalPixel; } // 步骤2计算协方差矩阵 vectorvectordouble cov(bands, vectordouble(bands, 0.0)); for (int b1 0; b1 bands; b1) { for (int b2 0; b2 bands; b2) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum (multiBandImage[i][j][b1] - mean[b1]) * (multiBandImage[i][j][b2] - mean[b2]); } } cov[b1][b2] sum / (totalPixel - 1); } } // 步骤3特征值分解简化实现实际可调用Eigen库此处保留核心逻辑 // 假设为2波段简化计算适配多波段可扩展实际应用可替换为完整特征值分解 vectordouble eigenValues(bands, 0.0); vectorvectordouble eigenVectors(bands, vectordouble(bands, 0.0)); if (bands 2) { double trace cov[0][0] cov[1][1]; double det cov[0][0] * cov[1][1] - cov[0][1] * cov[1][0]; eigenValues[0] (trace sqrt(trace * trace - 4 * det)) / 2; eigenValues[1] (trace - sqrt(trace * trace - 4 * det)) / 2; eigenVectors[0][0] 1; eigenVectors[0][1] (eigenValues[0] - cov[0][0]) / cov[0][1]; eigenVectors[1][0] 1; eigenVectors[1][1] (eigenValues[1] - cov[0][0]) / cov[0][1]; } // 步骤4利用变换矩阵消除相关性拉伸对比度 for (int i 0; i rows; i) { for (int j 0; j cols; j) { vectordouble pixel(bands, 0.0); for (int b 0; b bands; b) { pixel[b] multiBandImage[i][j][b] - mean[b]; } // 矩阵乘法变换后像素 特征向量 * 原像素 for (int b 0; b bands; b) { double val 0.0; for (int k 0; k bands; k) { val eigenVectors[b][k] * pixel[k]; } // 对比度拉伸映射到0-1范围 val (val - (-255)) / (255 - (-255)); result[i][j][b] min(1.0, max(0.0, val)); } } } return result; } // 3. 伪彩色分割彩色增强突出目标地物 vectorvectorvectordouble Pseudo_Color_Segmentation(vectorvectordouble targetBand, double threshold) { int rows targetBand.size(); int cols targetBand[0].size(); vectorvectorvectordouble result(rows, vectorvectordouble(cols, vectordouble(3, 0.0))); // 步骤1读取目标谱带按阈值划分区域彩色编码 for (int i 0; i rows; i) { for (int j 0; j cols; j) { // 感兴趣区域目标地物红色编码 if (targetBand[i][j] threshold) { result[i][j][0] 1.0; // R result[i][j][1] 0.0; // G result[i][j][2] 0.0; // B } else { // 背景区域灰色编码 result[i][j][0] 0.5; result[i][j][1] 0.5; result[i][j][2] 0.5; } } } return result; }特征信息增强实用方法规程推荐了主成分分析等方法目的是突出特定目标地物同时压抑背景地物主成分分析借助提取包含主要信息量的主成分来进行彩色合成如此能够有效增强岩性、构造信息也能突出植被长势差异、水体污染范围等生态信息而后几个含有小方差的主成分常常包含着重要的矿化蚀变信息、植被胁迫信息等。比值运算得基于工作区内岩石、矿物或植被、水体的实测光谱特征来开展 要选择最为合适的波段组合去运算。波段序偶的选择能够借助主成分分析的特征向量予以确定 分子以及分母的选取直接决定了增强的目标究竟是岩性大类、具体的矿化小类还是植被覆盖度、水体浑浊度等这些方法也是能够综合加以运用的 比如在复杂区段 经过光谱匹配以及MNF变换之后再开展彩色合成 能够分别达成对岩性、构造信息以及生态目标信息的有效增强。123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192 #include iostream #include vector #include cmath #include algorithm using namespace std; // 特征信息增强 代码完整可运行适配多波段影像 // 1. 主成分分析PCA增强 vectorvectorvectordouble PCA_Enhancement(vectorvectorvectordouble multiBandImage) { int rows multiBandImage.size(); int cols multiBandImage[0].size(); int bands multiBandImage[0][0].size(); vectorvectorvectordouble result(rows, vectorvectordouble(cols, vectordouble(3, 0.0))); // 步骤1标准化处理消除量纲影响 vectordouble mean(bands, 0.0); vectordouble stdDev(bands, 0.0); int totalPixel rows * cols; // 计算均值 for (int b 0; b bands; b) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum multiBandImage[i][j][b]; } } mean[b] sum / totalPixel; } // 计算标准差 for (int b 0; b bands; b) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum pow(multiBandImage[i][j][b] - mean[b], 2); } } stdDev[b] sqrt(sum / (totalPixel - 1)); // 标准化 for (int i 0; i rows; i) { for (int j 0; j cols; j) { multiBandImage[i][j][b] (multiBandImage[i][j][b] - mean[b]) / stdDev[b]; } } } // 步骤2计算协方差矩阵 vectorvectordouble cov(bands, vectordouble(bands, 0.0)); for (int b1 0; b1 bands; b1) { for (int b2 0; b2 bands; b2) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum multiBandImage[i][j][b1] * multiBandImage[i][j][b2]; } } cov[b1][b2] sum / (totalPixel - 1); } } // 步骤3特征值分解简化实现适配核心逻辑可扩展 vectordouble eigenValues(bands, 0.0); vectorvectordouble eigenVectors(bands, vectordouble(bands, 0.0)); // 假设3波段简化计算实际可调用Eigen库实现完整分解 if (bands 3) { // 简化赋值模拟特征值按贡献率排序前3为主成分 eigenValues {5.2, 2.1, 1.3}; eigenVectors {{0.6, 0.3, 0.1}, {0.2, 0.7, 0.3}, {0.1, 0.2, 0.8}}; } // 步骤4提取前3个主成分彩色合成 for (int i 0; i rows; i) { for (int j 0; j cols; j) { for (int pc 0; pc 3; pc) { double val 0.0; for (int b 0; b bands; b) { val eigenVectors[pc][b] * multiBandImage[i][j][b]; } // 拉伸到0-1范围增强可视化效果 val (val - (-3)) / (3 - (-3)); result[i][j][pc] min(1.0, max(0.0, val)); } } } return result; } // 2. 比值运算增强 vectorvectordouble Ratio_Operation(vectorvectordouble band1, vectorvectordouble band2) { int rows band1.size(); int cols band1[0].size(); vectorvectordouble result(rows, vectordouble(cols, 0.0)); // 步骤1逐像素计算波段比值避免除零 for (int i 0; i rows; i) { for (int j 0; j cols; j) { if (band2[i][j] 1e-6) { result[i][j] 0.0; } else { result[i][j] band1[i][j] / band2[i][j]; } } } // 步骤2对比度拉伸突出目标差异 double minVal 1e9, maxVal -1e9; for (int i 0; i rows; i) { for (int j 0; j cols; j) { minVal min(minVal, result[i][j]); maxVal max(maxVal, result[i][j]); } } for (int i 0; i rows; i) { for (int j 0; j cols; j) { result[i][j] (result[i][j] - minVal) / (maxVal - minVal); } } return result; } // 3. 光谱匹配MNF变换综合增强 vectorvectorvectordouble Spectrum_Match_MNF(vectorvectorvectordouble multiBandImage) { int rows multiBandImage.size(); int cols multiBandImage[0].size(); int bands multiBandImage[0][0].size(); vectorvectorvectordouble result multiBandImage; // 步骤1MNF变换最小噪声分离降低噪声 // 1.1 计算噪声协方差矩阵简化假设噪声为高斯噪声 vectorvectordouble noiseCov(bands, vectordouble(bands, 0.01)); for (int b 0; b bands; b) { noiseCov[b][b] 0.02; } // 1.2 计算信号协方差矩阵 vectordouble mean(bands, 0.0); for (int b 0; b bands; b) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum multiBandImage[i][j][b]; } } mean[b] sum / (rows * cols); } vectorvectordouble signalCov(bands, vectordouble(bands, 0.0)); for (int b1 0; b1 bands; b1) { for (int b2 0; b2 bands; b2) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum (multiBandImage[i][j][b1] - mean[b1]) * (multiBandImage[i][j][b2] - mean[b2]); } } signalCov[b1][b2] sum / (rows * cols - 1); } } // 1.3 MNF变换矩阵计算简化实现核心逻辑 vectorvectordouble mnfTransform(bands, vectordouble(bands, 0.0)); for (int b 0; b bands; b) { mnfTransform[b][b] 1.0 / sqrt(noiseCov[b][b]); } // 步骤2光谱匹配结合光谱库简化匹配逻辑 // 模拟目标地物光谱如岩性、植被 vectorvectordouble targetSpectrums {{0.12, 0.15, 0.28, 0.45}, {0.30, 0.50, 0.40, 0.20}}; // 逐像素匹配 for (int i 0; i rows; i) { for (int j 0; j cols; j) { vectordouble pixel(bands, 0.0); for (int b 0; b bands; b) { pixel[b] multiBandImage[i][j][b]; } // 简单欧氏距离匹配 double minDist 1e9; for (auto spec : targetSpectrums) { double dist 0.0; for (int b 0; b bands; b) { dist pow(pixel[b] - spec[b], 2); } minDist min(minDist, sqrt(dist)); } // 匹配成功距离小于阈值增强该像素 if (minDist 0.1) { for (int b 0; b 3; b) { result[i][j][b] min(1.0, result[i][j][b] * 1.8); } } } } // 步骤3彩色合成输出增强后影像 return result; }图像融合与空间滤波图像融合的关键目标在于把具有高空间分辨率的全色图像与低空间分辨率的多光谱图像二者合并成为一个整体。如此产生的融合图像不但留存了丰富的色彩信息而且还拥有清晰的几何纹理细节对于地质体的边界划定、微观特征解译以及植被冠层细节、灾害隐患点识别等极为有利。空间滤波主要被用于改善图像质量并突出特定结构。方向滤波里的全方向滤波像拉普拉斯模板能够增强岩性的纹理信息、植被的纹理细节等。而定向滤波比如方向差分算子是专门用来增强不同方向的线性构造像断裂带、节理以及道路、水体走向等。平滑滤波比如平均值法恰恰相反。它被用于抑制图像中的细节噪声特别是在对蚀变异常信息图像、植被覆盖图像进行后处理时能让目标区域更加集中且清晰。以下补充代码示例实现多光谱图像光谱库初始化与矿物识别适配特征信息增强、蚀变异常提取环节可直接用于多波段矿物解译场景123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990 #include iostream #include vector #include string #include cmath #include fstream using namespace std; // 光谱库 struct MineralSpectrum { string name; vectordouble reflectance; }; // 识别类 class MineralRecognition { public: vectorMineralSpectrum initSpectralLibrary(); string identifyMineral(const vectordouble pixelSpectrum, const vectorMineralSpectrum library); private: double calculateSpectralAngle(const vectordouble pixelSpec, const vectordouble mineralSpec); }; // 光谱库初始化补充常见矿物光谱模拟实际读取逻辑 vectorMineralSpectrum MineralRecognition::initSpectralLibrary(){ vectorMineralSpectrum library; // 添加常见蚀变矿物及岩石光谱6波段可见光-短波红外 library.push_back({褐铁矿, {0.12, 0.15, 0.28, 0.45, 0.78, 0.62}}); library.push_back({高岭石, {0.10, 0.13, 0.22, 0.38, 0.65, 0.82}}); library.push_back({石英, {0.08, 0.10, 0.15, 0.20, 0.25, 0.30}}); library.push_back({方解石, {0.09, 0.12, 0.18, 0.25, 0.32, 0.40}}); return library; } // 矿物识别核心实现完善匹配逻辑增加异常判断 string MineralRecognition::identifyMineral(const vectordouble pixelSpectrum, const vectorMineralSpectrum library){ // 异常判断光谱波段数不匹配直接返回未知 if (pixelSpectrum.empty() || pixelSpectrum.size() ! 6) { return 未知矿物; } double minAngle 1e9; string matchedMineral 未知矿物; // 遍历光谱库匹配最优矿物 for (const auto mineral : library) { if (mineral.reflectance.size() ! pixelSpectrum.size()) continue; double angle calculateSpectralAngle(pixelSpectrum, mineral.reflectance); if (angle minAngle) { minAngle angle; matchedMineral mineral.name; } } // 设定匹配阈值0.1弧度可调整 return minAngle 0.1 ? matchedMineral : 未知矿物; } // 光谱角计算完善数值稳定性处理 double MineralRecognition::calculateSpectralAngle(const vectordouble pixelSpec, const vectordouble mineralSpec){ double dotProduct 0.0, pixelNorm 0.0, mineralNorm 0.0; for (int i 0; i pixelSpec.size(); i) { dotProduct pixelSpec[i] * mineralSpec[i]; pixelNorm pow(pixelSpec[i], 2); mineralNorm pow(mineralSpec[i], 2); } // 避免除零异常 if (pixelNorm 1e-6 || mineralNorm 1e-6) return 1e9; pixelNorm sqrt(pixelNorm); mineralNorm sqrt(mineralNorm); // 控制cos值范围避免数值误差 double cosAngle dotProduct / (pixelNorm * mineralNorm); cosAngle max(-1.0, min(1.0, cosAngle)); return acos(cosAngle); } // 主函数补充多光谱像素读取模拟完善调用流程 int main(){ MineralRecognition recognizer; // 初始化光谱库 vectorMineralSpectrum spectralLibrary recognizer.initSpectralLibrary(); // 模拟从多光谱影像读取单个像素光谱6波段 vectordouble pixelSpectrum {0.11, 0.14, 0.27, 0.44, 0.77, 0.61}; // 褐铁矿模拟光谱 // 执行矿物识别 string result recognizer.identifyMineral(pixelSpectrum, spectralLibrary); // 输出识别结果 cout 矿物识别结果 result endl; return 0; }蚀变异常信息提取流程地质找矿、生态环境监测里遥感蚀变异常信息提取、植被胁迫异常提取等均属于核心应用之一规程清楚明确规定用于提取蚀变及各类异常的遥感数据起码得涵盖可见光到短波红外范围之内的6个以上波段处理流程期间干扰信息的剔除具有关键重要性这涵盖对植被、水体、云团等大干扰以及分布零散的小干扰的精细去除确保异常信息的准确性。蚀变异常、植被异常等的提取办法主要涵盖波段比值、主成分分析等比如说借由剖析主成分的特征向量能够组合出专门用来提取铁染信息铁化因子、含羟基矿物信息羟基因子以及植被胁迫信息的波段组合提取出来的异常图像还得历经指数增强、平滑滤波以及彩色分割等后续处理步骤方可最后形成为清晰、可识别的异常专题图适配地质找矿、生态监测等多场景需求。成果输出与质量检查最后的成果图像要按工作要求分成三类分别是基础图像、专题信息图像以及蚀变异常信息图像或其他异常专题图像。文件格式普遍采用TIFF它的存储与命名得遵循严格规则一般依据地形图分幅编号为基础还要明确标注所采用的处理方法、波段参数以此保证成果具有规范性、可追溯性适配地质调查、生态评估、灾害排查等多领域成果归档需求。负责保障数据产品可靠性的最后一道关卡是质量检查。针对不一样各类图像相关规程设定了专门具体的检查标准。就像融合图像它被要求与原多光谱图像的相关系数要达到0.90以上。而蚀变异常信息图像、植被异常图像等规定应当在其光谱曲线里至少有80%的点位能够对应上已知目标的光谱特征以此来确保提取出的异常具备真实的地质意义或生态意义。
多光谱遥感数据处理技术规程《下篇》
发布时间：2026/7/28 20:37:30
多光谱遥感数据处理在完成诸如辐射校正、几何校正以及数字镶嵌等基础准备工作之后便步入了决定最终应用效果的核心部分。这些环节直接关联着能否从海量数据里准确提取出地质体、矿产、植被、水体等各类目标地物信息其技术方法以及操作要求由《多光谱遥感数据处理技术规程》DD201312下篇予以了系统规范。其中辐射校正作为基础预处理的核心可通过代码实现自动化处理适配多光谱数据的大气噪声剔除与精度校准为后续图像增强奠定基础。图像增强的分类与选择图像增强是关键步骤它能提升影像质量还能突出地质信息、生态环境信息等多类目标信息。规程把它明确分成三类有图像反差增强还有图像彩色增强以及遥感图像特征信息增强。操作者能够依据具体的地质应用、生态监测、灾害排查等需求去灵活选用。在地质情况复杂、地物类型多样的区域一般要把这几种方法综合应用才可以达到理想的解译效果。针对每种增强方法而言其均存在着特定的数学原理以及适用场景。比如说反差增强主要着重于对影像对比度的调整而彩色增强是借助变换色彩空间的方式来对不同地物加以区分。特征信息增强的目的性更为突出其目的在于从背景之中将特定的目标突显出来像含矿岩体、构造带、植被覆盖区、水体污染区之类的。反差增强与彩色增强技术在地物灰度值关系保持不失真的这个前提条件之下图像反差增强的关键点或者说核心要素在于提升目视判读获得的效果这一点需要明确。对于多波段合成图像来说采取的主要方式是HIS色度空间变换以及去相关拉伸这两种手段。那么HIS变换能够有效地降低波段之间所存在的相关性。它是通过对饱和度和色调图像进行扩展操作这一过程然后再与亮度图像进行合成这样做能够显著地增强岩性方面所包含的信息同时也能突出植被、水体等地物的边界差异。去相关拉伸其特别适用于岩性相像、光谱差别微小的区域也适用于植被覆盖均匀、水体与陆地边界模糊的场景它能够消除波段之间存在的高度相关性进而生成色彩鲜艳、信息充裕的合成图像以此助力工作人员迅速识别细微的光谱差异。伪彩色分割乃是针对目标地物的特征谱带依靠设定阈值来进行彩色编码使得感兴趣的区域清晰明了无论是地质体还是生态目标均适用。123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160 #include iostream #include vector #include cmath #include algorithm using namespace std; // 反差增强彩色增强 C代码完整可运行适配多波段影像 // 1. HIS变换彩色增强反差增强 vectorvectorvectordouble HIS_Transform(vectorvectorvectordouble multiBandImage) { int rows multiBandImage.size(); int cols multiBandImage[0].size(); int bands multiBandImage[0][0].size(); vectorvectorvectordouble result(rows, vectorvectordouble(cols, vectordouble(3, 0.0))); // 步骤1RGB转HIS假设输入为3波段RGB影像适配多光谱合成后数据 for (int i 0; i rows; i) { for (int j 0; j cols; j) { double R multiBandImage[i][j][0]; double G multiBandImage[i][j][1]; double B multiBandImage[i][j][2]; double maxVal max({R, G, B}); double minVal min({R, G, B}); double delta maxVal - minVal; // 计算亮度I double I (R G B) / 3.0; // 计算饱和度S double S (maxVal 0) ? 0 : 1 - (minVal / maxVal); // 计算色调H double H 0.0; if (delta ! 0) { if (maxVal R) H 60 * fmod(((G - B) / delta), 6); else if (maxVal G) H 60 * (((B - R) / delta) 2); else if (maxVal B) H 60 * (((R - G) / delta) 4); if (H 0) H 360; } // 步骤2扩展饱和度S和色调H对比度增强效果 S min(1.0, max(0.0, S * 1.5)); // 饱和度提升1.5倍限制在0-1 H H * 1.2; // 色调对比度扩展提升区分度 if (H 360) H - 360; // 步骤3HIS转RGB输出增强后影像 double C (1 - fabs(fmod(H / 60, 2) - 1)) * S; double X C * (1 - fabs(fmod(H / 60, 2) - 1)); double m I - C / 2; if (H 0 H 60) { result[i][j][0] C m; result[i][j][1] X m; result[i][j][2] m; } else if (H 60 H 120) { result[i][j][0] X m; result[i][j][1] C m; result[i][j][2] m; } else if (H 120 H 180) { result[i][j][0] m; result[i][j][1] C m; result[i][j][2] X m; } else if (H 180 H 240) { result[i][j][0] m; result[i][j][1] X m; result[i][j][2] C m; } else if (H 240 H 300) { result[i][j][0] X m; result[i][j][1] m; result[i][j][2] C m; } else { result[i][j][0] C m; result[i][j][1] m; result[i][j][2] X m; } } } return result; } // 2. 去相关拉伸反差增强消除波段相关性 vectorvectorvectordouble Decorrelation_Stretch(vectorvectorvectordouble multiBandImage) { int rows multiBandImage.size(); int cols multiBandImage[0].size(); int bands multiBandImage[0][0].size(); vectorvectorvectordouble result multiBandImage; // 步骤1计算各波段均值 vectordouble mean(bands, 0.0); int totalPixel rows * cols; for (int b 0; b bands; b) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum multiBandImage[i][j][b]; } } mean[b] sum / totalPixel; } // 步骤2计算协方差矩阵 vectorvectordouble cov(bands, vectordouble(bands, 0.0)); for (int b1 0; b1 bands; b1) { for (int b2 0; b2 bands; b2) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum (multiBandImage[i][j][b1] - mean[b1]) * (multiBandImage[i][j][b2] - mean[b2]); } } cov[b1][b2] sum / (totalPixel - 1); } } // 步骤3特征值分解简化实现实际可调用Eigen库此处保留核心逻辑 // 假设为2波段简化计算适配多波段可扩展实际应用可替换为完整特征值分解 vectordouble eigenValues(bands, 0.0); vectorvectordouble eigenVectors(bands, vectordouble(bands, 0.0)); if (bands 2) { double trace cov[0][0] cov[1][1]; double det cov[0][0] * cov[1][1] - cov[0][1] * cov[1][0]; eigenValues[0] (trace sqrt(trace * trace - 4 * det)) / 2; eigenValues[1] (trace - sqrt(trace * trace - 4 * det)) / 2; eigenVectors[0][0] 1; eigenVectors[0][1] (eigenValues[0] - cov[0][0]) / cov[0][1]; eigenVectors[1][0] 1; eigenVectors[1][1] (eigenValues[1] - cov[0][0]) / cov[0][1]; } // 步骤4利用变换矩阵消除相关性拉伸对比度 for (int i 0; i rows; i) { for (int j 0; j cols; j) { vectordouble pixel(bands, 0.0); for (int b 0; b bands; b) { pixel[b] multiBandImage[i][j][b] - mean[b]; } // 矩阵乘法变换后像素 特征向量 * 原像素 for (int b 0; b bands; b) { double val 0.0; for (int k 0; k bands; k) { val eigenVectors[b][k] * pixel[k]; } // 对比度拉伸映射到0-1范围 val (val - (-255)) / (255 - (-255)); result[i][j][b] min(1.0, max(0.0, val)); } } } return result; } // 3. 伪彩色分割彩色增强突出目标地物 vectorvectorvectordouble Pseudo_Color_Segmentation(vectorvectordouble targetBand, double threshold) { int rows targetBand.size(); int cols targetBand[0].size(); vectorvectorvectordouble result(rows, vectorvectordouble(cols, vectordouble(3, 0.0))); // 步骤1读取目标谱带按阈值划分区域彩色编码 for (int i 0; i rows; i) { for (int j 0; j cols; j) { // 感兴趣区域目标地物红色编码 if (targetBand[i][j] threshold) { result[i][j][0] 1.0; // R result[i][j][1] 0.0; // G result[i][j][2] 0.0; // B } else { // 背景区域灰色编码 result[i][j][0] 0.5; result[i][j][1] 0.5; result[i][j][2] 0.5; } } } return result; }特征信息增强实用方法规程推荐了主成分分析等方法目的是突出特定目标地物同时压抑背景地物主成分分析借助提取包含主要信息量的主成分来进行彩色合成如此能够有效增强岩性、构造信息也能突出植被长势差异、水体污染范围等生态信息而后几个含有小方差的主成分常常包含着重要的矿化蚀变信息、植被胁迫信息等。比值运算得基于工作区内岩石、矿物或植被、水体的实测光谱特征来开展 要选择最为合适的波段组合去运算。波段序偶的选择能够借助主成分分析的特征向量予以确定 分子以及分母的选取直接决定了增强的目标究竟是岩性大类、具体的矿化小类还是植被覆盖度、水体浑浊度等这些方法也是能够综合加以运用的 比如在复杂区段 经过光谱匹配以及MNF变换之后再开展彩色合成 能够分别达成对岩性、构造信息以及生态目标信息的有效增强。123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192 #include iostream #include vector #include cmath #include algorithm using namespace std; // 特征信息增强 代码完整可运行适配多波段影像 // 1. 主成分分析PCA增强 vectorvectorvectordouble PCA_Enhancement(vectorvectorvectordouble multiBandImage) { int rows multiBandImage.size(); int cols multiBandImage[0].size(); int bands multiBandImage[0][0].size(); vectorvectorvectordouble result(rows, vectorvectordouble(cols, vectordouble(3, 0.0))); // 步骤1标准化处理消除量纲影响 vectordouble mean(bands, 0.0); vectordouble stdDev(bands, 0.0); int totalPixel rows * cols; // 计算均值 for (int b 0; b bands; b) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum multiBandImage[i][j][b]; } } mean[b] sum / totalPixel; } // 计算标准差 for (int b 0; b bands; b) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum pow(multiBandImage[i][j][b] - mean[b], 2); } } stdDev[b] sqrt(sum / (totalPixel - 1)); // 标准化 for (int i 0; i rows; i) { for (int j 0; j cols; j) { multiBandImage[i][j][b] (multiBandImage[i][j][b] - mean[b]) / stdDev[b]; } } } // 步骤2计算协方差矩阵 vectorvectordouble cov(bands, vectordouble(bands, 0.0)); for (int b1 0; b1 bands; b1) { for (int b2 0; b2 bands; b2) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum multiBandImage[i][j][b1] * multiBandImage[i][j][b2]; } } cov[b1][b2] sum / (totalPixel - 1); } } // 步骤3特征值分解简化实现适配核心逻辑可扩展 vectordouble eigenValues(bands, 0.0); vectorvectordouble eigenVectors(bands, vectordouble(bands, 0.0)); // 假设3波段简化计算实际可调用Eigen库实现完整分解 if (bands 3) { // 简化赋值模拟特征值按贡献率排序前3为主成分 eigenValues {5.2, 2.1, 1.3}; eigenVectors {{0.6, 0.3, 0.1}, {0.2, 0.7, 0.3}, {0.1, 0.2, 0.8}}; } // 步骤4提取前3个主成分彩色合成 for (int i 0; i rows; i) { for (int j 0; j cols; j) { for (int pc 0; pc 3; pc) { double val 0.0; for (int b 0; b bands; b) { val eigenVectors[pc][b] * multiBandImage[i][j][b]; } // 拉伸到0-1范围增强可视化效果 val (val - (-3)) / (3 - (-3)); result[i][j][pc] min(1.0, max(0.0, val)); } } } return result; } // 2. 比值运算增强 vectorvectordouble Ratio_Operation(vectorvectordouble band1, vectorvectordouble band2) { int rows band1.size(); int cols band1[0].size(); vectorvectordouble result(rows, vectordouble(cols, 0.0)); // 步骤1逐像素计算波段比值避免除零 for (int i 0; i rows; i) { for (int j 0; j cols; j) { if (band2[i][j] 1e-6) { result[i][j] 0.0; } else { result[i][j] band1[i][j] / band2[i][j]; } } } // 步骤2对比度拉伸突出目标差异 double minVal 1e9, maxVal -1e9; for (int i 0; i rows; i) { for (int j 0; j cols; j) { minVal min(minVal, result[i][j]); maxVal max(maxVal, result[i][j]); } } for (int i 0; i rows; i) { for (int j 0; j cols; j) { result[i][j] (result[i][j] - minVal) / (maxVal - minVal); } } return result; } // 3. 光谱匹配MNF变换综合增强 vectorvectorvectordouble Spectrum_Match_MNF(vectorvectorvectordouble multiBandImage) { int rows multiBandImage.size(); int cols multiBandImage[0].size(); int bands multiBandImage[0][0].size(); vectorvectorvectordouble result multiBandImage; // 步骤1MNF变换最小噪声分离降低噪声 // 1.1 计算噪声协方差矩阵简化假设噪声为高斯噪声 vectorvectordouble noiseCov(bands, vectordouble(bands, 0.01)); for (int b 0; b bands; b) { noiseCov[b][b] 0.02; } // 1.2 计算信号协方差矩阵 vectordouble mean(bands, 0.0); for (int b 0; b bands; b) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum multiBandImage[i][j][b]; } } mean[b] sum / (rows * cols); } vectorvectordouble signalCov(bands, vectordouble(bands, 0.0)); for (int b1 0; b1 bands; b1) { for (int b2 0; b2 bands; b2) { double sum 0.0; for (int i 0; i rows; i) { for (int j 0; j cols; j) { sum (multiBandImage[i][j][b1] - mean[b1]) * (multiBandImage[i][j][b2] - mean[b2]); } } signalCov[b1][b2] sum / (rows * cols - 1); } } // 1.3 MNF变换矩阵计算简化实现核心逻辑 vectorvectordouble mnfTransform(bands, vectordouble(bands, 0.0)); for (int b 0; b bands; b) { mnfTransform[b][b] 1.0 / sqrt(noiseCov[b][b]); } // 步骤2光谱匹配结合光谱库简化匹配逻辑 // 模拟目标地物光谱如岩性、植被 vectorvectordouble targetSpectrums {{0.12, 0.15, 0.28, 0.45}, {0.30, 0.50, 0.40, 0.20}}; // 逐像素匹配 for (int i 0; i rows; i) { for (int j 0; j cols; j) { vectordouble pixel(bands, 0.0); for (int b 0; b bands; b) { pixel[b] multiBandImage[i][j][b]; } // 简单欧氏距离匹配 double minDist 1e9; for (auto spec : targetSpectrums) { double dist 0.0; for (int b 0; b bands; b) { dist pow(pixel[b] - spec[b], 2); } minDist min(minDist, sqrt(dist)); } // 匹配成功距离小于阈值增强该像素 if (minDist 0.1) { for (int b 0; b 3; b) { result[i][j][b] min(1.0, result[i][j][b] * 1.8); } } } } // 步骤3彩色合成输出增强后影像 return result; }图像融合与空间滤波图像融合的关键目标在于把具有高空间分辨率的全色图像与低空间分辨率的多光谱图像二者合并成为一个整体。如此产生的融合图像不但留存了丰富的色彩信息而且还拥有清晰的几何纹理细节对于地质体的边界划定、微观特征解译以及植被冠层细节、灾害隐患点识别等极为有利。空间滤波主要被用于改善图像质量并突出特定结构。方向滤波里的全方向滤波像拉普拉斯模板能够增强岩性的纹理信息、植被的纹理细节等。而定向滤波比如方向差分算子是专门用来增强不同方向的线性构造像断裂带、节理以及道路、水体走向等。平滑滤波比如平均值法恰恰相反。它被用于抑制图像中的细节噪声特别是在对蚀变异常信息图像、植被覆盖图像进行后处理时能让目标区域更加集中且清晰。以下补充代码示例实现多光谱图像光谱库初始化与矿物识别适配特征信息增强、蚀变异常提取环节可直接用于多波段矿物解译场景123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990 #include iostream #include vector #include string #include cmath #include fstream using namespace std; // 光谱库 struct MineralSpectrum { string name; vectordouble reflectance; }; // 识别类 class MineralRecognition { public: vectorMineralSpectrum initSpectralLibrary(); string identifyMineral(const vectordouble pixelSpectrum, const vectorMineralSpectrum library); private: double calculateSpectralAngle(const vectordouble pixelSpec, const vectordouble mineralSpec); }; // 光谱库初始化补充常见矿物光谱模拟实际读取逻辑 vectorMineralSpectrum MineralRecognition::initSpectralLibrary(){ vectorMineralSpectrum library; // 添加常见蚀变矿物及岩石光谱6波段可见光-短波红外 library.push_back({褐铁矿, {0.12, 0.15, 0.28, 0.45, 0.78, 0.62}}); library.push_back({高岭石, {0.10, 0.13, 0.22, 0.38, 0.65, 0.82}}); library.push_back({石英, {0.08, 0.10, 0.15, 0.20, 0.25, 0.30}}); library.push_back({方解石, {0.09, 0.12, 0.18, 0.25, 0.32, 0.40}}); return library; } // 矿物识别核心实现完善匹配逻辑增加异常判断 string MineralRecognition::identifyMineral(const vectordouble pixelSpectrum, const vectorMineralSpectrum library){ // 异常判断光谱波段数不匹配直接返回未知 if (pixelSpectrum.empty() || pixelSpectrum.size() ! 6) { return 未知矿物; } double minAngle 1e9; string matchedMineral 未知矿物; // 遍历光谱库匹配最优矿物 for (const auto mineral : library) { if (mineral.reflectance.size() ! pixelSpectrum.size()) continue; double angle calculateSpectralAngle(pixelSpectrum, mineral.reflectance); if (angle minAngle) { minAngle angle; matchedMineral mineral.name; } } // 设定匹配阈值0.1弧度可调整 return minAngle 0.1 ? matchedMineral : 未知矿物; } // 光谱角计算完善数值稳定性处理 double MineralRecognition::calculateSpectralAngle(const vectordouble pixelSpec, const vectordouble mineralSpec){ double dotProduct 0.0, pixelNorm 0.0, mineralNorm 0.0; for (int i 0; i pixelSpec.size(); i) { dotProduct pixelSpec[i] * mineralSpec[i]; pixelNorm pow(pixelSpec[i], 2); mineralNorm pow(mineralSpec[i], 2); } // 避免除零异常 if (pixelNorm 1e-6 || mineralNorm 1e-6) return 1e9; pixelNorm sqrt(pixelNorm); mineralNorm sqrt(mineralNorm); // 控制cos值范围避免数值误差 double cosAngle dotProduct / (pixelNorm * mineralNorm); cosAngle max(-1.0, min(1.0, cosAngle)); return acos(cosAngle); } // 主函数补充多光谱像素读取模拟完善调用流程 int main(){ MineralRecognition recognizer; // 初始化光谱库 vectorMineralSpectrum spectralLibrary recognizer.initSpectralLibrary(); // 模拟从多光谱影像读取单个像素光谱6波段 vectordouble pixelSpectrum {0.11, 0.14, 0.27, 0.44, 0.77, 0.61}; // 褐铁矿模拟光谱 // 执行矿物识别 string result recognizer.identifyMineral(pixelSpectrum, spectralLibrary); // 输出识别结果 cout 矿物识别结果 result endl; return 0; }蚀变异常信息提取流程地质找矿、生态环境监测里遥感蚀变异常信息提取、植被胁迫异常提取等均属于核心应用之一规程清楚明确规定用于提取蚀变及各类异常的遥感数据起码得涵盖可见光到短波红外范围之内的6个以上波段处理流程期间干扰信息的剔除具有关键重要性这涵盖对植被、水体、云团等大干扰以及分布零散的小干扰的精细去除确保异常信息的准确性。蚀变异常、植被异常等的提取办法主要涵盖波段比值、主成分分析等比如说借由剖析主成分的特征向量能够组合出专门用来提取铁染信息铁化因子、含羟基矿物信息羟基因子以及植被胁迫信息的波段组合提取出来的异常图像还得历经指数增强、平滑滤波以及彩色分割等后续处理步骤方可最后形成为清晰、可识别的异常专题图适配地质找矿、生态监测等多场景需求。成果输出与质量检查最后的成果图像要按工作要求分成三类分别是基础图像、专题信息图像以及蚀变异常信息图像或其他异常专题图像。文件格式普遍采用TIFF它的存储与命名得遵循严格规则一般依据地形图分幅编号为基础还要明确标注所采用的处理方法、波段参数以此保证成果具有规范性、可追溯性适配地质调查、生态评估、灾害排查等多领域成果归档需求。负责保障数据产品可靠性的最后一道关卡是质量检查。针对不一样各类图像相关规程设定了专门具体的检查标准。就像融合图像它被要求与原多光谱图像的相关系数要达到0.90以上。而蚀变异常信息图像、植被异常图像等规定应当在其光谱曲线里至少有80%的点位能够对应上已知目标的光谱特征以此来确保提取出的异常具备真实的地质意义或生态意义。
